2017-12-20 20 views
0

私は郡のパーセルレベルのシェープファイルを持っています。私は同じ所有者だけでなく、マイル(約1610メートル)内のパーセルの数を計算することを目指しています。私は解決策を試してきましたが、以下は私のサンプルコードですが、かなり非効率的で、メモリ集約的です。私は公にデータをポストすることはできませんが、ここでいくつか作ったコードに問題がされています大規模なデータセットを持つ半径内のポイント数 - R

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は数十億の要素を持ち、おそらく私より多くのメモリを必要とします)。 スピードとメモリの両方の観点から、誰もこの問題に近づく良い方法についての考えを持っていますか?ヘルプは非常に高く評価されます。

+0

サンプルデータを生成した場合は、非常に役立ちます。参照:https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example/5965451#5965451 – RobertH

+0

@RobertH、ありがとう、私はそうした – winitheju

答えて

0

あなたは、私は、彼らが複数の荷物を持っている場合、それは他の所有者を複数回含まれてnumdifのあなたの計算は、正しくないと思わ適用

d <- d < 1609.34 
nt <- apply(d, 1, sum) 
nd <- apply(d, 1, function(i) length(unique(owner[i]))) - 1 

でカウントを行っている可能性があります。観測の多数を考えると

、私はこのルートを検討します:

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 

遅いですが、それは狂気の大行列を作成しません

私はまた、あなたのアプローチを前提としていることを追加する必要があります小包は考慮された距離に対して相対的に小さいので、重心を使用することはOKである。 これが当てはまらない場合は、計算コストを増やして、rgeos::gWithinDistanceでポリゴンの計算を行うことができます。

+0

ロバート、私は実際に複数の所有者が複数回OKですので、問題はありません。私は今旅行していますが、これを試し、次の日にどのように進むのかを教えてくれるでしょう。ありがとうございます – winitheju

+0

ロバート、これは働いた。しかし、それははるかに遅いです。 10,000ポイントは約2分、35,000ポイントは約25分かかりました!しかしこれは避けられないかもしれません。ご協力いただきありがとうございます – winitheju

関連する問題