speeding it up: plotting an area using the extent of multiple sf objects together

Question

I have three sf objects: one is multilinestring (creeks) another multipolygon (properties) and a third polygon (coastline). These are all big files, a creek might cross several properties, several creeks might be distributed around a large area, and all are within the coastline, which is large.

creeks: Simple feature collection with 1927331 features and 25 fields Geometry type: MULTILINESTRING

properties: Simple feature collection with 2 features and 25 fields Geometry type: MULTIPOLYGON

coastline: Simple feature collection with 8567 features and 1 field Geometry type: POLYGON

I can plot like this:

ggplot() +
  geom_sf(data = creeks) +
  geom_sf(data = properties)

but I clearly have no coastline, which I would like.

I can plot them like this:

ggplot() +
  geom_sf(data = creeks) +
  geom_sf(data = properties) +
  geom_sf(data = coastline, fill = NA)

and they plot to the scale of the coastline, which is too large.

I can plot them like this:

ggplot() +
  geom_sf(data = creeks) +
  geom_sf(data = properties) +
  geom_sf(data = st_crop(coastline, st_union(creeks, properties)), fill = NA)

and I get a map of appropriate size, but it takes several minutes to plot because, I suppose, of the union operation. Also, I find the boundary is too close to the creek and property data (no buffer).

I also tried this:

pl_area <- st_buffer(st_union(creeks, properties), 10000)
ggplot() +
  geom_sf(data = creeks) +
  geom_sf(data = properties) +
  geom_sf(data = st_crop(coastline, pl_area), fill = NA)

and this works, but it takes even more minutes!

I wonder if someone has a solution where I can plot my map at the right size/resolution and doesn't take forever? My first though was I might be able to set an extent using a union of bounding boxes or something.. but if so, I don't know how. Definietly not like this:

st_buffer(
  st_union(
    st_bbox(creeks), st_bbox(properties)),
  10000)

The spatial data I have don't have the same column names etc, so rbind doesn't work. I am open to the idea that I could align their column names and then rbind, but I was hoping for neater.

Does anyone have a ideas?

Answer 1

To union bboxes, you should convert those to geometries first. But you can also construct your own bbox by combining appropriate minimum and maximum values (it's simplest form is just a named vector of 4 values). Example bellow uses shapefile from sf package and takes 2 subsets to simulate a somewhat similar case.

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(ggplot2)
library(patchwork)

nc = read_sf(system.file("shape/nc.shp", package="sf"))
subset_1 <- dplyr::filter(nc, dplyr::between(CNTY_ID, 1800, 1850 ))
subset_2 <- dplyr::filter(nc, dplyr::between(CNTY_ID, 1900, 1950 ))

(bb_1 <- st_bbox(subset_1))
#>      xmin      ymin      xmax      ymax 
#> -81.74107  35.99472 -75.77316  36.58965
(bb_2 <- st_bbox(subset_2))
#>      xmin      ymin      xmax      ymax 
#> -82.96275  35.50681 -76.69376  36.25709

# convert bboxes to sfc, then union
bb_union <- st_union(st_as_sfc(bb_1),st_as_sfc(bb_2))
st_as_text(bb_union)
#> [1] "POLYGON ((-82.96275 35.50681, -76.69376 35.50681, -76.69376 36.01401, -75.77316 35.99472, -75.77316 36.58965, -81.74107 36.58965, -81.74107 36.28277, -82.96275 36.25709, -82.96275 35.50681))"

# or combine minimums of *xmin* & *ymin* pairs and maximums of *xmax* & *ymax* pairs:
(bb_merged <- c(pmin(bb_1[1:2], bb_2[1:2]), 
                pmax(bb_1[3:4], bb_2[3:4])))
#>      xmin      ymin      xmax      ymax 
#> -82.96275  35.50681 -75.77316  36.58965


p1 <- ggplot() +
  geom_sf(data = nc) +
  geom_sf(data = subset_1, fill = "darkred") +
  geom_sf(data = subset_2, fill = "darkgreen") +
  geom_sf(data = bb_union, linewidth = 1, fill = "gold", alpha = .2) +
  theme_minimal()

p2 <- ggplot() +
  geom_sf(data = st_crop(nc, bb_merged)) +
  geom_sf(data = subset_1, fill = "darkred") +
  geom_sf(data = subset_2, fill = "darkgreen") +
  ggtitle("st_crop(nc, bb_merged))") +
  theme_minimal()
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries

p3 <- ggplot() +
  geom_sf(data = st_crop(nc, bb_union)) +
  geom_sf(data = subset_1, fill = "darkred") +
  geom_sf(data = subset_2, fill = "darkgreen") +
  ggtitle("st_crop(nc, bb_union)") +
  theme_minimal()
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
  
p1/p2/p3

You could also consider simplification to speed up plotting (sf::st_simplify() , rmapshaper::ms_simplify() ). And you may need to adjust when working with different projections / dataums, either directly or indirectly -- note how bbox defined by 4 corners can cut into geometries in this example.