找到给定点的最近多边形

Question

我有一个 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 的行都是我需要分配最近多边形的行。

Answer 1

如果您可以转换为 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!

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

Answer 3

使用 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

找到给定点的最近多边形

find the nearest polygon for a given point

r

sp

sf