0
votes

I have the following puzzle regarding spatial data in R:

I have a dataset with street segments (with their respective starting and ending coordinates). I want to create a buffer of X meters around these points and then check whether a list of lat/lon points is within that buffer. Is there a way to do this in R?

I'm able to map the points and map the buffers using a combination of various packages: maptools, ggmap, rgdal, sp, and rgeos. But this procedure seems to only map the points and the buffers without allowing me to then check if the other coordinates I have are within the buffers. Ideally, I'd like to produce a vector of 1s and 0s describing whether the list of lat/lon points are within the buffer around the street segments.

Any ideas?

This is the code I've been using, but I get all missing values (and I know this shouldn't be the case). I've also tried using the gContains function from rgeos but it crashed R.

#Load shapefile in R and transform to appropriate CRS
shp <- readOGR(dsn="/Users/Maps/shapes", layer="shp")
shp_transf <- spTransform(shp_transf, CRS( "+init=epsg:21897" ))

#Create buffer around polygons
shp_buff <- gBuffer(shp_transf, width=40, byid=TRUE, quadsegs=10)

#Make my dataframe of lat/lon points into same projection as buffers
points <- SpatialPoints(points,proj4string=CRS(proj4string(shp_buff)))

#Use over function from SP pacakge
result <- as.integer(over(points, shp_buff)$OBJECTID)
1
What happens when you just run result <- over(points, shp_buff)?Jacob F
It creates an empty dataframe of the same dimensions as shp_buffop4
When you are declaring your points object, are the coordinates in the data frame already in the epsg:21897 CRS? The wording of the comment above that line kind of sounds like you are trying to transform them. Also, what happens when you run plot(shp_buff) then plot(points, add=T)?Jacob F
Hi Jacob - I run that to make sure that the points vector are in the same projection as the shapefile shp. When running plot(shp_buff) it creates a plot of the buffers but running plot(points, add=T) creates a null value.op4
My question is, are the points originally in the coordinate system you are defining them with or are you expecting them to be transformed by the SparialPoints function?Jacob F

1 Answers

1
votes

Turns out, the over() function wasn't working before because the shapefile with the buffers and the list of GPS coordinates were in a different CRS. I fixed it this way:

shp <- spTransform(shp, CRS("+proj=longlat +datum=WGS84")) proj4string(points) <- CRS("+proj=longlat +datum=WGS84")

Then, you run the over() function and it'll give you the number of times each point is in each buffer area:

x <- over(points, shp) table(x$ID)

181 304

118 8

So a particular GPS pulse came from segment ID 181 a total of 118 times and from segment ID 304 a total of 8 times.