2017-12-05 19 views
1

나는 그들의 주소를 기반으로 특정 개인에게 미립자 물질 노출을 할당하려고하는 연구에 착수했습니다. 경도와 위도 좌표가있는 두 개의 데이터 세트가 있습니다. 하나는 개인의 경우이고 다른 하나는 오후 노출 블록의 경우입니다. 가장 가까운 블록을 기반으로 각 피사체에 오후 노출 블록을 할당하려고합니다. R의 공간 가장 가까운 이웃 할당

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'와 같이) 사용하여 거리를 계산할 수 있습니다. 가장 가까운 이웃. –