如何在ggplot中应用多边形蒙版层

Question

我希望将密度图裁剪为仅着陆，同时保持sf。

这是一个简单的示例问题：

library(tidyverse)
library(sf)
library(albersusa)
library(ggthemes)
library(jsonlite)

dat <-
  fromJSON(
    "https://services1.arcgis.com/Hp6G80Pky0om7QvQ/arcgis/rest/services/Fortune_500_Corporate_Headquarters/FeatureServer/0/query?where=1%3D1&outFields=LATITUDE,LONGITUDE,NAME,PROFIT&outSR=4326&f=json"
  )

dat <- as.data.frame(dat$features$attributes)

top_50 <- dat %>%
  arrange(desc(PROFIT)) %>%
  head(50)

ggplot() +
  geom_sf(data = usa_sf()) +
  geom_density_2d_filled(aes(x = LONGITUDE, y = LATITUDE),
                         data = top_50,
                         alpha = .5) +
  xlim(-125,-66.5) +
  ylim(20, 50) +
  theme_map() +
  theme(legend.position = "none")

不确定我是否接近解决方案，但这是我一直在尝试的一些代码：

test <- (MASS::kde2d(
  top_50$LONGITUDE, top_50$LATITUDE,
  lims =  c(-125,-66.5, 20, 50)
))

ggpoly2sf <- function(poly, coords = c("long", "lat"),
                      id = "group", region = "region", crs = 4326) {
  sf::st_as_sf(poly, coords = coords, crs = crs) %>% 
    group_by(!! as.name(id), !! as.name(region)) %>% 
    summarize(do_union=FALSE) %>%
    sf::st_as_sf("POLYGON") %>% 
    ungroup() %>%
    group_by(!! as.name(region)) %>%
    summarize(do_union = TRUE) %>%
    ungroup()
}

v <- contourLines(test)
vv <- v
for (i in seq_along(v)) vv[[i]]$group <- i
vv <- do.call(rbind, lapply(vv, as.data.frame))

dsi_sf <- ggpoly2sf(vv, coords = c("x", "y"), region = "level") %>% st_as_sf()

usa <- usa_sf()

dsi_i_sf <- st_intersection(usa$geometry, dsi_sf)

ggplot() + 
  geom_sf(data=usa) +
  geom_sf(data=dsi_i_sf,color="red") +
  geom_density2d_filled(aes(x = LONGITUDE, y = LATITUDE),
                 data = top_50,alpha=.4) +
  xlim(-125,-66.5) +
  ylim(20, 50) +
  theme(legend.position = "none")

Answer 1

创建一个具有相同绘图尺寸的矩形：

rec_box <- data.frame(x=c(-125,-125,-66.5,-66.5,-125), y=c(20,50,50,20,20))

创建美国的轮廓并仅将纬度/经度点提取到数据框中：

outline <- map("usa", plot=FALSE)

outline <- data.frame(x=outline$x,y=outline$y)

将两者绑定在一起，创建一个中间有洞的多边形：

mask <- rbind(rec_box,outline)

添加geom_polygon() 以适当地绘制掩码数据和颜色：

  geom_polygon(data=mask,
                aes(x=x,y=y),color="white",fill="white")

一切相结合：

library(tidyverse)
library(sf)
library(albersusa)
library(ggthemes)
library(jsonlite)

dat <-
  fromJSON(
    "https://services1.arcgis.com/Hp6G80Pky0om7QvQ/arcgis/rest/services/Fortune_500_Corporate_Headquarters/FeatureServer/0/query?where=1%3D1&outFields=LATITUDE,LONGITUDE,NAME,PROFIT&outSR=4326&f=json"
  )

dat <- as.data.frame(dat$features$attributes)

top_50 <- dat %>%
  arrange(desc(PROFIT)) %>%
  head(50)

usa <- usa_sf()

outline <- map("usa", plot=FALSE)

outline <- data.frame(x=outline$x,y=outline$y)

rec_box <- data.frame(x=c(-125,-125,-66.5,-66.5,-125), y=c(20,50,50,20,20))

mask <- rbind(rec_box,outline)

ggplot() + 
  geom_sf(data = usa_sf()) +
  geom_density_2d_filled(aes(x = LONGITUDE, y = LATITUDE),
                         data = top_50,
                         alpha = .5) +
  xlim(-125,-66.5) +
  ylim(20, 50) +
  geom_polygon(data=mask,
                aes(x=x,y=y),color="white",fill="white") +
  theme_map() +
  theme(legend.position = "none")

真的很美。

Answer 2

对于美国上空插入 AK 和 HI 的遮罩层：

library(tidyverse)
library(sf)
library(albersusa) 
library(ggthemes)
library(jsonlite)
library(spatstat)


dat <-
  fromJSON(
    "https://services1.arcgis.com/Hp6G80Pky0om7QvQ/arcgis/rest/services/Fortune_500_Corporate_Headquarters/FeatureServer/0/query?where=1%3D1&outFields=LATITUDE,LONGITUDE,NAME,PROFIT&outSR=4326&f=json"
  )

dat <- as.data.frame(dat$features$attributes)

top_50 <- dat %>%
  arrange(desc(PROFIT)) %>%
  head(50)

usa <- usa_sf()

top50sf <- st_as_sf(top_50, coords = c("LONGITUDE", "LATITUDE")) %>%
  st_set_crs(4326) %>%
  st_transform(st_crs(usa))

# usa polygons combined
usa_for_mask <- usa_sf() %>%
  st_geometry() %>%
  st_cast('POLYGON') %>%
  st_union()

# bounding box of us & inset AK + HI,
# expand as needed 
us_bbox <- st_bbox(usa_for_mask) %>%
             st_as_sfc() %>% 
             st_as_sf()

us_mask <- st_difference(us_bbox, usa_for_mask)


ggplot() +
  geom_sf(data = usa) +
  geom_density_2d_filled(aes(x = LONGITUDE, y = LATITUDE),
                         data = top_50,
                         alpha = .5) +
  geom_sf(data = us_mask, fill = 'white') +
  xlim(-125,-66.5) +
  ylim(20, 50) +
  theme_map() +
  theme(legend.position = "none")

^{由reprex package (v1.0.0) 于 2021-04-05 创建}

您可以扩展边界框以消除绘图周围的紫色边框。

这可以满足您的要求，但几乎可以肯定在空间上不准确。它可以向普通观众传达一个观点，但不要根据它做出任何重大决定。

更准确的空间插值方法可以在这里找到：

https://rspatial.org/raster/analysis/4-interpolation.html

https://mgimond.github.io/Spatial/interpolation-in-r.html