2

I have two shapefiles—let’s call them shp1.shp and shp2.shp—and I want to create a new shapefile of overlaps between the two.

Essentially, I'm trying to determine which geographies from shp1 fall within shp2. Most shp1 polygons will contain shp2 polygons, and some shp2 polygons will fall within multiple shp1 polygons.

If I start with

library(sp) large_list <- over(shp1,shp2, returnList = TRUE)

that gets me a Large List of shared geographies. But how do I take that list and use it to do a spatial join, and create a new shapefile?

I'm relatively new to R (especially for GIS) and any help would be appreciated.

Aaron M
  • 125
  • 1
  • 8
  • 1
    Is [this](http://stackoverflow.com/a/23138742/980833) what you are looking for? (Note that even though it's distributed as part of the **raster** package, `union()` does take two `SpatialPolygons*` objects as arguments.) – Josh O'Brien Apr 22 '15 at 23:43
  • 1
    @JoshO'Brien I get the following error when I do that: Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, "rgeos_difference") : TopologyException: Input geom 0 is invalid: Self-intersection at or near point -120.76777118577075 37.538496877470358 at -120.76777118577075 37.538496877470358 – Aaron M Apr 22 '15 at 23:49
  • Can you zero in on the problem areas enough to produce a minimal reproducible example that you can share with us? Otherwise, unless someone has experienced that exact same error message and remembers what caused it, we won't be able to help. (I do have a hunch that it might be the same issue that caused me to introduce an `eps` argument in my other (non-accepted) answer to that same linked question). – Josh O'Brien Apr 22 '15 at 23:57
  • Definitely check out the [`rgeos` package](http://cran.r-project.org/web/packages/rgeos/), and use `gIntersects` for a basic test, or `gWithinDistance` if you need a distance threshold between the two. I'm not sure why others are suggesting `union()`. – Mike T Apr 22 '15 at 23:59
  • They are suggestion union because it directly addresses the question; although raster::intersect might be what Aaron is after here (don't worry, both use rgeos under the hood). The error message says that a polygon is not valid, but it is not clear whether that is in the input data or occurs after some processing. If the input is valid `gIsValid(spgeom, reason=TRUE)` then you should write the raster (and/or rgeos) package maintainer as this might be a bug. – Robert Hijmans May 03 '15 at 20:00

1 Answers1

0

You might find the answer in the function intersect, from package raster

library(raster)
newshape <- intersect(shape1, shape2)
rafa.pereira
  • 13,251
  • 6
  • 71
  • 109