I have rows of respondents with their residential census tract polygon and am trying to generate first-order aggregated areas. Meaning I want a given respondent’s census tract polygon to be merged with all the polygons touching its borders. A separate file contains all possible US census tract polygons.

st_union seems appropriate only for aggregating together ALL interesting polygons to a single mass polygon.

Using st_intersection and then aggregating by row ID seems viable, but I’m unsure how to do the aggregation procedure.

I hope to remain within the sf environment.

Similar threads I’ve quered: Dissolve only overlapping polygons in R

Joining polygons in R


If I understood correctly, given a polygon layer with N features, what you need is to generate a new set of N geometries where each geometry is the result of a union with it's neighbors.

A straightforward solution can be to use a for loop going over the 1:N polygons, applying st_union on the current set of neighbors for polygon i in each "round", as in -

geom = list()
for(i in 1:nrow(pol)) {
  geom[[i]] = st_union(pol[pol[i, ], ])
geom = do.call(c, geom)

Where pol is the original layer and geom is the new set of merged geometries.

Here is a reproducible example using the nc dataset -


# Sample data
nc = st_read(system.file("shape/nc.shp", package="sf"))

# Plot 1
plot(st_geometry(nc[25, ]), add = TRUE, col = "red")

# Merging each polygon with its neighbors
geom = list()
for(i in 1:nrow(nc)) {
  geom[[i]] = st_union(nc[nc[i, ], ])
geom = do.call(c, geom)

# Plot 2
plot(geom[25], add = TRUE, col = "red")

Plot 1 shows the original layer, with a specific polygon n=25 shown in red -

enter image description here

Plot 2 shows the layer with the new geometries, where the new polygon n=25 is merged with its neighbors -

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.