Improve handling of filter_geom inputs with "MULTIPOLYGON" geometry

[elipousson]

Is your feature request related to a problem? Please describe.

I got seriously tripped up today by the handling of the MULTIPOLYGON filter_geom objects when dealing with non-overlapping polygons that can't be merged with sf::st_union(). I originally wrote the code to handle that transformation so if it gave me trouble, it may still catch another user unaware.

Describe the solution you'd like

Casting the input MULTIPOLYGON to a MULTIPOINT and then using sf::st_concave_hull seems more fool-proof option. This reprex shows the differences between the two approaches:

nc <- sf::st_read(system.file("shape/nc.shp", package="sf"))
#> Reading layer `nc' from data source 
#>   `/Users/elipousson/Library/R/arm64/4.3/library/sf/shape/nc.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27

dare_co_geom <- sf::st_geometry(nc[56, ])

plot(dare_co_geom)

dare_co_poly <- dare_co_geom |> 
  sf::st_union() |> 
  sf::st_cast(to = "POLYGON")

dare_co_poly <- dare_co_poly[[1]]

plot(dare_co_poly, border = "blue", add = TRUE)

dare_co_hull <- dare_co_geom |> 
  sf::st_cast(to = "MULTIPOINT") |> 
  sf::st_concave_hull(ratio = 1, allow_holes = FALSE)

plot(dare_co_hull, border = "orange", add = TRUE)

ashe_co_geom <- sf::st_geometry(nc[1, ])

plot(ashe_co_geom)

ashe_co_poly <- ashe_co_geom |> 
  sf::st_union() |> 
  sf::st_cast(to = "POLYGON")

ashe_co_poly <- ashe_co_poly[[1]]

plot(ashe_co_poly, border = "blue", add = TRUE)

ashe_co_hull <- ashe_co_geom |> 
  sf::st_cast(to = "MULTIPOINT") |> 
  sf::st_concave_hull(ratio = 1, allow_holes = FALSE)

plot(ashe_co_hull, border = "orange", add = TRUE)

^{Created on 2024-03-11 with reprex v2.1.0}

Describe alternatives you've considered

I also looked at using sf::st_concave_hull but it can get a little funky depending on the density of nodes in the input geometry. I also considered suggesting that the handling of filter_geom inputs could be controlled by a secondary parameter or an option, but that seems excessive considering the relatively narrow use case. If it was a "settable" option, you could have a concave hull, convex hull, bounding box, or error always option – but again that seems excessive.

I've got a working version of this already so I'm going to open a PR. Happy to discuss on here if you prefer though.

[JosiahParry]

This appears to be working based on the example that you have. Do you have a reprex where the current implementation doesn't work?

nc <- sf::st_read(system.file("shape/nc.shp", package="sf"))
#> Reading layer `nc' from data source 
#>   `/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/sf/shape/nc.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
dare_co_geom <- sf::st_geometry(nc[56, ])

library(arcgis)
#> Attaching core arcgis packages:
#> → arcgisutils v0.2.0.9000
#> → arcgislayers v0.2.0

flayer <- arc_open("https://services.arcgis.com/P3ePLMYs2RVChkJx/ArcGIS/rest/services/USA_Counties_Generalized_Boundaries/FeatureServer/0")

counties <- arc_select(flayer, filter_geom = dare_co_geom)
#> ! `filter_geom` cannot be a "MULTIPOLYGON" geometry.
#> ℹ Using `sf::st_union()` and `sf::st_cast()` to create a "POLYGON" for
#>   `filter_geom`.

plot(counties$geometry)
plot(dare_co_geom, border = "red", lty = 3, lwd = 3, add = TRUE)

^{Created on 2024-03-21 with reprex v2.0.2}

[elipousson]

It depends on the data source. Here is an example with point data and a plot to show how the existing method gets only part of the data. I'm still not certain using a hull is the best or only way to fix this but I think it is an improvement.

library(sf)
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(arcgislayers)
library(ggplot2)

nc <- read_sf(system.file("shape/nc.shp", package="sf"))

dare_co_geom <- sf::st_geometry(nc[56, ])

dare_co_basemap <- ggplot() +
  geom_sf(data = dare_co_geom, fill = "lightgreen") +
  theme_void()

url <- "https://services.arcgis.com/NuWFvHYDMVmmxMeM/ArcGIS/rest/services/NCPedBikeCrashes/FeatureServer/1"

dare_co_data <- arcgislayers::arc_read(url, filter_geom = dare_co_geom)
#> ! `filter_geom` cannot be a "MULTIPOLYGON" geometry.
#> ℹ Using `sf::st_union()` and `sf::st_cast()` to create a "POLYGON" for
#>   `filter_geom`.

dare_co_poly <- dare_co_geom |> 
  sf::st_union() |> 
  sf::st_cast("POLYGON")

dare_co_poly <- dare_co_poly[1]

dare_co_hull <- dare_co_geom |> 
  sf::st_cast(to = "MULTIPOINT") |> 
  sf::st_concave_hull(ratio = 1, allow_holes = FALSE)

dare_co_hull_data <- arcgislayers::arc_read(url, filter_geom = dare_co_hull)

dare_co_basemap +
  geom_sf(data = dare_co_hull, fill = NA, linewidth = 1.25, color = "red") +
  geom_sf(data = dare_co_hull_data, color = "darkred", size = 1.5, alpha = 0.5) +
  geom_sf(data = dare_co_poly, fill = NA, linewidth = 1, color = "blue", alpha = 0.5) +
  geom_sf(data = dare_co_data, color = "darkblue")

^{Created on 2024-03-21 with reprex v2.1.0}