4

Given a data frame with latitudes and longitudes I want to add a column that simply contains the count of other points (of the same data frame) that are within a certain radius e.g. within 10 km of that specific point.

Example data:

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

Now given the radius variable I want a new column say "X" that for each point only contains the number of other points that are within "radius". I do not care about which points these are.

While this R - Finding closest neighboring point and number of neighbors within a given radius, coordinates lat-long topic and answer get close it does not solve the specific question of the simple count. This question is different as I need the counts of all points within the radius and not the points

Community
  • 1
  • 1
Kathi
  • 99
  • 1
  • 5
  • Pythagorean theorem, `distance² = (point1x - point2x)² + (point1y - point2y)²`, or if you want to save CPU and not calculate square root you can simply test if `10² < ((point1x - point2x)² + (point1y - point2y)²)` to know if its inside the radius. I don't know `R` syntax so I cant help you there, but I'm sure you can figure out from here. – Havenard Nov 29 '16 at 20:04
  • 1
    I am actually looking for geo distance and the R syntax is the key here and in particular how to get the count of all points within that radius. – Kathi Nov 29 '16 at 20:06
  • @Havenard we need to use some distance metric like Haversine to find distance between two points on earth (expressed by lat / lon), Euclidean distance metric will not work here. – Sandipan Dey Nov 29 '16 at 20:16

1 Answers1

5

Try this:

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
Pierre L
  • 28,203
  • 6
  • 47
  • 69
Sandipan Dey
  • 21,482
  • 2
  • 51
  • 63