我正在做多邊形分析
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包為每個點回傳最近的多邊形
uj5u.com熱心網友回復：
您可以st_join從sf包中使用：
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
轉載請註明出處，本文鏈接：https://www.uj5u.com/qukuanlian/488574.html
上一篇：R按組用NA替換最后n個值
下一篇：將表的輸出匯總到簡單的列