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_))

To clarify, by 'dissolve' do you mean 'merge' to single poly if w/in 20m & all have same properties, or 'disappear' them. And would be useful to `dput(mypolys[some_that_are_w/in_20m_and_some_not_w/in_20m, ]`, (maybe10), then copy `structure(...` to question above, between three backtics ``` above and below the `structure(...`, so everyone has data.. — Chris, Sep 30 '22 at 15:49
Hi Chris, thanks for your reply. Yes, I would like to merge them together into a single polygon. I tried to extract the data using dput in a dataset. Just added in the question. Haven't done this before so not sure if it is correct. — Lari Delazari, Sep 30 '22 at 16:23
Okay. Just a little problem with the `dput`, do the `shp_sub <- mypolys[some_that_are_w/in_20m_and_some_not_w/in_20m, ]` first, then `dput(shp_sub)`, then copy just the `structure(...` above. The above currently chokes because `shp_sub` is local to you, hence `not found` for others. And `library(sp)` above, though can be taken to `sf` so add tag, sp. — Chris, Sep 30 '22 at 17:11

score 2 · Answer 1 · answered Sep 30 '22 at 21:35

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).

Dissolving polygons by distance - R

1 Answers1

Linked