Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

A reverse clip saves only the part of your spatial object that is outside the bounds of another object, as opposed to a regular clip which saves the parts that are inside the other object. Here's how to do it in ArcMap.

The question is: how to do this in R?

Reproducible example (on Linux machines):

system("wget 'https://github.com/Robinlovelace/Creating-maps-in-R/archive/master.zip' -P /tmp/")
unzip("/tmp/master.zip", exdir = "/tmp/master")
uk <- readOGR("/tmp/master/Creating-maps-in-R-master/data/", "ukbord")
lnd <- readOGR("/tmp/master/Creating-maps-in-R-master/data/", "LondonBoroughs")
plot(uk)
plot(lnd, add = T, col = "black")

What I want here to do is to save all of the UK except for London. Visually, I want the black shape in the resulting image to be a hole.

enter image description here

share|improve this question
    
+1 for posting a reproducible example! –  Curlew Aug 6 at 16:07

1 Answer 1

up vote 6 down vote accepted

Seems to be a simple application of gDifference from the rgeos package:

> require(rgeos)
> ukhole = gDifference(uk, lnd)
Warning message:
In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_difference") :
  spgeom1 and spgeom2 have different proj4 strings
> plot(ukhole)

The projection warning is because the LondonBoroughs shapefile doesn't have a .prj file.

Just to make sure that's a hole and not an outline or another solid polygon:

> gArea(lnd) + gArea(ukhole) - gArea(uk)
[1] 0
share|improve this answer
    
Sooo simple, thanks for the fast response. Would be interested in looking at the source code of these functions to see what's going on under the hood. –  RobinLovelace Aug 6 at 16:19
    
Deep down they just call GEOS which is a C-code library of geometry functions trac.osgeo.org/geos –  Spacedman Aug 6 at 16:22
    
Interesting - and helps explain why it's reasonably fast I guess. Based on this page, seems it's not being actively developed can anyone confirm/refute this? svn.osgeo.org/geos/branches/3.4/ChangeLog –  RobinLovelace Aug 6 at 16:29
1  
For sure it is developed. See the timeline trac.osgeo.org/geos/timeline or mailing list archives lists.osgeo.org/pipermail/geos-devel –  user30184 Aug 6 at 16:40

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.