The behaviour of gIntersection is not to pass any intersected data by design:
Since there are no general matches between intersected spatial
objects, any arbitrary operations on attributes require assumptions
about unknown user intentions. This is why no data slots should be
passed through ...
... The design of gIntesection() is inentional, because only the user can
know what to do with attributes of entities that have their geometries
changed. Different users may make different assumptions, but there is
no general solution beyond passing through the IDs of the intersecting
geometries, as is done in the row.names() mechanism.
To my surprise, the raster package has a intersection function, which simply intersects and hands over the data as well.
The raster package has a few functions that extend rgeos by also
attempting to handle attribute data as well. In this case, see
raster::intersect And the list of functions here: ?"raster-package"
(section XIV)
The complete info I got on this:
http://r-sig-geo.2731867.n2.nabble.com/Intended-usage-of-gIntersection-td7587120.html
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.]
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.
Best Answer
There is no real difference other than the class of the returned object. The raster intersect function is a helper function that, for polygons, calls gIntersection from rgeos (not rgdal). I would recommend using raster's intersect functions because it will save you some steps in getting back to a SpatialPolygonsDataFrame object. One good way to explore these types of questions is to download the package source code and look at how it is structured and what it is doing.