问题描述
所以我们在不考虑投影问题的情况下完成了这项工作。问题是在哪里(以及如何)最好地添加重新投影,以便函数返回以公里为单位的值而不是当前度数:
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,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,r1_presence)
理想情况下,一个完美的答案是扩展 get_distances
函数以添加一个新选项,该选项执行适当的重新投影并以公里为单位返回值。
这里可能有几种不同的方法,我很好奇人们会想出什么。
解决方法
你问题的前提是错误的。通常不应将数据投影到计算距离;投影会失真,因此距离不会很精确。如果您的数据具有较大的空间范围,则尤其如此。
计算最近距离的方法可能会有所不同,例如取决于您拥有的点数。但在这个例子中,你可以只使用蛮力并计算所有距离
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[,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