How to apply a polygon mask layer in ggplot

user_4828653 0
Viewed 591

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

sf ggplot2 r
2 Answers

Create a rectangle of the same plot dimensions:

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

Create an outline of the US and extract only the lat/lon points into a dataframe:

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

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

Bind the two together to create a polygon with a hole in the middle:

mask <- rbind(rec_box,outline)

Add a geom_polygon() to plot the mask data and color appropriately:

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

Everything combined:

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

Really a thing of beauty.

user_48286530

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

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:

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

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

user_75473270
Related