21

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.

Performing reverse clip in ArcMap? shows how to do it in ArcMap.

How do I 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

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
RobinLovelace
  • 4,254
  • 5
  • 32
  • 47

3 Answers3

17

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
Spacedman
  • 63,755
  • 5
  • 81
  • 115
  • 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 06 '14 at 16:19
  • Deep down they just call GEOS which is a C-code library of geometry functions http://trac.osgeo.org/geos/ – Spacedman Aug 06 '14 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? http://svn.osgeo.org/geos/branches/3.4/ChangeLog – RobinLovelace Aug 06 '14 at 16:29
  • 1
    For sure it is developed. See the timeline http://trac.osgeo.org/geos/timeline or mailing list archives http://lists.osgeo.org/pipermail/geos-devel/ – user30184 Aug 06 '14 at 16:40
12

Answer for Simple Features:

sf package draws on Geometry Engine Open Source, and as such can access the list of commands such as st_within etc.

One such command, st_difference, will do the job:

require(sf)

# make a square simple feature
s <- rbind(c(1,1), c(1,5), c(5,5), c(5,1), c(1,1))
s.sf <-st_sfc(st_polygon(list(s)))
s.pol = st_sf(ID = "sq", s.sf)

# make a smaller triangle simple feature
s2 <- rbind(c(2,2), c(3,4), c(4,2), c(2,2))
s2.sf <-st_sfc(st_polygon(list(s2)))
s2.pol = st_sf(ID = "tr", s2.sf)

# find the 'difference', i.e. reverse of st_intersection
t <- st_difference(s.pol,s2.pol)

plot(t)

# have a look at the new geometry, a well known text format with exterior followed by hole
st_geometry(t)[[1]]
POLYGON((1 1, 1 5, 5 5, 5 1, 1 1), (2 2, 4 2, 3 4, 2 2))

also see towards the bottom of this article

can also be done by coercing Sp to sf with st_as_sf. Heed the warnings as attributes can be tricky to manage!

Sam
  • 1,591
  • 2
  • 15
  • 37
9

Bit late to the party, but there is a simple way to do this with mask using the 'inverse' argument;

ukhole <- mask(uk, lnd, inverse = TRUE)