Spatial 1 to 1 join by proximity in sf and R?

Question

I'm trying to join a Simple Feature of POINT geometry type to another Simple Feature of POLYGON geometry type using st_join and nngeo::st_nn as the person asking this other question also did. The difference is that I want a 1:1 match from both sides (i.e. one point per polygon).

Above is a picture of how my data looks like when mapped. My code looks like:

matched<-st_join(buildings, addresses, join=nngeo::st_nn, maxdist=50,k=1)

Instead of a 1 to 1 match, I keep getting more than one point merged to the same polygon or more than one polygon merged to the same point depending on how order the arguments to st_join, with repeated points or polygons accordingly.

What I want this to do is to match each polygon to the closest point which also has that polygon as the closest to itself (i.e., to match polygons to points when they are mutually the closest to each other), and to randomise matching when more than one point or polygon meet those criteria.

Is there any way of doing that with sf or sp in r? I have thought about extracting the distance matrices and doing a two-sided matching between the two sets of features (like a Gale- Shapley algorithm), but I want to give it a shot with what is already available before putting many hours into coding that.

Answer 1

Here is an example, with 100 polygons and 150 points, using similar syntax of st_join. I'm never getting more than one match when using k=1 (there can be zero matches when the nearest feature is further than maxdist).

If you can please post a reproducible example where you are getting more than one match will be happy to look into it.

library(sf)

## Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0

library(nngeo)

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

## Reading layer `nc' from data source `/home/michael/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
## proj4string:    +proj=longlat +datum=NAD27 +no_defs

nc = st_transform(nc, 2163)
pol = st_centroid(nc)

## Warning in st_centroid.sf(nc): st_centroid assumes attributes are constant
## over geometries of x

pol = st_buffer(pol, 10000)
pol$pol_id = 1:nrow(pol)
pol = pol[, "pol_id"]
set.seed(2)
pnt = st_sample(nc, 150)
pnt = st_as_sf(pnt)
pnt$point_id = 1:nrow(pnt)

# Plot sample data
plot(st_geometry(pol))
plot(pnt, add = TRUE, col = "red")

# Find nearest point per polygon
matched = st_join(pol, pnt, join = nngeo::st_nn, maxdist = 5000, k = 1, progress = FALSE)
matched

## Simple feature collection with 100 features and 2 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 1434397 ymin: -972247.3 xmax: 2163573 ymax: -638980.7
## epsg (SRID):    2163
## proj4string:    +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
## First 10 features:
##    pol_id point_id                       geometry
## 1       1       36 POLYGON ((1653386 -769725.7...
## 2       2       NA POLYGON ((1684474 -755795.8...
## 3       3       43 POLYGON ((1724447 -755458.3...
## 4       4       76 POLYGON ((2125485 -648980.7...
## 5       5      139 POLYGON ((2006808 -681505.5...
## 6       6       NA POLYGON ((2044114 -678049, ...
## 7       7       NA POLYGON ((2108020 -654859.7...
## 8       8       90 POLYGON ((2066626 -661943.5...
## 9       9      141 POLYGON ((1947392 -700792, ...
## 10     10      134 POLYGON ((1763898 -747402.5...

plot(matched)

# Find nearest polygon per point
matched = st_join(pnt, pol, join = nngeo::st_nn, maxdist = 5000, k = 1, progress = FALSE)
matched

## Simple feature collection with 150 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 1439778 ymin: -951909.3 xmax: 2156285 ymax: -653444.2
## epsg (SRID):    2163
## proj4string:    +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
## First 10 features:
##    point_id pol_id                         x
## 1         1     NA POINT (1564713 -844541.5)
## 2         2     NA POINT (1963419 -909022.2)
## 3         3     29 POINT (1863993 -756996.4)
## 4         4     NA POINT (1551746 -893225.5)
## 5         5     66 POINT (1521780 -917867.4)
## 6         6      8 POINT (2064407 -674241.4)
## 7         7     67 POINT (1846009 -845801.1)
## 8         8     NA POINT (2008213 -946539.7)
## 9         9     21 POINT (2079893 -694914.8)
## 10       10     54 POINT (1932235 -798927.7)

plot(matched)