Filter and exclude a level in stat_density2d

I have code and data like that to map some locations in the USA:

library(sf)

library(tidyverse)

set.seed(123) # for reproducibility
Latitude = round(runif(5000, 31, 46), digits = 5)
Longitude = round(runif(5000, -120, -80), digits = 5)
df = data.frame(Latitude, Longitude) 


world <- ne_countries(scale = "medium", returnclass = "sf")
usa = filter(world,admin =="United States of America")

# Plot
ggplot()  
  geom_sf(data = usa, fill = "blue",color = "black",alpha=.9)  
  coord_sf(
    xlim = c(-119, -74),
    ylim = c(22, 51),
   default_crs = sf::st_crs(4326),
    crs = st_crs("ESRI:102003"),
    expand = TRUE,
    lims_method = "box",
    label_axes = list(
      bottom = "E", top = "E",
      left = "N", right = "N"
    )) stat_density2d_filled(data = df, aes(x = Longitude, y = Latitude, fill=(..level..),
         alpha = (..level..)),  geom = "polygon")

I have the output like this:

I am looking to exclude some levels like (0,0002,0004], and (0,0004,0006]. Additionally, I would like to manually color the levels. I have this code scale_fill_manual(values = c("khaki1", "khaki2", "khaki2", "orange", "orange", "red", "red", "red", "red")) , however would like to use scale_fill_gradientn and color palettes in the RColorBrewer package like YlOrRd.

CodePudding user response：

I don't have the usa basemap layer, but does this do what you want with your density layer?

library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(ggplot2)
library(tidyverse)

set.seed(123) # for reproducibility
Latitude = round(runif(5000, 31, 46), digits = 5)
Longitude = round(runif(5000, -120, -80), digits = 5)
df = data.frame(Latitude, Longitude) 



breaks <- seq(from = 0.0006, to = 0.0020, by = 0.0002)

# Plot
ggplot()  
  # geom_sf(data = usa, fill = "blue",color = "black",alpha=.9)  
  coord_sf(
    xlim = c(-119, -74),
    ylim = c(22, 51),
    default_crs = sf::st_crs(4326),
    crs = st_crs("ESRI:102003"),
    expand = TRUE,
    lims_method = "box",
    label_axes = list(
      bottom = "E", top = "E",
      left = "N", right = "N"
    ))  
  stat_density2d_filled(data = df, 
                        aes(x = Longitude, y = Latitude, 
                            fill = after_stat(level), 
                            alpha = after_stat(level)),  
                        geom = "polygon",
                        breaks = breaks)  
  scale_fill_brewer(palette = "YlOrRd")

^{Created on 2023-01-31 by the reprex package (v2.0.1)}