find the nearest polygon for a given point

Question

I have a SpatialPointsDataFrame and a SpatialPolygons. I want to check for each point in SpatialPointsDataFrame, which polygon in the SpatialPolygons does it lie.

I can do this using sp::over to achieve this:

However for cases where some of the points in SpatialPointsDataFrame lie either on the edges or outside the polygon and in such cases I want to assign the nearest polygon from the SpatialPolygons. Here's the sample dataset:

set.seed(1)
library(raster)
library(rgdal)
library(rgeos)
  
p <- shapefile(system.file("external/lux.shp", package="raster"))
p2 <- as(1.5*extent(p), "SpatialPolygons")
proj4string(p2) <- proj4string(p)
pts <- spsample(p2, n=10, type="random")
  
## Plot to visualize
plot(p, col=colorRampPalette(blues9)(12))
plot(pts, pch=16, cex=.5,col="red", add = TRUE)

over(pts, p)
  
ID_1       NAME_1 ID_2     NAME_2 AREA
1     1     Diekirch    3    Redange  259
2    NA         <NA>   NA       <NA>   NA
3    NA         <NA>   NA       <NA>   NA
4    NA         <NA>   NA       <NA>   NA
5    NA         <NA>   NA       <NA>   NA
6    NA         <NA>   NA       <NA>   NA
7     3   Luxembourg   10 Luxembourg  237
8     3   Luxembourg    8   Capellen  185
9     2 Grevenmacher    6 Echternach  188
10   NA         <NA>   NA       <NA>   NA
  
  

All rows with NA are the ones I need to assign the nearest polygon.

Labanna

Answer 1

If you're ok to convert to sf objects, you can find the nearest polygon to each point outside the polys using something on this lines:

set.seed(1)
library(raster)
library(rgdal)
library(rgeos)
library(sf)
library(mapview)
  
p <- shapefile(system.file("external/lux.shp", package="raster"))
p2 <- as(1.5*extent(p), "SpatialPolygons")
proj4string(p2) <- proj4string(p)
pts <- spsample(p2, n=10, type="random")
  
## Plot to visualize
plot(p, col=colorRampPalette(blues9)(12))
plot(pts, pch=16, cex=.5,col="red", add = TRUE)

# transform to sf objects
psf   <- sf::st_as_sf(pts) %>% 
    dplyr::mutate(ID_point = 1:dim(.)[1])
polsf <- sf::st_as_sf(p)

# remove points inside polygons
in_points  <- lengths(sf::st_within(psf,polsf))
out_points <- psf[in_points == 0, ]

# find nearest poly
nearest <- polsf[sf::st_nearest_feature(out_points, polsf) ,]  %>% 
    dplyr::mutate(id_point = out_points$ID)
nearest

> Simple feature collection with 6 features and 6 fields
> geometry type:  POLYGON
> dimension:      XY
> bbox:           xmin: 5.810482 ymin: 49.44781 xmax: 6.528252 ymax: 50.18162
> geographic CRS: WGS 84
>   ID_1       NAME_1 ID_2           NAME_2 AREA                       geometry id_point
> 1    2 Grevenmacher    6       Echternach  188 POLYGON ((6.385532 49.83703...        1
> 2    1     Diekirch    1         Clervaux  312 POLYGON ((6.026519 50.17767...        2
> 3    3   Luxembourg    9 Esch-sur-Alzette  251 POLYGON ((6.039474 49.44826...        5
> 4    2 Grevenmacher    7           Remich  129 POLYGON ((6.316665 49.62337...        6
> 5    3   Luxembourg    9 Esch-sur-Alzette  251 POLYGON ((6.039474 49.44826...        7
> 6    2 Grevenmacher    6       Echternach  188 POLYGON ((6.385532 49.83703...        9
> 

#visualize to check
mapview::mapview(polsf["NAME_2"]) + mapview::mapview(out_points)

HTH !

Drain

Answer 2

The sf package's st_join() function can assign the points to the nearest polygon.

set.seed(1)
library(raster)
library(rgdal)
library(rgeos)
library(sf)

p <- shapefile(system.file("external/lux.shp", package="raster"))
p2 <- as(1.5*extent(p), "SpatialPolygons")
proj4string(p2) <- proj4string(p)
pts <- spsample(p2, n=10, type="random")

pts <- st_as_sf(pts)
p <- st_as_sf(p)

pts <- st_join(pts, p, join = st_nearest_feature)

Petta

Answer 3

With lon/lat data (as in your example), you can use geosphere::dist2Line

library(geosphere)
dist2Line(pts, p)
#         distance      lon      lat ID
# [1,]  1161.79335 5.864012 49.50125 10
# [2,]    64.55319 5.985080 49.45904 10
# [3,]  5929.42723 6.190536 49.97124  4
# [4,]  8295.91091 6.516485 49.72418  8
# [5,]  7471.54277 5.863943 50.06754  1
# [6,]  5522.13076 6.528252 49.80857  6
# [7,] 28518.21197 6.524343 49.81309  6
# [8,] 25964.73248 6.120430 50.16320  1
# [9,]  1602.34368 6.269915 49.67434  8
#[10,] 18250.14130 5.859116 50.06171  1

Jessamine