我有一些海上的地理点和一张瑞典海域的地图。我想知道每个点位于哪个海域,对于海域之外的点,到最近海域的距离,包括最近海域的ID。
library(terra)
#> terra 1.7.71
map <- vect("/vsizip//vsicurl/https://www.smhi.se/polopoly_fs/1.140307!/Havsomr_SVAR_2016_3b.zip")
pts <- data.frame(id = 11:25,
lon = c(15.29416656, 15.10333347, 15.27083302, 15.27250004, 15.1291666,
15.31333351, 15.24166679, 15.30716705, 15.32333374, 15.38833332,
15.41416645, 15.40666676, 15.27166653, 15.26083374, 15.1916666),
lat = c(56.1566658, 56.15000153, 56.15250015, 56.13999939, 56.16383362,
56.15333176, 56.1558342, 56.17416763, 56.12366486, 56.14250183,
56.15916824, 56.17666626, 56.17666626, 56.17300034, 56.1570015)) |>
vect(geom = c("lon", "lat"), crs = "EPSG:4326") |>
project("EPSG:3006")
正如您所看到的,当您运行这个简单的代码时,
to_id
始终是NA
并且属于逻辑类型。我想要将 "HID"
中的 map
列放入 to_id
列(或地图中的其他 ID)。我可以通过某种方式得到它吗?
npol <- nearest(pts, map[, "HID"], centroids = FALSE)
npol
#> class : SpatVector
#> geometry : points
#> dimensions : 15, 5 (geometries, attributes)
#> extent : 506419.7, 525716.6, 6219890, 6225817 (xmin, xmax, ymin, ymax)
#> coord. ref. : SWEREF99 TM (EPSG:3006)
#> names : from_id from_x from_y to_id distance
#> type : <int> <num> <num> <logical> <num>
#> values : 1 5.183e+05 6.224e+06 <NA> 0
#> 2 5.064e+05 6.223e+06 <NA> 0
#> 3 5.168e+05 6.223e+06 <NA> 0
不幸的是,
?terra::nearest
的文档没有提供nearest()
的任何示例,尤其是centroids = FALSE
。就我个人而言,如果该点位于多边形内,我预计距离为零,否则对应于到多边形最近节点的距离 - 但我无法确定这意味着什么。显然这似乎不符合我的期望。
但是,对于这些数据,我们可以简单地使用
terra::distance()
计算完整的距离矩阵,然后选择我们感兴趣的信息:相关多边形的索引(您可以使用它来索引来自 "HID"
的 map
列) )和距离本身:
distm <- distance(pts, map[, "HID"])
pts[, "to_id"] <- apply(distm, MARGIN = 1, FUN = which.min)
pts[, "to_HID"] <- map[pts$to_id, ]$HID
pts[, "distance"] <- apply(distm, MARGIN = 1, FUN = min)
pts
#> class : SpatVector
#> geometry : points
#> dimensions : 15, 4 (geometries, attributes)
#> extent : 506419.7, 525724.3, 6219890, 6225817 (xmin, xmax, ymin, ymax)
#> coord. ref. : SWEREF99 TM (EPSG:3006)
#> names : id to_id to_HID distance
#> type : <int> <int> <chr> <num>
#> values : 11 279 560940-151740 0
#> 12 288 560850-150580 0
#> 13 279 560940-151740 0
根据此输出,您的 15 个点中有 3 个在岸上。
编辑:
只是出于好奇而使用
{sf}
来重现这一点:
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
map_sf <- st_as_sf(map)
pts_sf <- st_as_sf(pts)
st_nearest_feature(pts_sf, map_sf)
#> [1] 279 288 279 371 289 279 327 279 371 244 280 280 279 279 327
st_distance(pts_sf, map_sf) |> apply(MARGIN = 1, FUN = min)
#> [1] 0.000000 0.000000 0.000000 0.000000 58.765524 0.000000 0.000000
#> [8] 27.636855 0.000000 0.000000 7.755143 0.000000 0.000000 0.000000
#> [15] 0.000000
最后一次使用
{nngeo}
:
library(nngeo)
pts_nn <- st_nn(pts_sf, map_sf, returnDist = TRUE)
lapply(pts_nn, unlist)
#> $nn
#> [1] 279 288 279 371 289 279 327 279 371 244 280 280 279 279 327
#>
#> $dist
#> [1] 0.000000 0.000000 0.000000 0.000000 58.765524 0.000000 0.000000
#> [8] 27.636855 0.000000 0.000000 7.755143 0.000000 0.000000 0.000000
#> [15] 0.000000