4

I have a shapefile with polygons and another one with points. I would like to create a new shapefile with those polygons that don't contain any point inside. I tried the function sp::over in R and it returns a vector indicating which polygons have points, but I have been unable to create the new shapefile from there.

 library(sp) 
 polygons <-over(basins,points_proj)
7

Think of the over function returning a query, now you have to do something with it. One easy approach is to pass an is.na() evaluation to a bracket index of the the SpatialPolygonsDataFrame object. This will return a subset of polygons that do not intersect the points.

First, create some example data.

require(sp)
data(meuse)
coordinates(meuse) = ~x+y
sr1=Polygons(list(Polygon(cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409, 
  180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676, 
  332618, 332413, 332349)))),'1')
sr2=Polygons(list(Polygon(cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437, 
  179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683, 
  331133, 331623, 332152, 332357, 332373)))),'2')
sr3=Polygons(list(Polygon(cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875, 
  179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
  c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004, 
  329783, 329665, 329720, 329933, 330478, 331062, 331086)))),'3')
sr4=Polygons(list(Polygon(cbind(c(180304, 180403,179632,179420,180304),
  c(332791, 333204, 333635, 333058, 332791)))),'4')
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf=SpatialPolygonsDataFrame(sr, data.frame(row.names=c('1','2','3','4'), PIDS=1:4))

plot(srdf)
  plot(meuse, pch=20, add=TRUE)

Now, we can subset the polygons to ones not intersecting points.

polys.sub <- srdf[is.na(sp::over(srdf, sp::geometry(meuse))), ] 
plot(meuse, pch=20)
  plot(polys.sub, add=T)

By adding a "is not" evaluation (using !), you can invert the subset to include polygons that DO intersect the points.

polys.sub <- srdf[!is.na(sp::over(srdf, sp::geometry(meuse))), ] 
plot(meuse, pch=20)
  plot(polys.sub, add=T)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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