Dissolving polygons by distance - R

Question

I'm quite new to geospatial analysis in R and need help with dissolving polygons based on their location.

I have a shapefile with about 1500 polygons. Each polygon has a property code. There are a total of 4 properties in the shape attribute table. I needed to dissolve all polygons that have the same property code if they are within a distance of up to 20 meters from each other. Any insights on how can I do this?

Thanks in advance!

data subset:

> dput(SpatialPolygons(list(Polygons(list(Polygon(shp_sub)), 1))), control="all")
new("SpatialPolygons", polygons = list(new("Polygons", Polygons = list(
    new("Polygon", labpt = c(662107.09446353489, 670702.83409774397
    ), area = 15642909.541120354, hole = FALSE, ringDir = 1L, 
        coords = structure(c(662914.16988108843, 660819.24410834303, 
        661686.14174023282, 655590.40998587455, 662663.06165560661, 
        654986.60489466891, 656554.61134749465, 662744.22298060439, 
        653763.25605930295, 654452.19243265921, 662809.59039208689, 
        660147.03013036621, 662176.09769302572, 660630.43945095874, 
        664144.72288613662, 662914.16988108843, 663103.61422573915, 
        666180.53174116358, 665542.97683514666, 665529.15833193157, 
        664023.14232105017, 663416.57999173785, 664675.82179644238, 
        663715.85020686744, 663411.2359451541, 663617.32660092064, 
        663030.70651094906, 676394.80790592649, 675262.54310305626, 
        672662.00823758543, 671630.53395956859, 663103.61422573915
        ), .Dim = c(16L, 2L)))), plotOrder = 1L, labpt = c(662107.09446353489, 
670702.83409774397), ID = "1", area = 15642909.541120354)), plotOrder = 1L, 
    bbox = structure(c(653763.25605930295, 663030.70651094906, 
    664144.72288613662, 676394.80790592649), .Dim = c(2L, 2L), .Dimnames = list(
        c("x", "y"), c("min", "max"))), proj4string = new("CRS", 
        projargs = NA_character_))

Answer 1

This is an interesting problem, and I can't think of an easier solution than via iteration over rows (kinda like running cursor over a table in SQL if you get my meaning - not the first choice, but something you keep in mind if regular stuff fails).

I could not make your data work for me, so allow me to use the well known & much loved North Carolina shapefile. I have created a fake property code to have something to filter by, and a technical id to identify each polygon.

library(sf)
library(dplyr)

set.seed(1234) # or what not...

shape <- st_read(system.file("shape/nc.shp", package="sf")) %>%   # included with sf package
    mutate(id = 1:n(), # create a technical id for rows
           property = sample(LETTERS[1:4], nrow(.), replace = T)) %>% 
    select(id, property)

plot(shape["property"])

This is my start point; I will aim to dissolve polygons that

• are within 20 meters of each other
• share the same property code (color on map)

To do so I will iterate over the rows of shape data frame, first creating a vector of id's of offending polygons (meeting the two conditions = distance and property code), secondly storing their geometry in the row being evalued, and thirdly remove the offending polygons other than the row being evalued.

Note that this will decrease the number of rows in the data frame, so it is not practical to iterate via a lapply call or a for cycle. This made me resort to a while cycle, relying on overflow in the last call (when the last row has been reached the id will overflow and return NA; this will end the while cycle).

i <- shape$id[1] # start with the first id

while (!is.na(i)) {

  # get id's of offending polygons - within 20 meters and sharing property code
  offending <- st_join(select(shape[shape$id == i, ], -c(property, id)),
                       shape,
                       join = st_is_within_distance,
                       dist = 20) %>% 
    filter(property == shape$property[shape$id == i]) %>% 
    pull(id)
  
  # dissolve geometry of offending rows into i-th row
  shape$geometry[shape$id == i] <- st_union(st_geometry(shape[shape$id 
                                                              %in% offending, ])) %>% 
    st_make_valid()
  
  # eliminate offending polygons (those that offend and are not i-th row)
  shape <- shape %>% 
    filter(!id %in% setdiff(offending, shape$id[shape$id == i]))
  
  # iterate to next row - when this overflows while cycle will end
  i <- shape$id[which(shape$id == i)+1]
  
}

plot(shape["property"])

Perhaps not the most efficient piece of code, but it does the job (when other approaches fail).