2016-11-29 76 views
0

給定一個包含緯度和經度的數據幀,我想添加一個簡單包含某個半徑範圍內其他點(同一個數據幀)的點數的列,例如在特定點的10公里內。計算某個半徑內的點數

示例數據:

set.seed(1) 
radius<-10 
lat<-runif(10,-90,90) 
long<-runif(10,-180,180) 
id<-1:10 
dat<-cbind(id,lat,long) 

     id  lat   long 
[1,] 1 -42.20844 -105.8491530 
[2,] 2 -23.01770 -116.4395691 
[3,] 3 13.11361 67.3282248 
[4,] 4 73.47740 -41.7226614 
[5,] 5 -53.69725 97.1429112 
[6,] 6 71.71014 -0.8282728 
[7,] 7 80.04155 78.3426630 
[8,] 8 28.94360 177.0861941 
[9,] 9 23.24053 -43.1873354 
[10,] 10 -78.87847 99.8802797 

現在給出的半徑可變我希望有一個新的列說,「X」爲每個點只包含有不到「半徑」,其他的點數。我不關心這些是哪一點。

雖然這R - Finding closest neighboring point and number of neighbors within a given radius, coordinates lat-long主題和答案關閉它並不能解決簡單計數的具體問題。這個問題是不同的,因爲我需要的半徑內的所有點的次數,而不是點

+0

勾股定理,'distance²=(point1x - point2x)²+(point1y - point2y) ²',或者如果你想保存CPU並且不計算平方根,你可以簡單地測試'10 2 <((point 1 x - point 2 x)2 +(point 1 y - point 2 y)2)'來知道它是否在半徑內。我不知道'R'的語法,所以我無法幫助你,但我相信你可以從這裏弄明白。 – Havenard

+0

我實際上是在尋找地理距離,R語法是這裏的關鍵,特別是如何獲得該半徑內所有點的數量。 – Kathi

+0

@Havenard我們需要使用像Haversine這樣的距離度量來找到地球上兩個點之間的距離(用經緯度表示),Euclidean距離度量在這裏不起作用。 –

回答

2

試試這個:

library(geosphere) 
cbind(dat, X=rowSums(distm (dat[,3:2], 
     fun = distHaversine)/1000 <= 10000)) # number of points within distance 10000 km 

     id  lat   long X 
[1,] 1 -42.20844 -105.8491530 5 
[2,] 2 -23.01770 -116.4395691 5 
[3,] 3 13.11361 67.3282248 5 
[4,] 4 73.47740 -41.7226614 6 
[5,] 5 -53.69725 97.1429112 4 
[6,] 6 71.71014 -0.8282728 6 
[7,] 7 80.04155 78.3426630 6 
[8,] 8 28.94360 177.0861941 5 
[9,] 9 23.24053 -43.1873354 6 
[10,] 10 -78.87847 99.8802797 4 
+0

申請是非常低效的,在這種情況下 –

+0

distm的輸出是什麼指標?米或公里? – Kathi

+1

@Kathi以米爲單位。 –