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

本文介绍了如何在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( "services1.arcgis/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( "services1.arcgis/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软件包(v1.0.0)创建于2021-04-05 sup>

Created on 2021-04-05 by the reprex package (v1.0.0)

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

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:

rspatial/raster/analysis/4-interpolation.html

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