r – How to Select Polygons by Points in R Using Overlay Techniques

overlayr

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)

Best Answer

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)

Related Solutions

[GIS] Checking if Lng and Lat fall inside polygons from ESRI Shapefile

You will want to import your shape file. You can try package rgdal, specifically function readOGR.

After reading into R your csv file, you may want to coerce it to a SpatialPoints class. Make sure your projections are identical.

After you have the polygon and SpatialPoints, it should be just a matter of using function over from package sp.

I recently asked a similar question. I hope the data in the self contained example should be sufficient to work through the second part of your question.

R – Calculate Distance Between Points and Nearest Polygon

Let's set up some sample data:

library(sp)
library(spdep)
example(columbus)
plot(columbus)

We'll now get some points. Click 20 times on the map, some inside and some outside polygons:

pts = locator(20,type="p")

Convert to Spatial data type:

spts = SpatialPoints(pts)

Now use rgeos to compute point-polygon distances, and take the minimum for each point:

library(rgeos)
> apply(gDistance(spts, columbus,byid=TRUE),2,min)
         1          2          3          4          5          6          7 
0.25947439 0.03898401 1.27156515 1.50490316 0.00000000 0.00000000 0.00000000 
         8          9         10         11         12         13         14 
0.00000000 0.00000000 0.31312329 0.26742466 0.53934325 0.07764322 0.03909773 
        15         16         17         18         19         20 
0.11156343 0.29243322 0.08872334 0.00000000 0.00000000 0.00000000

There you go. That's the minimum distance from each of the 20 points to any of the polygons. The zeroes are for points inside one of the polygons.

Note you should use data in a projected spatial coordinate system and not lat-long.

If you have 17,000 points then you might want to do this in smaller subsets to avoid creating a 17,000 x 5,000 matrix if you have 5,000 polygons. You didn't say.

Best Answer

Related Solutions

[GIS] Checking if Lng and Lat fall inside polygons from ESRI Shapefile

R – Calculate Distance Between Points and Nearest Polygon

Related Question