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.

I am going crazy here, so anyone who helps me i will be ever grateful to them.

SOS

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)
  • If R crashes with a fatal error then that's probably a bug - you should try to create a minimal example and report it. It would help if we had some data we could all use to reproduce your problem. – Spacedman Aug 9 at 14:22
  • Also, please only state packages you actually use in your examples - I don't think you use tidyverse here and that breaks your script for people who don't have it installed. – Spacedman Aug 9 at 14:24
  • Great! Your code uses up_booths_1 and up_villages_1 which you don't create, but if I do pi <- intersect(up_booths_transform_buffer1km, up_villages_transform) then it doesn't crash or error, and I get the intersections that I suspect you want. Can you check that it fails for you when you intersect the buffer with the villages? Can you edit and add version numbers of R, sp, rgeos, raster? Have you tried intersecting subsets of your data - running out of memory can cause crashes. Try the first ten buffers with all the villages, for example. – Spacedman Aug 9 at 14:52
  • You'll get that error if you've not done library(raster) because without that it uses the default intersect function, not the clever geographical one from the raster package. – Spacedman Aug 9 at 15:20
  • Thank you so much for helping me, you are godsend. I am including raster library, adn have even removed all the extra ones. I still get the same error. Just now i added raster::intersect and then my R crashed with fatal error :( – Tini Tan Goya Aug 9 at 15:23
up vote 0 down vote accepted

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.

  • Thanks very much! this works great. but now i have another issue, when i run this on the whole data i get "Error in RGEOSBinPredFunc(spgeom1, spgeom2, byid, func) : rgeos_binpredfunc_prepared: maximum returned dense matrix size exceeded". I tried to use split to create a list by breaking booths and villages per district - but its super slow. I am putting this up as another question, i would be grateful if you could look at it @spacedman – Tini Tan Goya Aug 10 at 8:47
  • here is the new question - gis.stackexchange.com/questions/292388/… – Tini Tan Goya Aug 10 at 8:58
  • The new question is using sp and rgeos again. Why aren't you using sf like in my answer? – Spacedman Aug 10 at 10:02

Your Answer

 

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.