找到给定点的最近多边形答案

【问题标题】：find the nearest polygon for a given point找到给定点的最近多边形
【发布时间】：2020-12-28 11:17:06
【问题描述】：

我有一个SpatialPointsDataFrame 和一个SpatialPolygons。我想检查SpatialPointsDataFrame 中的每个点，它位于SpatialPolygons 中的哪个多边形。

我可以使用sp::over 来实现这一点：

但是对于SpatialPointsDataFrame 中的某些点位于边缘的情况或在多边形之外，在这种情况下，我想从 SpatialPolygons。这是示例数据集：

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

所有带有 NA 的行都是我需要分配最近的多边形的行。

【问题讨论】：

标签： r sf sp

【解决方案1】：

对于 lon/lat 数据（如您的示例），您可以使用 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

【讨论】：

【解决方案2】：

sf 包的st_join() 函数可以将点分配给最近的多边形。

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)

【讨论】：

【解决方案3】：

如果您可以转换为sf 对象，您可以使用以下内容找到距多边形外每个点最近的多边形：

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！

【讨论】：