2

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

enter image description here

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.

  • 1
    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... – Spacedman Aug 23 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? – Kenji Aug 23 at 11:02
  • 1
    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? – Spacedman Aug 23 at 20:57
  • MatchingR. I'll have a look at that! Thank you. – Kenji Aug 24 at 4:36
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")

Sample data

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

Nearest point per polygon

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

Nearest polygon per point

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.