Loop polygon clipping with gIntersection

2
votes

Question

I am trying to perform a number of polygon clips using the gIntersection function with R in a loop. I can obtain the correct clips and re-enter data manually (so I can turn the resulting SpatialPolygons object back into a SpatialPolygonsDataFrame object). What I can't do is get this working in a loop with for() or apply().

At the moment this isn't a problem. I have 9 English regions (with London), so it's not a huge challenge to set each clip up manually. But, I want to eventually clip LSOAs in LADs, which essentially means setting up >400 clips.

So, my question is, how do I turn my manual clips into a working loop?

Minimal Reproducible Example

To keep things simple, let's use the English regions (n = 9). For each of the 9 regions, I'm going to clip the counties. The following code loads the appropriate shapefiles and reprojects them as British National Grid:

require(rgdal)
require(rgeos)
# English counties shapefile (~ 10MB zipped)
download.file(
  "https://dl.dropboxusercontent.com/s/6o0mi28pjo1kh9k/england-counties.zip",
  "ec", method = "wget")
unzip("ec")
ec <- readOGR("england-counties", "england_ct_2011")
proj4string(ec) <- CRS("+init=epsg:27700")

# English regions (~6MB zipped)
download.file(
  "https://dl.dropboxusercontent.com/s/p69m0vk2fh4xe3o/england-regions-2011.zip",
  "er", method = "wget")
unzip("er")
er <- readOGR("england-regions-2011", "England_gor_2011")
proj4string(er) <- CRS("+init=epsg:27700")

You should be left with two objects, er (English regions) and ec (English counties). Both are SpatialPolygonsDataFrame objects.

Taking the first region - East of England E12000006 - let's clip the counties and turn the result back in to a SpatialPolygonsDataFrame object:

ee <- gIntersection(ec, er[er$CODE == "E12000006", ],
                    byid = T, drop_not_poly = T)
row.names(ee) <- as.character(gsub(" 0", "", row.names(ee)))
# gIntersection adds ' 0' to each row.name?
ee <- SpatialPolygonsDataFrame(ee, ec@data[row.names(ee), ])

A plot of ee confirms this worked:

Plot of East of England region with counties

As you can see, this is a nice workflow for just a few shapes, but I really want to loop through all regions and, ultimately, many more polygons.

What I've Tried

I'm not very good with apply() loops, so what I've tried so far is a for() loop (which I know is relatively slow, but still quicker than typing everything out!):

regions <- as.character(er$CODE)  # length = 9 as expected
for(i in 1:length(regions)){
  as.name(paste0(regions[i], "c")) <- 
    gIntersection(ec, er[er$CODE == regions[1], ], byid = T, drop_not_poly = T)
}

Rather than the expected behaviour I get the following error:

Error in as.name(paste0(regions[1], "c")) <- gIntersection(ec, er[er$CODE ==  : 
could not find function "as.name<-"

I also tried wrapping the object name in an eval() but get the following error:

Error in eval(as.name(paste0(regions[1], "c"))) <- gIntersection(ec, er[er$CODE ==  :
could not find function "eval<-"

What am I missing?

In addition to the gIntersection, I would like to re-create a SpatialPolygonsDataFrame object if possible. I've tried the following code, having done one gIntersection manually, but again it doesn't work:

for(i in 1:length(regions)){
   row.names(as.name(paste0(regions[i], "c"))) <- as.character(gsub(" 0", "",
     row.names(as.name(paste0(regions[i], "c")))))
}

I get the following error:

Error in `rownames<-`(x, value) : 
  attempt to set 'rownames' on an object with no dimensions

I'm also not sure how to increment the " 0", as this increases by one for each new region (" 1", " 2", etc.)

Again, setting the first example up manually I also can't perform the final SpatialPolygonsDataFrame step:

for(i in 1:length(regions)){    
    as.name(regions[i]) <- SpatialPolygonsDataFrame(regions[i],
      ec@data[row.names(regions[i], )])
}

For this I get the following error:

Error in stopifnot(length(Sr@polygons) == nrow(data)) :
trying to get slot "polygons" from an object of a basic class ("character") with no
slots

Where I've looked

The following SO examples are related by do not seem to help, or at least I can't see how I would make them apply to this example:

Thanks for taking the time to read this.

1
rloopsgispolygonclip

1 Answers

1
votes

Does this help?

ee <- lapply(regions, function(x)  
  gIntersection(ec, er[er$CODE == x, ], byid = TRUE, drop_not_poly = TRUE))

This gives you a list of SpatialPolygonsDataFrames, one for each region. Which you can access in the usual way, e.g.

ee[[1]]
plot(ee[[1]])  # to plot the first region with counties

Edit

Your orignial code should work with a sligh modification (see blow). 

res <- list()
for (i in 1:length(regions)) {
  ee <- gIntersection(ec, er[er$CODE == regions[i], ],
                      byid = TRUE, drop_not_poly = TRUE)
  row.names(ee) <- as.character(gsub(paste0(" ", i-1), "", row.names(ee)))
  ee <- SpatialPolygonsDataFrame(ee, ec@data[row.names(ee), ])
  res[[i]] <- ee
}

If that solves the problem, then the problem was, that row names of ee always incremented by one and you did not account for this.