r - Creating buffers around points and merging with SpatialPolygonsDataFrame to create a list of matched polygons and area of intersection?

I have data on polling booths as points and data on village boundaries as polygons. I want to create buffers of 1km around the polling booths and (a) save the ID of all the villages that intersect with the buffer (b) in another column have the amount of area that intersects in Km2.

I have tried everything but it doesn't work. I am sharing my code below. When i use intersect on Spatialpolygons my R crashes with fatal error. I tried to rasterise but i get some strange output that i do not understand.

sample dataset - https://drive.google.com/open?id=1Ha47ZqwG7bCuxExJk6zOs8KoqDiqD5A4

rm(list=ls())
##LOADING DATASETS REQUIRED FOR THIS
load("sampledata.Rdata")

#booths data - dataframe - upbooths_raw
#booths data - spatial file - up_booths_sample
#village data - up_villages_sample
#loading libraries that i might use
library(rgeos)
library(raster)



up_booths<-up_booths_raw
#up_booths_raw<-subset(upid2017matched_gis2,district_name_17=="Jhansi")
up_booths<-subset(up_booths, is.na(latitude.x)==FALSE & is.na(longitude.x)==FALSE)

#project the booths
coordinates(up_booths) <- ~longitude.x+latitude.x

#make same coordinates copy village CRS to booths

up_booths@proj4string <- up_villages_sample@proj4string

#new CRS which is good for Uttar Pradesh India
newcrs <- CRS("+proj=utm +zone=42N +datum=WGS84 +units=km")

#transforms booths and villages to new crs
up_booths_transform <- spTransform(up_booths, newcrs)
up_villages_transform<- spTransform(up_villages_sample,newcrs )

#create booths with 1km buffer

up_booths_transform_buffer1km <- gBuffer(up_booths_transform, byid = TRUE, width = 1, capStyle="ROUND")

#runinng this crashes R into fatal error
pi <- intersect(up_booths_transform_buffer1km, up_villages_transform)

Answer

It looks like something is broken in your sp or raster or possibly rgeos setup. Upgrade everything to the latest versions and try again.

Or use the sf packages for your data objects. Here is a script that loads your sample data and converts to sf and does the transform, buffer, area calculation and intersection:

First use these packages and load the data:

library(sf)
library(sp)

load("sampledata.Rdata")

Convert data frame to sf using these as coordinates:

booths = st_as_sf(up_booths_raw,coords=c("longitude.y","latitude.y"))

Convert sp to sf:

villages = st_as_sf(up_villages_sample)

villages has a coordinate system, take it to the booths:

st_crs(booths) = st_crs(villages)

At this point you should plot the villages and the booths to make sure nothing has gone wrong. Then proceed:

Transform to metric coordinate system:

booths_km = st_transform(booths, "+proj=utm +zone=42N +datum=WGS84 +units=km")
villages_km = st_transform(villages, "+proj=utm +zone=42N +datum=WGS84 +units=km")

Create 1km buffer zone:

booths_buffer = st_buffer(booths_km, 1)
booths_villages = st_intersection(booths_buffer, villages_km)

And that's our intersections. How many?

dim(booths_villages)
## [1] 3034   41

3034 - more than the number of booths and villages, but that's okay because one booth buffer can overlap multiple villages.

Compute the area of each intersection:

booths_villages$area = st_area(booths_villages)

And have a look at some columns:

booths_villages[1:4, c("booth_id_12","NAME","area")]


>     booths_villages[1:4, c("booth_id_12","NAME","area")]
Simple feature collection with 4 features and 3 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 1502.222 ymin: 2904.236 xmax: 1504.222 ymax: 2906.236
epsg (SRID):    32642
proj4string:    +proj=utm +zone=42 +datum=WGS84 +units=km +no_defs
      booth_id_12     NAME          area                       geometry
60586           8 Lohagarh 3.140157 km^2 POLYGON ((1504.222 2905.236...

60587           9 Lohagarh 3.140157 km^2 POLYGON ((1504.222 2905.236...
60588          10 Lohagarh 3.140157 km^2 POLYGON ((1504.222 2905.236...
60589          11 Lohagarh 3.140157 km^2 POLYGON ((1504.222 2905.236...

That's an id from the booth tab, a name from the village table, the area of intersection and the polygon of intersection. I couldn't find a unique identifier in the booths table so you might want to add a row number to booths and do it again.

If this fails with a crash then there's possibly something wrong with your installation of the underlying geographic processing libraries.