11

I am trying to intersect two SpatialPolygonsDataFrames and get a SpatialPolygonsDataFrame as the result. Unfortunately, using the gIntersection function from rgeos (which works impressively quickly to intersect the polygons), I cannot seem to retrieve the associated dataframes. Consider the following example:

> fracPoly <- gIntersection( toSingle, fromSingle )
> class(toSingle)
[1] "SpatialPolygonsDataFrame"
> class(fromSingle)
[1] "SpatialPolygonsDataFrame"
> class(fracPoly)
[1] "SpatialPolygons"

I can write a wrapper function which handles the transfer of data.frames, but it will be a minor pain to get all the checking right and before I did I was hoping someone could either confirm that there's no better way or point me towards another function (or option for gIntersection) which would allow me to retain the associated data.frames.

Update

On further reflection, this may may be very deliberate behavior by gIntersection. After all, of the two SPDFs, whose data.frame do you pass along? So I may have to write a wrapper which merges the two.

Ari B. Friedman
  • 2,234
  • 5
  • 21
  • 25
  • 1
    How are you loading your polgyons - readShape* or readOGR? I get weird behaviour with gIntersection depending on which library I've used to load SPDFs, which I've not gotten to the bottom of. – Simbamangu Oct 06 '12 at 07:27
  • @Simbamangu Interesting. I'm using readShapePoly and then merging in a data.frame.... – Ari B. Friedman Oct 06 '12 at 10:01
  • gIntersection SHOULD give the data.frames merged for the overlapping areas - if I run Vector|Geoprocessing|Intersect in QGIS, the output is a set of merged attributes for the overlap, and doesn't QGIS use the same GEOS library for spatial operations like this? – Simbamangu Oct 15 '12 at 07:30
  • Well if it should then I'm stumped. I already rewrote my code to merge the data.frames by hand and munge them back in, so I'm not going to spend any more time on this for now. But the hint that the function reading in the shapefile matters is helpful. – Ari B. Friedman Oct 15 '12 at 13:55
  • Actually after updating my packages both functions act the same (was getting odd errors with readShape* loaded objects before); annoying that the QGIS code is seemingly more advanced than the rgeos. Could really use it in R, though - any chance you could share that code? – Simbamangu Oct 15 '12 at 14:30
  • @Simbamangu As I recall I wrote code specific to the one variable I needed rather than merging both. But I'll take a look and see if I can generalize it. Should just be a merge of the @data slots, with the appropriate all. option. – Ari B. Friedman Oct 15 '12 at 16:04
  • 1
    This clearly is not an answer, but I don't have enough points to leave a comment... Ari, I was wondering if you would share your chunk of code to extract the variables following rgeos overlay functions. I am having trouble wrapping my head around a good way to retain the original polygon id's from various rgeos operations such as gUnion... – jed.a.long Feb 28 '13 at 21:55
  • @jed.a.long I'll try to go find the code I wrote and release it soon. Thanks for reminding me. – Ari B. Friedman Mar 01 '13 at 02:07
  • @jed.a.long It's been awhile, but I think the code you asked for is here: https://github.com/gsk3/taRifx.geo/blob/2fc3fdd35713a028bef8eefa4bef9fdc446b1b4c/R/R-GIS.R#L518 – Ari B. Friedman Mar 08 '13 at 02:10

3 Answers3

10

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

Bernd V.
  • 3,179
  • 25
  • 49
2

For some project i had the same need. Much more than keeping the data.frame, we had to put in place further code to manage output geometry type, and proceed to some cleaning (e.g. clean geometry collections), to have some complete intersection geoprocess. In case you still need to do such a task in R, you can try the RFigisGeo package:

#install RFigisGeo
require(devtools) 
install_github("RFigisGeo", "openfigis")
require(RFigisGeo)

#compute intersection
result <- getIntersection(features1, features2)
eblondel
  • 921
  • 6
  • 9
0

For those (like me) for whom the above answers did not work, the link here explains that you can do this precise thing with raster's intersect.

How do I retain all attribute data when clipping two polygons in R?

I used this to crop a SpatialPointsDataFrame with a SpatialPolygons shapefile; it creates a cropped/clipped version of the SpatialPointsDataFrame, maintaining the original data.

Leah Bevis
  • 195
  • 2
  • 7