使用R中的terra，返回最接近点的多边形。

Question

我正在进行点在多边形分析

library(terra)
library(rnaturalearth)
  
crdref <- "+proj=longlat +datum=WGS84"

lonlat<- structure(c(-123.115684, -81.391114, -74.026122, -122.629252, 
                      -159.34901, 7.76101, 48.080979, 31.159987, 40.621058, 47.50331, 
                       21.978049, 36.90086), .Dim = c(6L, 2L), 
                      .Dimnames = list(NULL,c("longitude", "latitude")))

pts <- vect(lonlat, crs = crdref)
world_shp <- rnaturalearth::ne_countries()
world_shp <- terra::vect(world_shp, crs = crdref)
world_shp <- terra::project(world_shp, crdref)

plot(world_shp)
points(pts, col = "red", pch = 20)      

所有这些点都位于多边形的边缘上，因此当我尝试提取每个点所在的多边形时，会得到一个NA。

e <- terra::extract(world_shp, pts) 
e$sovereignt
NA
  

有没有一种方法可以使用terra包返回每个点的最近多边形？

- 89_Simple

---

我在 > pts <- vect(lonlat, crs = crdref) 上遇到了错误 Error in h(simpleError(msg, call)) : error in evaluating the argument 'x' in selecting a method for function 'vect': object 'lonlat' not found` - Mossa

抱歉，我已经修改了问题中的错误。 - 89_Simple

Answer 1

你可以使用sf包中的st_join函数：

library(rnaturalearth)
library(tidyverse)
library(sf)
library(rgdal)

devtools::install_github("ropensci/rnaturalearthhires")
library(rnaturalearthhires)


# create points dataframe
points <- tibble(longitude = c(123.115684, -81.391114, -74.026122, -122.629252,
                               -159.34901, 7.76101),
                 latitude = c(48.080979, 31.159987, 40.621058, 47.50331,
                              21.978049, 36.90086)) %>% 
  mutate(id = row_number())

# convert to spatial
points <- points %>% 
  st_as_sf(coords = c("longitude", "latitude"),
           crs = 4326)

# load Natural Earth data
world_shp <- rnaturalearth::ne_countries(returnclass = "sf",
                                         scale = "large")


# Plot points on top of map (just for fun)
ggplot() +
  geom_sf(data = world_shp) +
  geom_sf(data = points, color = "red", size = 2)


# check that the two objects have the same CRS
st_crs(points) == st_crs(world_shp)


# locate for each point the nearest polygon
nearest_polygon <- st_join(points, world_shp,
                           join = st_nearest_feature)

nearest_polygon

- NorthNW

---

请尝试在 st_make_valid 后绘制 world_shp 的结果。 - Mossa

亲爱的主啊。那些多边形绝对是非常无效的。 - NorthNW

我已根据@Mossa的评论编辑了我的答案。这是ne_countries()函数的错误，可以通过指定returnclass和scale来修复。这又需要安装rnaturalearthhires包。 - NorthNW

好的，谢谢。我之前做过类似的解决方案，但想检查一下 terra 包中是否有可用的内容。不过还是谢谢。稍后我会接受这个答案。 - 89_Simple

terra包主要用于栅格数据。您试图将多边形视为栅格数据（使用extract）。我认为您完全可以采用sf方法，没有任何理由不这样做。 - NorthNW

我喜欢这个，投了一票；非常感谢！ - Mossa

Answer 2

您可以使用terra::nearest

示例数据

library(terra)
v <- vect(system.file("ex/lux.shp", package="terra"))
pts <- matrix(c(5.812,6.577,5.864,50.126,49.774,49.488), ncol=2)
pts <- vect(pts, crs = crs(w))

解决方案

pv = as.points(v)
n = nearest(pts, pv)

i = pv$NAME_2[n$to_id]
i
#[1] "Clervaux"         "Echternach"       "Esch-sur-Alzette"

插图

plot(v, border="gray", lwd=3)
text(v, "NAME_2", cex=.6)
points(pts)
lines(pts, n, col="blue", lwd=3)
text(pts, i, cex=.6, col="red", halo=TRUE, adj=c(.5,1.5))

对于上面的例子，你也可以这样做

n = nearest(pts, v)

但在大多数情况下，您需要的是到边界的距离，而不是质心

n = nearest(pts, v, centroids=FALSE)

当前版本返回错误的ID（对于节点而不是多边形），因此无法使用。

- Robert Hijmans