[GIS] Subset a SpatialPolygonsDataFrame by ID in R

I've used the gIntersect function from the rgeos package to determine which polygons intersect in two SpatialPolygonsDataFrames (spdfs). The result is a logical matrix with rownames comprising the polygon IDs from sfdf #1 and colnames the polygon IDs from spdf #2. I'd like to use this matrix to subset the two original spdfs, omitting those polygons that have no overlap, before using the function gIntersection. I can easily convert the logical matrix into two vectors of IDs (one for each spdf) that represent those polygons with at least some overlap.

My question is – can I use a vector of polygon IDs to subset a spdf? If so, will it correctly subset values in the @data slot? I know I can generally subset spdfs using this syntax:

spdf[spdf@data$myVariable[someVector], ]

This won't work for IDs, as they are buried much deeper inside the spdf list structure.

Best Answer

The matrix ids returned from gIntersect should correspond to the rownames in each source sp object. You should be able to just index the rownames position in order to subset the data.

r <- c(1,5,3,9,10)
sp.polys <- sp.polys[r,]

Related Solutions

[GIS] Intersection between two SpatialPolygonsDataFrame using gIntersection function in R

I cannot as yet offer a complete solution, but the immediate issue is that the parcel object has invalid geometry. If you can fix the geometry, we should be able to get the intersections to work. As an example of problematic geometry, this expression is fine:

gIsValid(hr_kud, reason = TRUE)

But this one will consistently crash R:

gIsValid(parcel, reason = TRUE)

So, one or more of the 20,069 polygons in parcel appears to be broken in some way. I don't know how to validate geometry within R but Quantum GIS does have a tool for this. I assume that you have the original shapefile or other data source from which you created parcel but you can recreate this by using the writeOGR function from the rgdal library like so:

library(rgeos)
library(rgdal)

setwd("your_directory_name_here")

load("Data.RData")
class(parcel)
class(hr_kud)
object.size(hr_kud)
object.size(parcel)
parcel@proj4string
hr_kud@proj4string

gIsValid(hr_kud, reason = TRUE)
#gIsValid(parcel, reason = TRUE) # warning! this will crash R

writeOGR(parcel, ".", driver = "ESRI Shapefile", layer = "poly")

This should produce 4 files:

poly.shp
poly.shx
poly.dbf
poly.prj

In QGIS, first load your shapefile by going to the 'Layer' menu and clicking on 'Add Vector Layer..." then select the poly.shp file or the folder in which the file resides. Then go to the 'Vector' menu, click on the 'Geometry Tools' menu item then on the 'Check geometry validity' item. Then click OK on the next dialog box. With more than 20k polys this took quite some time (about half an hour I think), but the result was this: [EDIT: I upgraded to QGIS 1.8 Lisboa with two results - a massive increase in validation speed and the discovery of other errors, as per the revised screenshot below.]

QGIS geometry check result

So, it looks like you have a number of issues there. Before you close the dialog box, click on the error in the dialog box list and QGIS will show you the offending polygon. I'm not experienced in this kind of GIS but the tips in this tutorial might help. Google QGIS and "self intersecting" for more ideas.

[GIS] R apply function in turn to each subset of a spatialPolygonsDataFrame

apply is not really the right tool, since this is a data frame (of sorts) not an array.

How about

gArea(parcProj, byid = TRUE)

Best Answer

Related Solutions

[GIS] Intersection between two SpatialPolygonsDataFrame using gIntersection function in R

[GIS] R apply function in turn to each subset of a spatialPolygonsDataFrame

Related Question