13

I have these two polygons:

library(sp); library(rgeos); library(maptools)

coords1 <- matrix(c(-1.841960, -1.823464, -1.838623, -1.841960, 55.663696,
                    55.659178, 55.650841, 55.663696), ncol=2)
coords2 <- matrix(c(-1.822606, -1.816790, -1.832712, -1.822606, 55.657887,
                    55.646806, 55.650679, 55.657887), ncol=2)
p1 <- Polygon(coords1)
p2 <- Polygon(coords2)
p1 <- Polygons(list(p1), ID = "p1")
p2 <- Polygons(list(p2), ID = "p2")
myPolys <- SpatialPolygons(list(p1, p2))
spdf1 = SpatialPolygonsDataFrame(myPolys, data.frame(variable1 = c(232,
                                                                   242), variable2 = c(235, 464), row.names = c("p1", "p2")))
proj4string(spdf1) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
plot(spdf1, col="red")

coords1a <- matrix(c(-1.830219, -1.833753, -1.821154, -1.830219, 55.647353,
                     55.656629, 55.652122, 55.647353), ncol=2)
p1a <- Polygon(coords1a)
p1a <- Polygons(list(p1a), ID = "p1a")
myPolys1 <- SpatialPolygons(list(p1a))
spdf2 = SpatialPolygonsDataFrame(myPolys1, data.frame(variable1 = c(2),
                                                      variable2 = c(3), row.names = c("p1a")))
proj4string(spdf2) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
plot(spdf2, col="yellow", add=T)

enter image description here

I want to cut out parts of spdf1 that are intersected by spdf2. However, I want spdf1 to remain as a SpatialPolygonsDataFrame and to retain any information contained within spdf1@data.

I¹ve tried gDifference as follows, which cuts out parts of spdf1 that are intersected by spdf2, but then converts spdf1 to SpatialPolygons (i.e. discarding information contained in spdf1@data).

gDifference(spdf1, spdf2, byid=T)

How can I cut into spdf1 with spdf2 but retain data contained in spdf1@data? I have checked and tried this similar question without how to overlay a polygon over SpatialPointsDataFrame and preserving the SPDF data?

9

The simplest approach would be

  library(raster)
  x <- spdf1 - spdf2

  # or, more formally
  y <- erase(spdf1,  spdf2)

See ?'raster-package' (section XIV) for more functions that deal with polygon overlay. These functions use the base-functions of rgeos under the hood, in 'user-level' (as opposed to 'developer-level') functions.

  • What do you mean by "in 'user-level' (as opposed to 'developer-level') functions"? – luciano Nov 9 '15 at 21:16
  • rgeos provides geometric operations but does not deal with the attributes of the data. Thus using these functions requires a lot of work to keep everything together. The raster functions simplify this and the behave like similar functions in GIS software, – Robert Hijmans Nov 9 '15 at 21:27
  • Yeah, but this can lead to Error in SpatialPolygonsDataFrame(part2, x@data[match(row.names(part2), : row.names of data and Polygons IDs do not match – jebyrnes Jan 12 '18 at 16:51
  • That would be a bug. Can you provide an example of that? – Robert Hijmans Jan 13 '18 at 18:31
4

A workaround would be to re-add the attributes after doing the clip, while converting from SpatialPolygons to SpatialPolygonsDataFrame.

sp3 <- gDifference(spdf1, spdf2, byid = TRUE)
row.names(sp3) <- row.names(spdf1)

spdf3 <- SpatialPolygonsDataFrame(sp3, data = spdf1@data)

spdf3@data

   variable1 variable2
p1       232       235
p2       242       464

plot(spdf3, col="red")

enter image description here

  • This looks like a problem I have, only the output clip in my particular instance produces rownames from spdf1 that do not exist (as a simple gsub to get rid of the 2nd digit in row.names should take care of the rownames match, no?) – jebyrnes Jan 12 '18 at 17:11
  • Error in SpatialPolygonsDataFrame(clip, data = as.data.frame(all_spdfs_together@data)) : Object length mismatch: clip has 1718 Polygons objects, but as.data.frame(all_spdfs_together@data) has 86 rows – jebyrnes Jan 12 '18 at 17:18
  • Sure - I have a bunch of polygons of things in the ocean. Some are incorrectly placed on land, or overlap with land. I want to clip that out so I only have areas that are in the ocean. I have a coastline shapefile to compare against. Here's the code - gist.github.com/jebyrnes/c2e8d2b6c82849dad3a813d952ab8bb0 – jebyrnes Jan 12 '18 at 20:12
  • 1
    nevermind - the raster::erase solution works (it hadn't before for some odd reason) – jebyrnes Jan 12 '18 at 20:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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