2017-12-05 9 views
1

私は特定の個人にそのアドレスに基づいて粒状物質の曝露を割り当てようとしている研究に取り組んでいます。私は経度と緯度の座標を持つ2つのデータセットを持っています。 1つは個人用、1つは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]) 
} 

arrows plot

+0

ありがとう、それは多少役に立ちます。私は、 "r"オブジェクトに含まれる情報から、最も近いブロックIDを持つsubjedtIDと一致するデータセットに問題が起こっていると思います。 – afossa

+0

私はそれを追加しました: 'data.frame(subject = subjectID、block = blockID [r])' – RobertH

0

大きなデータセットをお持ちの場合は、this answerの@ user3507085で説明されているように、非常に効率的な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") 
+0

Egeさん、ありがとうございました。データセットは非常に大きい。私はこのバージョンでもプレイします。地理的距離に変換する際に正しい楕円体(球形の地球モデル)を使用していることを確認しながら、何かを検討してください。これは、共通のWGS楕円体を使用するpointDistance関数の利点です。 – afossa

+0

私はこのアプローチを 'nabor'で使用して各点の最近隣を見つけて、別の関数(' pointDistance'または 'geosphere :: distGeo')を使って距離を計算することができます最も近い隣人。 –

関連する問題