Nearest adjacent polygon, distance and closest point to other polygons

Question

I have multiple polygons in a dataset and I would like to:

• Identify the nearest polygon to each polygon and what the distance between them is
• Calculate the coordinates of where the nearest parts of the two polygons are (so I can draw a line and visually check the distances)
• If the distance is 800 metres of less, join the polygons together to make multipart polygons

This code does half of my first ask and I know st_distance can do the latter. I was hoping for a solution that wouldn't need for a matrix of every distance between every polygon to be generated.

library(sf)
library(dplyr)

download.file("https://drive.google.com/uc?export=download&id=1-I4F2NYvFWkNqy7ASFNxnyrwr_wT0lGF" , destfile="ProximityAreas.zip")
unzip("ProximityAreas.zip")
Proximity_Areas <- st_read("Proximity_Areas.gpkg") 

Nearest_UID <- st_nearest_feature(Proximity_Areas)

Proximity_Areas <- Proximity_Areas %>% 
        select(UID) %>% 
        mutate(NearUID = UID[Nearest_UID])

Is there a method of producing two outputs 1) an appended Proximity_Areas file that included the distance and XY coorindates for the nearest points for the UID and Neatest_UID and 2) a file that looks similar to the original Proximity_Areas file, just with merged polygons if the criteria is met?

Answer 1

Once you have created index of nearest neighbors you can calculate the connecting lines via a sf::st_nearest_points() call.

An interesting aspect is that if you make the call on geometries (not sf, but sfc objects) you do the calculation pairwise (i.e. not in a matrix way).

The call will return linestrings, which is very helpful since you can calculate their length and have two of your objectives (nearest points & distance) at a single call...

lines <- Proximity_Areas %>% 
  st_geometry() %>% # extact geometry
  # create a line to nearest neighbour as geometry
  st_nearest_points(st_geometry(Proximity_Areas)[Nearest_UID], pairwise =T) %>%
  # make sf again (so it can hold data)
  st_as_sf() %>% 
  # add some data - start, finish, lenght
  mutate(start = Proximity_Areas$UID,
         end = Proximity_Areas$UID[Nearest_UID],
         distance = st_length(.)) 

glimpse(lines)  
# Rows: 39
# Columns: 4
# $ x        <LINESTRING [m]> LINESTRING (273421.5 170677..., LINESTRING (265535.1 166136..., LINESTRING (265363.3 1…
# $ start    <chr> "U001", "U002", "U003", "U004", "U005", "U006", "U007", "U008", "U009", "U010", "U011", "U012", "…
# $ end      <chr> "U026", "U010", "U013", "U033", "U032", "U014", "U028", "U036", "U011", "U008", "U028", "U030", "…
# $ distance [m] 317.84698 [m], 579.85131 [m], 529.67907 [m], 559.96441 [m], 0.00000 [m], 80.54011 [m], 754.94311 [m…

mapview::mapview(lines)

The part about joining close objects together is a bit tricky, since you don't know how many polygons you will end up with - you can have a polygon A that is far from C, but will end up merged since both are close to B. This does not vectorize easily and you are likely to end up running a while loop. For a possible approach consider this related answer Dissolving polygons by distance - R