Getting spatial polygon intersections, excluding borders

Question

I am trying to find which City Council districts represent which zip codes in NYC by using the sf library and leaflet. I have spatial polygons for both council districts and zip codes. My problem is that whenever I try to find the intersection between the polygons, I get extra zip codes that have overlapping borders with the district borders. There is no portion of the district that falls into those zip codes beyond the border, but using a spatial intersection includes those bordering zip codes.

Is there a spatial manipulation function that would return just the intersection of the polygons without the places where the polygons only share a border?

Here is a visual of my issue (see the overlapping border of zip code (in grey) and council district (in black):

And my code:

library(sf)

#NYC Zip Codes
nyc_zips <- read_sf('ZIP_CODE_040114.shp')
#Left joining zip code boundaries with total count of drivers in each zip
ny_zips_drivers <- left_join(nyc_zips, driver_zips, by = 'ZIPCODE')

#City Council
cc_shp <- read_sf('nycc_22a/nycc.shp')

#Spatial Intersection of council districts and zip codes
cc_intersection <- st_intersection(x = cc_shp, y = ny_zips_drivers)

Try to use a negative buffer. The issue here is that your need to be mindful of your CRS, and specifically the units. Assuming that your shop is projected (i.e. no lonlat), try `cc_intersection <- st_intersection(x = cc_shp, y = st_buffer(ny_zips_drivers, -100))`. After you can use the ids of the result for filtering out the non-buffered polygons — dieghernan, Mar 30 '22 at 20:44
This worked! Thank you so much. Didn't think to add a buffer originally but now it all works. — phillipwong, Apr 04 '22 at 16:24

Jindra Lacko · Answer 1 · 2022-03-30T21:08:58.290

I believe you will find useful the model argument of sf::st_intersection(). It works only for unprojected CRS (in the context of projected CRS GEOS will be used, which works in context of DE-9IM which does not contain this feature).

It can be used to specify whether polygons contain (model = "closed") or do not contain (model = "open") their border.

For an example consider this piece of code, built on top of the well known and much loved nc.shp that ships with {sf}:

library(sf)
library(dplyr)

shape <- st_read(system.file("shape/nc.shp", package="sf"))   # included with sf package
  
# will S2 be used? at all?? If not >> st_transform to EPSG:4326 as a good default
st_is_longlat(shape)
# [1] TRUE

#  County Mecklenburg (as in Charlotte of Mecklenburg-Strelitz)
strelitz <- shape[shape$NAME == 'Mecklenburg', ] 

# closed behaviour - including borders = co. Mecklenburg + all its NC neighbours
st_intersection(shape, strelitz, model = "closed") %>% 
  pull(NAME)
# [1] "Iredell"     "Lincoln"     "Mecklenburg" "Cabarrus"    "Gaston"      "Union"   

# open behaviour - excluding borders = co. Mecklenburg only
st_intersection(shape, strelitz, model = "open") %>% 
  pull(NAME)
# [1] "Mecklenburg"

Getting spatial polygon intersections, excluding borders

1 Answers1