This document introduces the nngeo package. The nngeo package includes functions for spatial join of layers based on k-nearest neighbor relation between features. The functions work with spatial layer object defined in package sf, namely classes sfc and sf.
CRAN version:
install.packages("remotes")
remotes::install_github("michaeldorman/nngeo")
GitHub version:
install.packages("nngeo")
The nngeo package comes with three sample datasets:
citiestownswaterThe cities layer is a point layer representing the location of the three largest cities in Israel.
cities
#> Simple feature collection with 3 features and 1 field
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 34.78177 ymin: 31.76832 xmax: 35.21371 ymax: 32.79405
#> geographic CRS: WGS 84
#> name geometry
#> 1 Jerusalem POINT (35.21371 31.76832)
#> 2 Tel-Aviv POINT (34.78177 32.0853)
#> 3 Haifa POINT (34.98957 32.79405)
The towns layer is another point layer, with the location of all large towns in Israel, compiled from a different data source:
towns
#> Simple feature collection with 193 features and 4 fields
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 34.27 ymin: 29.56 xmax: 35.6 ymax: 33.21
#> geographic CRS: WGS 84
#> First 10 features:
#> name country.etc pop capital geometry
#> 12 'Afula Israel 39151 0 POINT (35.29 32.62)
#> 17 'Akko Israel 45606 0 POINT (35.08 32.94)
#> 40 'Ar'ara Israel 15841 0 POINT (35.1 32.49)
#> 41 'Arad Israel 22757 0 POINT (35.22 31.26)
#> 43 'Arrabe Israel 20316 0 POINT (35.33 32.85)
#> 52 'Atlit Israel 4686 0 POINT (34.93 32.68)
#> 103 'Eilabun Israel 4296 0 POINT (35.4 32.83)
#> 104 'Ein Mahel Israel 11014 0 POINT (35.35 32.72)
#> 105 'Ein Qiniyye Israel 2101 0 POINT (35.15 31.93)
#> 112 'Ilut Israel 6536 0 POINT (35.25 32.72)
The water layer is an example of a polygonal layer. This layer contains four polygons of water bodies in Israel.
water
#> Simple feature collection with 4 features and 1 field
#> geometry type: POLYGON
#> dimension: XY
#> bbox: xmin: 34.1388 ymin: 29.45338 xmax: 35.64979 ymax: 33.1164
#> geographic CRS: WGS 84
#> name geometry
#> 1 Red Sea POLYGON ((34.96428 29.54775...
#> 2 Mediterranean Sea POLYGON ((35.10533 33.07661...
#> 3 Dead Sea POLYGON ((35.54743 31.37881...
#> 4 Sea of Galilee POLYGON ((35.6014 32.89248,...
Figure \ref{fig:layers} shows the spatial configuration of the cities, towns and water layers.
plot(st_geometry(water), col = "lightblue")
plot(st_geometry(towns), col = "grey", pch = 1, add = TRUE)
plot(st_geometry(cities), col = "red", pch = 1, add = TRUE)
st_nn functionThe main function in the nngeo package is st_nn.
The st_nn function accepts two layers, x and y, and returns a list with the same number of elements as x features. Each list element i is an integer vector with all indices j for which x[i] and y[j] are nearest neighbors.
For example, the following expression finds which feature in towns[1:5, ] is the nearest neighbor to each feature in cities:
nn = st_nn(cities, towns[1:5, ], progress = FALSE)
#> lon-lat points
nn
#> [[1]]
#> [1] 4
#>
#> [[2]]
#> [1] 3
#>
#> [[3]]
#> [1] 2
This output tells us that towns[4, ] is the nearest among the five features of towns[1:5, ] to cities[1, ], etc.
st_connect functionThe resulting nearest neighbor matches can be visualized using the st_connect function. This function builds a line layer connecting features from two layers x and y based on the relations defined in a list such the one returned by st_nn:
l = st_connect(cities, towns[1:5, ], ids = nn)
#> Calculating nearest IDs
#> Calculating lines
l
#> Geometry set for 3 features
#> geometry type: LINESTRING
#> dimension: XY
#> bbox: xmin: 34.78177 ymin: 31.26 xmax: 35.22 ymax: 32.94
#> geographic CRS: WGS 84
#> LINESTRING (35.21371 31.76832, 35.22 31.26)
#> LINESTRING (34.78177 32.0853, 35.1 32.49)
#> LINESTRING (34.98957 32.79405, 35.08 32.94)
Plotting the line layer l gives a visual demonstration of the nearest neighbors match, as shown in Figure \ref{fig:st_connect}.
plot(st_geometry(towns[1:5, ]), col = "darkgrey")
plot(st_geometry(l), add = TRUE)
plot(st_geometry(cities), col = "red", add = TRUE)
text(st_coordinates(cities)[, 1], st_coordinates(cities)[, 2], 1:3, col = "red", pos = 4)
text(st_coordinates(towns[1:5, ])[, 1], st_coordinates(towns[1:5, ])[, 2], 1:5, pos = 4)
![Nearest neighbor match between \texttt{cities} (in red) and \texttt{towns[1:5, ]} (in grey)](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAWgAAAFoCAIAAAD1h/aCAAAACXBIWXMAAAsSAAALEgHS3X78AAARMElEQVR4nO3df1DUdf7A8RfLz+QQwaA0OZCmU+sq7bQMEe80sxsPx+Q8bK5/tJ93UzMl2NBQf2jN5HHGqKVOjDN+u8au5urm6zalTcEyJ6VjX1pFKU0hsSSRatnFhV1Zdr9/4CkQVi/57H5gP8/H+Aewn4HXH8zT92s/LMSEQiEBAA2b2QMAGH0IBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcuBK7d++ePn16ampqfn7+kSNHzB4HkUY4oPb1118vX768oqLi9OnTeXl5xcXFZk+ESIsJhUJmz4BR5o033tixY8f7778vIm63Oy0t7fvvvx83bpzZcyFyCAeG1tnZ2dLS0tvbO2HChMzMzP4Peb1ev9+fnp4uIrW1tQ888EBTU5NJY8IcrCoYQnt7+6effpqZmZmVldXS0nLs2LH+jyYnJ/dVw26333fffZWVlSaNCdPEmT0ARqLPPvssLy8vPj5eRNLT0x0Oxw033GCzXfpvxu12P/roowcPHnzrrbfmzJlj3qQwBycODCEUCsXHx7/zzjuvv/66iKSkpHR1dV18tKenZ9GiRenp6Q0NDVTDmggHhhAbG+v1etevX3/PPfeEQiGPxzNmzJiLj+7atSsQCFRWVgaDQb/f7/f7TRwVpmBVwRBuuumm5557Lj8///z58/v27cvJyem/p9TX19fX1yclJV38SEdHR2pqqhmTwhzcVcEQenp65syZU1VVlZSUdO2113KrFYMQDgzh5ZdfDoVCjz/+uNmDYIQiHBisu7t7wYIFDocjMTHR7FkwQvHkKAbbuHHjww8/HOFqzJs3L+a/CgsLI/mlcQV4chQDuN1uu92+d+/eCH/d5ubmpqama665RkRiY2Mj/NWhRTgwQEVFRUlJSVyc8d8YPp+vo6Pjqquu+uH9l76HcnNzDf+iCBNWFVzS3t7+8ccfFxUVGf6Zm5qaDhw44HK5jh8//tFHHwWDwf6Pnjx5Mi4ububMmWlpaQsXLmxubjZ8ABiLcOCS559/vqysLCYmxthP6/V6W1tbCwoKpk2bNnPmzOuuu+748eP9L/B4PHl5eW+++ebZs2dvvPHGFStWGDsADMddFVzQ0tLy4IMPfvDBB4Z/5q+++ur8+fPXX39937s9PT0HDhy43M+qnzp1Kjs7+7vvvut7HR1GJk4cuGDt2rXPPPNMOD5zYmJid3f3xXd9Pt+gWzb79++vq6vrezs+Pj42NpY7wSMc4YCIyBdffOFyuebNmzfgo9u2SXa2pKVJYaG0tV3xJ7/66qu//fbbtrY2Eenq6nI6nYOeBz1z5kxRUdHRo0eDweDGjRsXLlyYnJx8xV8OEcCqAhGR4uLisrKyGTNmXPpQY6MsWCAOh0yYIA88IGlpsn37FX9+v9//+eefd3R0JCQkTJkyZfz48YMu2Lx580svveTxeObOnbtly5a++7IYsQiHtZ06JV98cdjvf2Hnzr5X0F9it8t//iMbNoiI/OtfsnWrOBymzIgRiJ/jsLCnn5Zjx+S22zq3bHn5978f/OiSJbJkibhccuiQbNsm3OlAPzzHYVXV1XLunPz733W//e3/FBamp6XJe+8NcZnTKU8+KceOyaxZER8RIxfhsKq6OikqEpGKiopnn31W/vhH+e99jQHmzxenUyorOXGgP8JhVRkZfTdKdu7cmZWVJWfOSEbGgAtefFGqqi68feed0twsgUDEp8QIRTisqqhIXnpJGhtTUlLk88+lslL+9KcBF0yeLJs2yZEj4vHIhg2SlydheAELRim+Fazqmmtkxw55/nk5dUomTZLt2+W66wZcsGyZnDghS5eKyyWzZ8trr5k0KEYibscCUGNVAaBGOACoEQ4AaoQDgBrhAKBGOACoEQ4AaoQDgBrhAKBGOACoEQ4AaoQDgBrhgMk2b97c2tpq9hTQIRwwTSAQqK6uLi8v7+zsNHsW6PD7OBAuoVCotbXV6/WmpqYO+ecOVq1a5XQ6fT5f5GfDMPH7OBAWwWCwrq4uIyMjNTW1vb3d5/PdcccdQ145YcKE2traKVOmRHhCDAerCsKiubl50qRJ06ZNmzhx4q233pqUlNQ2jL8Fh5GGcCAs3G53Rr/ffpyRkdHR0WHiPDAW4UBYJCcnezyei+96PB7+HGw0IRwIi9zc3KNHj/Y9OXry5Mm2traJEyeaPRQMQzgQFgkJCXPnzu3o6GhsbOzp6cnPz7fZ+GaLHtxVAaDGfwIA1AgHADXCAUCNcABQIxwA1AgHADXCgSu3bdu27OzstLS0wsJCXopiKYQDV6ixsXHt2rV79uz58ssvExISysvLzZ4IkUM48BN6enq8Xu8Pf1Cwqanp/vvvnzZt2rhx41asWNHU1GTKeDAFPzmKywoGg/X19T6fLzEx8dy5czfffHP/F7z2cblchw4dWrduXXFx8SOPPGLKnIg8fgMYLuvYsWPjx4/Pzc0VkZ6enr179xYUFMTFDfiecTqdJSUlZ8+enTVrlkljwgSsKris9vb2nJwcj8dTW1sbHx+fmZnpcrkGXTN//nyn01lZWblixQpThoQpCAcuy2azBYPBkpISr9crIoFAoP9x48UXX6yqqup7+84772xubg4EAuYMiogjHLisrKys7du3BwKBxYsXd3R0uFyu1NTUi49Onjx506ZNR44c8Xg8GzZsyMvLG7TFIIrx5Cguy+PxFBQUrFu3buzYsUlJSbfccsuYMWP6X1BRUVFVVeVyuWbPnr1169bs7GyzRkWEEQ5c1kMPPbR48eKlS5eaPQhGHFYVDK26urq3t5dqYEicODAEj8dz99137969Oy0tzexZMBJx4sAQSkpKysrKqAYuh3BgsOrq6kAgwJKCH8GqggFYUvBzcOLAACwp+DkIBy5hScHPxKqCC/qWlPfeey89Pd3sWTDSceLABaWlpWVlZVQDPwfhgIhIdXX1+fPnWVLwM7GqgCUFapw4wJICNcJhdTU1NSwp0GJVsTSv17to0SK73c5xAyqcOCxt9erVJSUlVANahMO6ampq/H7/vffea/YgGH1YVSyKJQXDwYnDolhSMByEw4pqamp8Ph9LCq4Yq4rlsKRg+DhxWE5JSQlLCoaJcFhLTU1Nd3c3SwqGiVXFQlhSYBROHBbCkgKjEA6rqKmp6erqYkmBIVhVLIElBcbixGEJpaWlq1evphowCuGIfjU1NV6vd9myZWYPgujBqhLlWFIQDpw4ohxLCsKBcESzmpqac+fOsaTAcKwqUYslBeHDiSNqsaQgfAhHdHI4HCwpCB9WlSjk9Xrvuusuu92ekZFh9iyITpw4otCaNWvWrFlDNRA+hCPaOByOzs5OlhSEFatKVGFJQWRw4ogqLCmIDMIRPRwOh8fjYUlBBLCqRAmWFEQSJ44o8dRTT5WWlg6oxu7dMn26pKZKfr4cOWLeaIhChCMaOBwOt9tdVFR06UNffy3Ll0tFhZw+LXl5Ulxs3nSIQqwqo1YoJLW10tzcff31859+eteuXZmZmZcefeMN2bFD3n9fRMTtlrQ0+f57GTfOrGERZThxjE6BgCxbJnv2SGLi/z322OspKQOqISKFhfLPf1542+mUyZOpBgxEOEanV1+V3/1O/va32kmTXpk+ffKvfy3V1QMuSE6Wvpe32e1y331SWWnKmIhWcWYPgCvidMpf/iIif//731999VX55BOpr5cFCwZc43bLo4/KwYPy1lsyZ445cyJKEY7RKStLTpyQm2569913RUSOH5df/nLABT09smiR/OY30tAg8fGmzIgoRjhGp5UrZelS+cUv5JZbpK5O3n5b9uwZcMGuXRIISGWlBIPi94uIJCaaMimiEs9xjE6ZmfL22/Lhh/LXv0pjo9jtctVVAy6or5f6eklKuvTP7TZpVkQhbscCUOPEAUCNcABQIxwA1AgHADXCAUCNcABQIxwA1AgHADXCAUCNcABQIxwA1AgHADXCAUCNcABQIxwA1AgHADXCAUCNcABQIxwA1AgHADXCAUCNcABQIxwA1AgHADXCAUCNcABQIxwA1AgHADXCAUCNcCBSNm+W1lazh4AxCAfCLxCQ6mopL5fOTrNHgTHizB4AFrBqlTid4vOZPQcMQzgwbC0tsmWLfPONzJ4tDz0kCQmDL/jHP0REJkyI/GgIE1YVDM+JE/LnP8sf/iAvvCBxcbJ8uQSDZs+EsCMcGJ5Nm2TTJikokEmT5JFHJDdXDhwweyaEHeHA8LS0yK9+5fP5Vq5cKSIydaqcPGnySAg/woHhufVWqa0NBoNnz54VEamtlenTzZ4JYceToxie0lJZsiTh8OE5p0/LqlWSmytTp5o9E8KOcGB4UlPlww9j7PYUv18ee0xuu+2yV37zTQTHQnixqmDY4uNj7r33fydO/LFqILoQDhjAZrMFuQtrJYQDgBrhAKBGOACoEQ4AaoQDgBrhAKBGOACoEQ4AaoQDgBrhAKBGOACoEQ4AaoQDgBrhAKBGOACoEQ4AaoQDgBrhAKBGOACoEQ4AaoQDgBrhAKBGOACoEQ4YJhQKmT0CIoRwwBixsbH8TSbrIBwwRmxsbG9vr9lTIEIIB4zBX4G0FMIBY3DisBTCAWMQDkshHDAGq4qlEA4YgxOHpRAOGINwWArhgDFYVSyFcMAYnDgshXDAGITDUggHjMGqYimEA8bgxGEphAPGIByWQjhgDJvNRjisg3DAGLys3lIIB4zBqmIphAPGIByWQjhgDG7HWgrhgDE4cVgK4YAxCIelEA4Yg1XFUggHjMGJw1IIB4zBD4BZCuGAMfgBMEshHDAGq4qlEA4Yg1XFUggHjMGqYimEA8b4yVWlq6vriSeeiNg8CCvCAWP85Kry5JNPvvLKKxGbB2FFOGCMH19V7HZ7Y2NjJOdBWBEODFdvb29DQ0NTU1NDQ8Phw4d/eO44c+ZMeXn5tm3bTBkP4UA4MFwHDx4cO3bs1KlTp06dmpycfOjQof6PhkKhlStXrl+/PjMz06wJYTjCgWEJhUKdnZ05OTl9q0pubq7H4+l/wdatW7OzsxcvXmzWhAgHwoFhCYVCNptNRObOnXv77beLSExMTP8nO/bv379z586UlJTc3Fyfz5eSkrJv3z7TxoVB4sweAKObzWaz2Wwej2fKlCki4na74+Li+lLS57XXXut7o62tLScnp7Oz05xBYSjCgeGaMWPGJ598MmbMGBHp7u6eNWuW2RMh7GJCoZDZMyAaeL3emJiYvnwg6hEOAGo8OQpAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFAjXAAUCMcANQIBwA1wgFA7f8BU/NPX/K6O1kAAAAASUVORK5CYII= "Nearest neighbor match between \texttt{cities} (in red) and \texttt{towns[1:5, ]} (in grey)")
The st_nn can also return the complete logical matrix indicating whether each feature in x is a neighbor of y. To get the dense matrix, instead of a list, use sparse=FALSE.
nn = st_nn(cities, towns[1:5, ], sparse = FALSE, progress = FALSE)
#> lon-lat points
nn
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] FALSE FALSE FALSE TRUE FALSE
#> [2,] FALSE FALSE TRUE FALSE FALSE
#> [3,] FALSE TRUE FALSE FALSE FALSE
k>0It is also possible to return any k-nearest neighbors, rather than just one. For example, setting k=2 returns both the 1st^ and 2nd^ nearest neighbors:
nn = st_nn(cities, towns[1:5, ], k = 2, progress = FALSE)
#> lon-lat points
nn
#> [[1]]
#> [1] 4 3
#>
#> [[2]]
#> [1] 3 1
#>
#> [[3]]
#> [1] 2 5
nn = st_nn(cities, towns[1:5, ], sparse = FALSE, k = 2, progress = FALSE)
#> lon-lat points
nn
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] FALSE FALSE TRUE TRUE FALSE
#> [2,] TRUE FALSE TRUE FALSE FALSE
#> [3,] FALSE TRUE FALSE FALSE TRUE
Using returnDist=TRUE the distances list is also returned, in addition the the neighbor matches, with both components now comprising a list:
nn = st_nn(cities, towns[1:5, ], k = 2, returnDist = TRUE, progress = FALSE)
#> lon-lat points
nn
#> $nn
#> $nn[[1]]
#> [1] 4 3
#>
#> $nn[[2]]
#> [1] 3 1
#>
#> $nn[[3]]
#> [1] 2 5
#>
#>
#> $dist
#> $dist[[1]]
#> [1] 56364.74 80742.62
#>
#> $dist[[2]]
#> [1] 53968.63 76186.87
#>
#> $dist[[3]]
#> [1] 18265.72 32476.24
Finally, the search for nearest neighbors can be limited to a search radius using maxdist. In the following example, the search radius is set to 50,000 meters (50 kilometers). Note that no neighbors are found within the search radius for cities[2, ]:
nn = st_nn(cities, towns[1:5, ], k = 2, maxdist = 50000, progress = FALSE)
#> lon-lat points
nn
#> [[1]]
#> integer(0)
#>
#> [[2]]
#> integer(0)
#>
#> [[3]]
#> [1] 2 5
The st_nn function can also be used as a geometry predicate function when performing spatial join with sf::st_join. For example, the following expression spatially joins the two nearest towns[1:5, ] features to each cities features, using a search radius of 50 km:
cities1 = st_join(cities, towns[1:5, ], join = st_nn, k = 2, maxdist = 50000)
#> lon-lat points
Here is the resulting layer:
cities1
#> Simple feature collection with 4 features and 5 fields
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 34.78177 ymin: 31.76832 xmax: 35.21371 ymax: 32.79405
#> geographic CRS: WGS 84
#> name.x name.y country.etc pop capital geometry
#> 1 Jerusalem <NA> <NA> NA NA POINT (35.21371 31.76832)
#> 2 Tel-Aviv <NA> <NA> NA NA POINT (34.78177 32.0853)
#> 3 Haifa 'Akko Israel 45606 0 POINT (34.98957 32.79405)
#> 3.1 Haifa 'Arrabe Israel 20316 0 POINT (34.98957 32.79405)
Here is another example, finding the 10-nearest neighbor towns features for each cities feature:
x = st_nn(cities, towns, k = 10)
#> lon-lat points
l = st_connect(cities, towns, ids = x)
The result is visualized in Figure \ref{fig:cities_towns}.
plot(st_geometry(l))
plot(st_geometry(cities), col = "red", add = TRUE)
plot(st_geometry(towns), col = "darkgrey", add = TRUE)
Nearest neighbor search also works for non-point layers. The following code section finds the 20-nearest towns features for each water body in water[-1, ].
nn = st_nn(water[-1, ], towns, k = 20, progress = FALSE)
#> lines or polygons
Again, we can calculate the respective lines for the above result using st_connect. Since one of the inputs is line/polygon, we need to specify a sampling distance dist, which sets the resolution of connecting points on the shape exterior boundary.
l = st_connect(water[-1, ], towns, ids = nn, dist = 100)
#> Calculating nearest IDs
#> Calculating lines
The result is visualized in Figure \ref{fig:water_towns}.
plot(st_geometry(water[-1, ]), col = "lightblue", border = "grey")
plot(st_geometry(towns), col = "darkgrey", add = TRUE)
plot(st_geometry(l), col = "red", add = TRUE)
