# Using sf to combine polygons that share borders?

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 -

``````library(sf)

# Sample data

# Plot 1
plot(st_geometry(nc))
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)
plot(geom, add = TRUE, col = "red")
``````

Plot 1 shows the original layer, with a specific polygon n=25 shown in red - Plot 2 shows the layer with the new geometries, where the new polygon n=25 is merged with its neighbors - • Your `pol` is not defined (as in, running your code produces an error). – CoderGuy123 Oct 2 '19 at 1:40
• The first code section is to illustrate the concept, it's not reproducible. The second code section is a specific (reproducible) example. – Michael Dorman Oct 2 '19 at 5:43

After hours of trying, I made a simple function that merges units by some grouping variable, which is usually what one wants to do.

``````> library(sf)
> library(plyr)
Reading layer `nc' from data source `/home/emil/R/x86_64-pc-linux-gnu-library/3.6/sf/shape/nc.shp' using driver `ESRI Shapefile'
Simple feature collection with 100 features and 14 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
epsg (SRID):    4267
> #what we want to merge
> plot(st_geometry(nc))
> plot(st_geometry(nc[1:2, ]), add = TRUE, col = "red")
> plot(st_geometry(nc[-(1:2), ]), add = TRUE, col = "green")
`````` ``````> #merge units as desired by group variable
> st_union_by = function(geo, group) {
+   # browser()
+   geo
+   group
+
+   y2 = list()
+   #loop over by groups and merge units
+   for (i in unique(group)) {
+     #which units
+     z = geo[group == i]
+
+     #merge
+     y = Reduce(st_union, z)
+     y2[[i]] = y
+   }
+
+   st_sfc(y2)
+ }
> nc\$x = c(rep(1, 2), rep(2, 98))
> nc\$x
 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
> #merge
> nc2_geo = st_union_by(nc\$geometry, nc\$x)
> #plot results
> plot(nc2_geo)
> plot(st_geometry(nc2_geo), add = TRUE, col = "red")
> plot(st_geometry(nc2_geo), add = TRUE, col = "green")
`````` ``````> #make a new sf df
> nc2 = plyr::ddply(nc, "x", function(d) {
+   data.frame(
+     AREA = sum(d\$AREA)
+   )
+ })
> #make sf df like originally
> st_geometry(nc2) = nc2_geo
> nc2
Simple feature collection with 2 features and 2 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
epsg (SRID):    NA
proj4string:    NA
x   AREA                       geometry
1 1  0.175 POLYGON ((-81.47276 36.2343...
2 2 12.451 MULTIPOLYGON (((-81.76536 3...
> plot(nc2)
`````` 