#r #sf #sp
#r #sf #sp
Вопрос:
У меня есть 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 — это те, которые мне нужны, чтобы назначить ближайший полигон.
Ответ №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!
Ответ №2:
st_join() Функция пакета sf может назначить точки ближайшему полигону.
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:
С данными 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