对于表中的每个观测值,根据纬度和经度 (R) 计算 x 米内的表中其他观测值的数量
Posted
技术标签:
【中文标题】对于表中的每个观测值,根据纬度和经度 (R) 计算 x 米内的表中其他观测值的数量【英文标题】:For each observation in a table, count number of other observations in table within x metres based on latitude and longitude (R) 【发布时间】:2020-06-19 06:44:25 【问题描述】:我有一个带有经度和纬度的位置数据框,它看起来像:
set.seed(211)
latitude <-runif(5000, min=50, max=55)
longitude <- runif(5000, min=-2, max=0)
location_id <- seq(1,5000)
reprex <- data.frame(location_id, latitude, longitude)
对于每个 location_id,我需要计算列表中距离该位置 10 英里(约 16000 米)范围内的其他位置的数量。
我正在考虑在某种 for 循环(或者可能是应用函数)中为此使用 geosphere::distGeo(),但我只是无法弄清楚如何对其进行编码,以便它比较列表中的所有其他项目使用当前项目并计算有多少在某个阈值内,记录该值,然后移动到下一行。
有人知道怎么写吗?
【问题讨论】:
这里的计算速度是一个重要问题吗?针对您的问题,最简单的代码选项可能不会是最快的。如果您需要速度,当然值得开发一种算法,将您的点分为几类(想想地图上的网格模式)。如果您的网格步长为 10 英里,那么您只需查看同一组或相邻组中的元素,而不是浏览每个点的整个地图。根据您拥有的点数及其密度,此优化的影响可能或多或少重要。 【参考方案1】:distGeo
函数可以做到这一点,但你需要一个循环。注意坐标的第一列必须是经度。
lst <- vector(50, mode="list")
for(i in 1:50)
dist <- distGeo(p1=reprex[i,c(3,2)], p2=reprex[-i,c(3,2)])
lst[[i]] <- sum(dist<16000)
reprex$n <- unlist(lst)
table(unlist(lst))
0 1 2
34 10 6
所以 50 个点中有 34 个在 10 英里(约 16,000 米)内没有任何其他点,10 个只有 1 个,6 个有 2 个。
【讨论】:
【参考方案2】:fields
中的 rdist.earth
函数似乎对此很有用,例如:
library(fields)
dist.matrix <- rdist.earth(reprex[-1])
colSums(dist.matrix<10)
在这种情况下,唯一的技巧是在逻辑矩阵上使用colSums
来计算TRUE
值的数量。
注意,miles 是默认值,km 可以与参数 miles=FALSE
一起使用。
【讨论】:
【参考方案3】:将循环隐藏在(仍然很慢)apply
中并解开纬度和经度(它们通常是相反的),您可以尝试类似
set.seed(211)
latitude <-runif(5000, min=50, max=55)
longitude <- runif(5000, min=-2, max=0)
location_id <- seq(1, 5000)
reprex <- data.frame(location_id, latitude, longitude)
library(geosphere)
within10m <- function(p1, p2, dist=16000)
sum(geosphere::distGeo(p1, p2) <= dist)
matpoints <- as.matrix(reprex[, 3:2])
reprex$neighbours <-
apply(matpoints, 1, within10m, p2=matpoints) - 1
head(reprex)
# location_id latitude longitude neighbours
# 1 1 51.17399 -1.1489713 48
# 2 2 54.52623 -1.8554624 39
# 3 3 54.84852 -0.3014742 56
# 4 4 51.72104 -1.8644226 50
# 5 5 51.32793 -0.7417923 56
# 6 6 50.07346 -0.8939857 36
【讨论】:
【参考方案4】:最终我在这里使用了答案,因为它非常优雅并且避免了循环:Calculating number of points within a certain radius
我使用了代码:
library(geosphere) # for distHaversine() and distm() functions
reprex <- cbind(reprex, # appends to the dataset...
count_nearby=rowSums( # ... a column which counts the rows in the dataset...
distm (reprex[,3:2], fun = distHaversine) # ... where the distance between current and other rows...
<= 16000)-1 # ... is less than 16000 metres. Take one away because it counts itself!
) #close the cbind brackets!
【讨论】:
以上是关于对于表中的每个观测值,根据纬度和经度 (R) 计算 x 米内的表中其他观测值的数量的主要内容,如果未能解决你的问题,请参考以下文章