具有大型数据集的半径内的点数- R
Number of points within a radius with large datasets- R
我有县地块级别的 shapefile,我的目标是计算一英里(约 1610 米)内的地块总数,以及同一所有者的地块数量。我已经完成了一个解决方案,下面是我的示例代码,但它效率很低并且占用大量内存。我不能公开 post 数据,但这是一些编造代码的问题:
library(rgdal)
library(rgeos)
library(geosphere)
nobs<-1000 # number of observations
nowners<-50 # number of different owners
crs<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
long<-runif(nobs,min=-94.70073, max=-94.24141) #roughly adair county in iowa
lat<-runif(nobs,min=41.15712,max=41.50415) #roughly adair county in iowa
coords<-cbind(long,lat)
owner<-sample(1:nowners,nobs, replace=T) # give id's to owners
df<-as.data.frame(owner)
centroids<-SpatialPointsDataFrame(coords,df,proj4string = CRS(crs)) # make spatial dataframe
d<-distm(centroids) # distance from centroids to other centroid
numdif<-matrix(0,length(owner)) #vectors of 0s to be replaced in loop
numtot<-matrix(0,length(owner))
for (i in 1:length(owner)) {
same_id<-df$owner[i]==owner ## identify locations with same owner ID
numdif[i]<-as.numeric(sum(d[i,]<1609.34 & same_id==F)) #different parcel owners
numtot[i]<-as.numeric(sum(d[i,]<1609.34)) #total parcels
}
生成的 "numdif" 和 "numtot" 向量给出了我想要的:分别具有不同所有者和总数的相邻地块数量的向量。然而,这个过程对于我的 "nobs." 大得多的县来说是非常耗时和内存密集型的一些县有 50-75,000 个观测值(因此生成的矩阵 m 有数十亿个元素,可能需要比我更多的内存有)。
从速度和内存的角度来看,有没有人想过更好的方法来解决这个问题?非常感谢帮助。
您可以在申请中完成计数
d <- d < 1609.34
nt <- apply(d, 1, sum)
nd <- apply(d, 1, function(i) length(unique(owner[i]))) - 1
我认为您对 numdif 的计算不正确,因为如果其他所有者有多个地块,它会多次包含其他所有者。
鉴于大量观察,我会考虑这条路线:
d <- lapply(1:nrow(coords), function(i) which(distGeo(coords[i, ,drop=FALSE], coords) < 1609.34))
ntot <- sapply(d, length)
ndif <- sapply(d, function(i) length(unique(owner[i]))) - 1
速度较慢,但不会创建疯狂的大矩阵
我还应该补充一点,您的方法假设包裹相对于所考虑的距离较小,因此使用质心就可以了。
如果不是这种情况,可以使用 rgeos::gWithinDistance
对多边形进行计算,但计算成本会增加。
我有县地块级别的 shapefile,我的目标是计算一英里(约 1610 米)内的地块总数,以及同一所有者的地块数量。我已经完成了一个解决方案,下面是我的示例代码,但它效率很低并且占用大量内存。我不能公开 post 数据,但这是一些编造代码的问题:
library(rgdal)
library(rgeos)
library(geosphere)
nobs<-1000 # number of observations
nowners<-50 # number of different owners
crs<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
long<-runif(nobs,min=-94.70073, max=-94.24141) #roughly adair county in iowa
lat<-runif(nobs,min=41.15712,max=41.50415) #roughly adair county in iowa
coords<-cbind(long,lat)
owner<-sample(1:nowners,nobs, replace=T) # give id's to owners
df<-as.data.frame(owner)
centroids<-SpatialPointsDataFrame(coords,df,proj4string = CRS(crs)) # make spatial dataframe
d<-distm(centroids) # distance from centroids to other centroid
numdif<-matrix(0,length(owner)) #vectors of 0s to be replaced in loop
numtot<-matrix(0,length(owner))
for (i in 1:length(owner)) {
same_id<-df$owner[i]==owner ## identify locations with same owner ID
numdif[i]<-as.numeric(sum(d[i,]<1609.34 & same_id==F)) #different parcel owners
numtot[i]<-as.numeric(sum(d[i,]<1609.34)) #total parcels
}
生成的 "numdif" 和 "numtot" 向量给出了我想要的:分别具有不同所有者和总数的相邻地块数量的向量。然而,这个过程对于我的 "nobs." 大得多的县来说是非常耗时和内存密集型的一些县有 50-75,000 个观测值(因此生成的矩阵 m 有数十亿个元素,可能需要比我更多的内存有)。 从速度和内存的角度来看,有没有人想过更好的方法来解决这个问题?非常感谢帮助。
您可以在申请中完成计数
d <- d < 1609.34
nt <- apply(d, 1, sum)
nd <- apply(d, 1, function(i) length(unique(owner[i]))) - 1
我认为您对 numdif 的计算不正确,因为如果其他所有者有多个地块,它会多次包含其他所有者。
鉴于大量观察,我会考虑这条路线:
d <- lapply(1:nrow(coords), function(i) which(distGeo(coords[i, ,drop=FALSE], coords) < 1609.34))
ntot <- sapply(d, length)
ndif <- sapply(d, function(i) length(unique(owner[i]))) - 1
速度较慢,但不会创建疯狂的大矩阵
我还应该补充一点,您的方法假设包裹相对于所考虑的距离较小,因此使用质心就可以了。
如果不是这种情况,可以使用 rgeos::gWithinDistance
对多边形进行计算,但计算成本会增加。