对于地理点的全球数据集（在 lat/long 中）如何找到最近的邻居来解释我们星球的球形性质

Question

所以我们在不考虑投影问题的情况下完成了这项工作。问题是在哪里（以及如何）最好地添加重新投影，以便函数 returns 以千米为单位的值而不是当前度数：

library(raster)
library(purrr)
library(sf)

#example presence data from model
r1 <- raster(nrow=360, ncol=720)
crs(r1) <- "+proj=longlat +datum=WGS84 +no_defs"
values(r1) <- rbinom(ncell(r1), 2, 0.01)

r1_points <- rasterToPoints(r1)
r1_df <- data.frame(r1_points)
r1_presence <- r1_df %>% dplyr::filter(layer==1)

#example survey data
survey_points <- cbind(rnorm(50) * 5 + 10, rnorm(50) + 50)
pt2 <- st_multipoint(cbind(survey_points[,1], survey_points[,2]))

#distance between each modelled presence (pt1) and survey point (pt2)
get_distances <- function(i, pt2, df) {
  pt1 <- st_multipoint(cbind(df[i, 1], df[i, 2]), dim = "XY")
  a <- st_nearest_points(pt1, pt2)
  return(st_length(a))
}

#loop for all modelled presences
output <- map_dbl(1:nrow(r1_presence), get_distances, pt2, r1_presence)

理想情况下，一个完美的答案是扩展 get_distances 函数以添加一个新选项，该选项可以进行适当的重新投影和 returns 以千米为单位的值。

这里可能有几种不同的方法，我很好奇人们会想出什么。

Answer 1

你提问的前提是错误的。你通常不应该投影你的数据来计算距离；投影会扭曲，因此距离不会很精确。如果您的数据具有较大的空间范围，则尤其如此。

计算最近距离的方法可能会有所不同，例如取决于您有多少个点。但是在这个例子中你可以使用蛮力计算所有的距离

library(raster)

r1 <- raster(nrow=360, ncol=720)
crs(r1) <- "+proj=longlat +datum=WGS84 +no_defs"
set.seed(1)
values(r1) <- rbinom(ncell(r1), 2, 0.01)
r1_presence <- rasterToPoints(r1, fun=function(x)x==1)
survey_points <- cbind(rnorm(50) * 5 + 10, rnorm(50) + 50)

d <- pointDistance(r1_presence[,1:2], survey_points, lonlat=TRUE)
mind <- apply(d, 1, min)

head(mind)
#[1] 4268674 4261209 4258182 4254560 4220200 4218188

使用terra v 1.2-0（目前是开发版，你可以install from github）你可以

library(terra)
#terra version 1.2.0
from <- vect(r1_presence[,1:2], crs="+proj=longlat +datum=WGS84")
to <- vect(survey_points, crs="+proj=longlat +datum=WGS84")
x <- nearest(from, to)
x
# class       : SpatVector 
# geometry    : points 
# dimensions  : 5158, 7  (geometries, attributes)
# extent      : -2.271833, 22.16701, 48.51506, 51.98087  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
# names       : from_id from_x from_y to_id  to_x  to_y  distance
# type        :   <num>  <num>  <num> <num> <num> <num>     <num>
# values      :       1 -171.2  89.75    25 8.752 51.98 4.269e+06
#                     2 -128.2  89.75    25 8.752 51.98 4.261e+06
#                     3 -119.8  89.75    25 8.752 51.98 4.258e+06

对于地理点的全球数据集（在 lat/long 中）如何找到最近的邻居来解释我们星球的球形性质

For a global dataset of geographic points (in lat/long) how to find nearest neighbors accounting for the spherical nature of our planet

r

geospatial