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

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.

• You could loop over the data, matching the nearest point/poly pair then removing them from the data for the next iteration. But you might end up with a final point/poly a long distance away. If that's okay then fine. But to find the set of pairings that minimises the total inter-pair distance sounds a bit like a travelling salesperson problem in complexity terms... Aug 23, 2019 at 10:57
• I know that the points represent one of the buildings in the same 100x100m grid cell. Your solution of iterating and removing would not lead to long distances if I restricted the matches to the same cell. Do you think that this would be feasible on datasets with over 20000 points and polygons? Aug 23, 2019 at 11:02
• I've just tried playing with the Gale-Shapley functions from the `matchingR` package using a negative distance matrix as the utility and am getting interesting and feasible results - any use? Aug 23, 2019 at 20:57
• MatchingR. I'll have a look at that! Thank you. Aug 24, 2019 at 4:36

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
``````
``````## 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
``````
``````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)
``````