Search code examples
rr-sfr-sp

find the nearest polygon for a given point

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)

enter image description here

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.

Solution

  • 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 !