R中的空间最近邻分配

Posted 2023-03-13

【中文标题】R中的空间最近邻分配

【英文标题】：Spatial nearest neighbor assignment in R

【发布时间】：2018-05-19 16:23:23

【问题描述】：

我正在进行一项研究，试图根据具体个人的地址将颗粒物暴露量分配给他们。我有两个带有经度和纬度坐标的数据集。一个 if 用于个人，一个 if 用于 pm 曝光块。我想根据最接近的块为每个主题分配一个 pm 曝光块。

library(sp)
library(raster)
library(tidyverse)

#subject level data
subjectID<-c("A1","A2","A3","A4")

subjects<-data.frame(tribble(
~lon,~lat,
-70.9821391,    42.3769511,
-61.8668537,    45.5267133,
-70.9344039,    41.6220337,
-70.7283830,    41.7123494
))

row.names(subjects)<-subjectID

#PM Block Locations 
blockID<-c("B1","B2","B3","B4","B5")

blocks<-data.frame(tribble(
~lon,~lat,
-70.9824591,    42.3769451,
-61.8664537,    45.5267453,
-70.9344539,    41.6220457,
-70.7284530,    41.7123454,
-70.7284430,    41.7193454
))

row.names(blocks)<-blockID

#Creating distance matrix
dis_matrix<-pointDistance(blocks,subjects,lonlat = TRUE)

###The above code doesnt preserve the row names. Is there a way to to do 
that?

###I'm unsure about the below code
colnames(dis_matrix)<-row.names(subjects)
row.names(dis_matrix)<-row.names(blocks)

dis_data<-data.frame(dis_matrix)

###Finding nearst neighbor and coercing to usable format 
getname <-function(x) 
row.names(dis_data[which.min(x),])


nn<-data.frame(lapply(dis_data,getname)) %>% 
gather(key=subject,value=neighbor)

此代码为我提供了有意义的输出，但我不确定其有效性和效率。任何有关如何改进和修复此代码的建议表示赞赏。我还收到错误消息：

Warning message:
attributes are not identical across measure variables;
they will be dropped 

我无法确定其来源。

感谢您的观看！

【参考方案1】：

这里有一些示例数据，说明如何使用pointDistance：

library(raster)

#subject level data
subjectID <- c("A1","A2","A3","A4")
subxy <- matrix(c(-65, 42, -60, 4.5, -70, 20, -75, 41 ), ncol=2, byrow=TRUE)
#PM Block Locations 
blockID <- c("B1","B2","B3","B4","B5")
blockxy <- matrix(c(-68, 22, -61, 25, -70, 31, -65, 11,-63, 21), ncol=2, byrow=TRUE)

# distance of all subxy to all blockxy points
d <- pointDistance(subxy, blockxy, lonlat=TRUE)

# get the blockxy record nearest to each subxy record
r <- apply(d, 1, which.min)
r
#[1] 3 4 1 3

所以这些对是：

p <- data.frame(subject=subjectID, block=blockID[r])
p

#  subject block
#1      A1    B3
#2      A2    B4
#3      A3    B1
#4      A4    B3

说明它的工作原理：

plot(rbind(blockxy, subxy), ylim=c(0,45), xlab='longitude', ylab='latitude')
points(blockxy, col="red", pch=20, cex=2)
points(subxy, col="blue", pch=20, cex=2)
text(subxy, subjectID, pos=1)
text(blockxy, blockID, pos=1)
for (i in 1:nrow(subxy)) 
    arrows(subxy[i,1], subxy[i,2], blockxy[r[i],1], blockxy[r[i],2])

【讨论】：

谢谢，这有点帮助。我认为我遇到问题的地方是从“r”对象中包含的信息到与最接近块 ID 的 subjedtID 匹配的数据集。

我补充说：data.frame(subject=subjectID, block=blockID[r])

【参考方案2】：

如果您有一个大数据集，您可能希望使用 @user3507085 在 this answer 中解释的非常高效的 nabor 包。由于该问题已作为题外话关闭，因此我已将答案复制粘贴在下面，因此它在此线程中“保持活力”。我不知道这是否被认为是不好的做法，如果需要，我很乐意删除/编辑（注意knn 给出的距离不是地理距离，但我想他们可能是通过包括 arcsin 在内的简单变换转换为球面距离）：

lonlat2xyz=function (lon, lat, r) 

lon = lon * pi/180
lat = lat * pi/180
if (missing(r)) 
    r <- 6378.1
x <- r * cos(lat) * cos(lon)
y <- r * cos(lat) * sin(lon)
z <- r * sin(lat)
return(cbind(x, y, z))


lon1=runif(100,-180,180);lon2=runif(100,-180,180);lat1=runif(100,-90,90);lat2=runif(100,-90,90)

xyz1=lonlat2xyz(lon1,lat1)
xyz2=lonlat2xyz(lon2,lat2)

library(nabor)

out=knn(data=xyz1,query = xyz2,k=20)

library(maps)

map()
points(lon1,lat1,pch=16,col="black")
points(lon2[1],lat2[1],pch=16,col="red")
points(lon1[out$nn.idx[1,]],lat1[out$nn.idx[1,]],pch=16,col="blue")

【讨论】：

感谢 Ege，在这里考虑效率肯定是有帮助的。数据集非常大。我也会玩这个版本。尽管它确保您使用正确的椭球体（地球的球形模型）来转换为地理距离，但需要考虑一些事情。这就是使用普通 WGS 椭球的 pointDistance 函数的优点。

我认为您可以将这种方法与nabor 一起使用来找到每个点的最近邻居，然后您可以使用另一个函数（如pointDistance 或geosphere::distGeo）来准确计算到最近邻居的距离。