2

I am using the sf package in R along with the arcpullr package to pull in data from an ArcGIS REST service and work with it as an sf object. I have run into an issue with a MULTIPOLYGON where sf is only buffering a part of the MULTIPOLYGON (i.e., it buffers one polygon but only tiny slivers of the other). I have not been able to replicate the problem when buffering the example found here.

Here is an MRE (sorry, you'll have to install arcpullr if you don't have it).

# install.packages("arcpullr")
library(arcpullr)
library(sf) # redundant since loaded with arcpullr, but here for brevity
library(ggplot2)

tax_parcel_url <- paste0(
  "https://mapservices.legis.wisconsin.gov/arcgis/rest/services/", 
  "WLIP/Parcels/FeatureServer/0"
)
parcel <- 
  get_spatial_layer(tax_parcel_url, where = "PARCELID = 'HA-11'") %>%
  st_transform(crs = 3071)
parcel_buffer <- st_buffer(parcel, dist = 10)

# map of parcel
ggplot(data = parcel) + geom_sf() # this is correct

# map of parcel and buffer - buffer "misses" part of multipolygon
ggplot(data = parcel_buffer) + 
  geom_sf(color = "blue") + 
  geom_sf(data = parcel, color = "red", alpha = 0.3) + 
  theme_bw()
AndrewGB
  • 16,126
  • 5
  • 18
  • 49
code_cowboy
  • 596
  • 5
  • 18
  • 1
    Hi! I just want to point out that I get a "Timeout was reached" error when running get_spatial_layer. Can you try to rerun the code? Or maybe share the object in other ways? – agila Dec 09 '21 at 15:17
  • Thanks @agila. Mine runs in 0.41 seconds. It's possible there was a server issue when you ran it. Try it again? – code_cowboy Dec 15 '21 at 02:29

1 Answers1

4

You can pipe in st_make_valid to parcel, which will correct the geometry with the sf object (i.e., make the geometry valid). This will allow for the buffer to plot correctly on all parts of the multipolygon.

library(arcpullr)
library(sf)
library(ggplot2)

parcel <-
  arcpullr::get_spatial_layer(tax_parcel_url, where = "PARCELID = 'HA-11'") %>%
  sf::st_transform(crs = 3071) %>%
  sf::st_make_valid()

parcel_buffer <- sf::st_buffer(parcel, dist = 10)

ggplot(data = parcel_buffer) +
  geom_sf(color = "blue") +
  geom_sf(data = parcel,
          color = "red",
          alpha = 0.3) +
  theme_bw()

Output

enter image description here

Mapview

You can see that the buffer is clear in mapview as well.

library(mapview)

mapview::mapview(parcel_buffer) + mapview::mapview(parcel)

enter image description here

AndrewGB
  • 16,126
  • 5
  • 18
  • 49