Create buffer and count points in R

Question

I asked this question previously but did not get a response, so I'll try and do a better job this time around!

I want to analyze spatial density of gas station points using R. I need to create a buffer (let's say 1,000m) around the gas stations and count the number of gas stations within the buffer. I'll then need to play around with buffer distances to see what's a reasonable buffer to see something interesting. I won't post the entire shape file because it's fairly messy, but this is what the data look like:

all <- readShapePoints("sbc_gas.shp") 
all.df <- as(all, "data.frame")
head(all)

OBJECTID Fuellocati               Name   Latitude      Longitude     
      1      34828     WORLD OIL #104    34.44190      -119.8304    
      2      48734  STOP AND SHOP GAS    34.41962      -119.6768    
      3      51276 EL RANCHERO MARKET    34.41911      -119.7162    
      4      52882  EDUCATED CAR WASH    34.44017      -119.7439    
      5      74038           CIRCLE K    34.63925      -120.4406    
      6     103685    7-ELEVEN #23855    34.40506      -119.5296

I was able to create a buffer around the points with the following code, but now how do I count the number of points within the buffer?

require(sp)
require(rdgal)
require(geosphere)

coordinates(all) <- c("Longitude", "Latitude")
pc <- spTransform(all, CRS( "+init=epsg:3347" ) ) 
distInMeters <- 1000
pc100km <- gBuffer(pc, width=100*distInMeters, byid=TRUE )
# Add data, and write to shapefile
pc100km <- SpatialPolygonsDataFrame(pc100km, data=pc100km@data )
writeOGR( pc100km, "pc100km", "pc100km", driver="ESRI Shapefile" )

plot(pc100km)

I'm open to other ways to go about this.

Can you give more info about "can't get either to work?" Are you getting any error output? Can you enter debug mode and inspect what the variables are holding? — Blake Regalia
For starters, it looks like distm isn't defined. Did you mean to put distInMeters there instead? — Blake Regalia
distm should be okay. It's a function from the geosphere package that takes the lat & long of points and creates a matrix of the distance between the points. I'll play around with this though. Maybe it's the issue. — JAG2024
I think if you're trying to "count the number of points within the buffer" what you really need to do is create a nested for loop and compute the distance between every pairwise combination, then set a condition when the distance is less than the buffer radius. What you have now just creates buffer geometries about the input points. — Blake Regalia

JAG2024 JAG2024 · Accepted Answer · 2017-03-01T09:03:04

Came up with a (simple) solution for my own question using geosphere:

cbind(coordinates(all.df), X=rowSums(distm (coordinates(all.df)[,1:2], fun = distHaversine) / 1000 <= 10)) # number of points within distance 10 km
    Longitude Latitude  X
0   -119.8304 34.44190 25
1   -119.6768 34.41962 29
2   -119.7162 34.41911 34
3   -119.7439 34.44017 39
4   -120.4406 34.63925 13
5   -119.5296 34.40506  7
6   -120.4198 34.93860 26
7   -119.8221 34.43598 30

Create buffer and count points in R

1 Answers