如何使用 sf 包计算到 R 中最近多边形的距离

Question

我想计算一个点到最近的多边形的距离。

library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(rnaturalearth)
library(ggplot2)

pts <- st_as_sf(
    data.frame(longitude = 3.7, latitude = 52.3667),
    coords = c("longitude", "latitude"),
    crs = 4326
  )

ne_land <- ne_download(
  category = "physical",
  type = "land",
  returnclass = "sf",
  scale = "small"
)
#> OGR data source with driver: ESRI Shapefile 
#> Source: "/tmp/Rtmplq0wPd", layer: "ne_110m_land"
#> with 127 features
#> It has 3 fields

此图显示了点相对于陆地多边形的位置。

ggplot() +
  geom_sf(data = ne_land) +
  geom_sf(data = pts) +
  coord_sf(xlim = c(2, 5), ylim = c(50, 54))

然后我试图通过以下方式获取到最近多边形的距离：

nearest <- st_nearest_feature(pts, ne_land)
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar
nearest
#> [1] 113

dist <-  st_distance(pts, ne_land[nearest, ], by_element=TRUE)
dist
#> 0 [m]

如图所示，dist 等于 0 米。根据之前的图表，这不是我所期望的结果。任何帮助表示赞赏。 编辑： 根据最后一个答案，我尝试了其他一些位置和更高的陆地分辨率，也得到了 0 米的距离。

ne_land <-
  rnaturalearth::ne_download(
    category = "physical",
    type = "land",
    returnclass = "sf",
    scale = "large"
  )
#> OGR data source with driver: ESRI Shapefile 
#> Source: "/tmp/Rtmplq0wPd", layer: "ne_10m_land"
#> with 10 features
#> It has 3 fields

pts <- st_as_sf(
  data.frame(longitude = 4.1833, latitude = 52.1833),
  coords = c("longitude", "latitude"),
  crs = 4326
)

st_distance(pts, st_union(ne_land))
#> although coordinates are longitude/latitude, st_union assumes that they are planar
#> Units: [m]
#>      [,1]
#> [1,]    0

ggplot() + 
  geom_sf(data = st_nearest_points(pts, st_union(ne_land))) + 
  geom_sf(data = ne_land) +
  geom_sf(data = pts) +
  coord_sf(xlim = c(2, 5), ylim = c(50, 54))
#> although coordinates are longitude/latitude, st_union assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar

^{由reprex package (v2.0.0) 于 2021-04-06 创建}

Answer 1

您需要在ne_land 数据上使用st_union() 来找到最近的单个点。

st_distance(pts, st_union(ne_land))

although coordinates are longitude/latitude, st_union assumes that they are planar
Units: [m]
         [,1]
[1,] 36319.29

ggplot() + 
 geom_sf(data = st_nearest_points(pts, st_union(ne_land))) + 
 geom_sf(data = ne_land) +
 geom_sf(data = pts) +
 coord_sf(xlim = c(2, 5), ylim = c(50, 54))

编辑： 以下距离与上述距离有很大不同，因为 pts 已更改为已编辑问题。

我不知道为什么它不适用于大型自然地球数据，但这里有一个解决方法：

st_nearest_points(pts, ne_land) %>% 
  st_length()

#although coordinates are longitude/latitude, st_nearest_points assumes #that they are planar
#11481.03 [m]

edit2： 这也有效，但返回不同的距离。这可能与使用的投影 (crs) 或将多边形简化为线有关。

st_combine(ne_land) %>% st_cast('MULTILINESTRING') %>% st_distance(pts)
Units: [m]
         [,1]
[1,] 10136.87