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

问题描述

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

I'm looking to crop the density plot to only land while keeping to sf.

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

Here's a simple example problem:

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")

不知道我是否正在接近解决方案，但这是我一直在尝试的一些代码:

Not sure if I'm getting close to a solution but here's some of the code I've been trying:

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")

推荐答案

在美国使用AK&HI插图:

For a mask layer over the US with AK & HI inset:

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")

您可以展开边界框以摆脱绘图周围的紫色边框.

You can expand the bounding box to get rid of the purple border around the plot.

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

This does what you're asking for, but almost certainly isn't spatially accurate. It can get a point across to a general audience, but don't make any big decisions based on it.

可以在此处找到更准确的空间插值方法:

More accurate spatial interpolation methods can be found here:

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

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

这篇关于如何在ggplot中应用多边形蒙版层的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！