0

我有一个 csv 文件，其中包含 17,305 个池塘的池塘区域和纬度和经度坐标。对于每个池塘，我想确定其 1 公里范围内所有池塘的坐标。我是 R 新手，所以我想我可以调整一些最近邻代码。我在 Crawley 的 The R Book 中找到了这个循环：

x<-runif(100)
y<-runif(100)

par(pty="s")
plot(x,y,pch=16)

distance<-function(x1, y1, x2, y2) sqrt((x2 − x1)^2 + (y2 − y1)^2)

r<-numeric(100)
nn<-numeric(100)
d<-numeric(100)
for (i in 1:100) {
for (k in 1:100) d[k]<-distance(x[i],y[i],x[k],y[k])
r[i]<-min(d[-i])
nn[i]<-which(d==min(d[-i]))
}   

for (i in 1:100) lines(c(x[i],x[nn[i]]),c(y[i],y[nn[i]]))

我对其进行了调整并在化石中使用了 deg.dist 函数，该函数使用 Haversine 公式而不是使用毕达哥拉斯。

install.packages("fossil")
library(fossil)

Pond_A<-read.csv("C:\\ PondArea_data\\Pond_areas.csv")

r<-numeric(17305)
nn<-numeric(17305)
d<-numeric(17305)
for (i in 1:17305){
for (k in 1:17305) d[k]<-with(Pond_A,deg.dist(Longitude[i],Latitude[i],Longitude[k],Latitude[k]))
  r[i]<-min(d[-i])
  nn<-which(d<=1)
}

这似乎给了我最后一个池塘 1 公里内所有池塘的身份。但是尽我所能，我无法弄清楚如何为所有池塘找到答案。如果有人能给我一个解决方案并解释它为什么起作用，我将不胜感激。

谢谢，

艾丹

4

1 回答 1

0

您可以使用 rgeos 包中的 gWithinDistance 创建一个布尔矩阵。row/col 值表示 sp 对象的行名。然后，您可以将矩阵强制转换为数据框并分配回 sp 对象。对于这个例子，我使用 sp 包中的 meuse 数据。 

require(sp)
require(rgeos)
data(meuse)
  coordinates(meuse) <- ~x+y

# Create boolean matrix where TRUE is distance condition is |nnd <= d| TRUE else FALSE
d=200
DistMat <- gWithinDistance(meuse, meuse, dist=d, byid=TRUE)  

# Turn self-evaluation values to NA 
diag(DistMat) <- NA

# Join back to data
cids <- colnames(DistMat)
  DistMat <- as.data.frame(DistMat)
    names(DistMat) <- paste("NID", cids, sep=".")
      meuse@data <- data.frame(meuse@data, DistMat) 
        str(meuse@data)
于 2013-08-28T23:32:52.270 回答