How to clip WorldMap with polygon in R?

Question

I have imported a world map dataset from www.GADM.org using the R package raster. I would like to clip it to a polygon I create to reduce the size of the map. I can retrieve the data and I can create the polygon no problem, but when I use the 'gIntersection' command I get an obscure error message.

Any suggestions on how to clip my World Map dataset?

library(raster)
library(rgeos)

## Download Map of the World ##
WorldMap <- getData('countries')

## Create the clipping polygon
clip.extent <- as(extent(-20, 40, 30, 72), "SpatialPolygons")
proj4string(clip.extent) <- CRS(proj4string(WorldMap))

## Clip the map
EuropeMap <- gIntersection(WorldMap, clip.extent, byid = TRUE)

Error Message:

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : Geometry collections may not contain other geometry collections In addition: Warning message: In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : spgeom1 and spgeom2 have different proj4 strings

Answer 1

You don't need to use PBS (I have learnt this the hard way, as the r-sig-geo link posted by @flowla was a question originally posted by me!). This code shows how to do it all in rgeos, with thanks to various different postings from Roger Bivand. This would be the more canonical way of subsetting, without resorting to coercion to PolySet objects.

The reason for the error message is that you can't gIntersection on a collection of SpatialPolygons, you need to do them individually. Find out which ones you want using gIntersects. I then subset each country polygon using gIntersection. I had a problem passing a list of SpatialPolygons objects back to SpatialPolygons, which turns the croppped shapefiles into SpatialPolygons, which was because not all cropped objects were of class SpatialPolygons. Once we excluded these everything works fine.

# This very quick method keeps whole countries
gI <- gIntersects(WorldMap, clip.extent, byid = TRUE )
Europe <- WorldMap[which(gI), ]
plot(Europe)


#If you want to crop the country boundaries, it's slightly more involved:
# This crops countries to your bounding box
gI <- gIntersects(WorldMap, clip.extent, byid = TRUE)
out <- lapply(which(gI), function(x){ 
        gIntersection(WorldMap[x,], clip.extent)
   })

# But let's look at what is returned
table(sapply(out, class))
#   SpatialCollections    SpatialPolygons 
#                    2                 63 


# We want to keep only objects of class SpatialPolygons                 
keep <- sapply(out, class)
out <- out[keep == "SpatialPolygons"]


# Coerce list back to SpatialPolygons object
Europe <- SpatialPolygons(lapply(1:length(out), function(i) {
          Pol <- slot(out[[i]], "polygons")[[1]]
          slot(Pol, "ID") <- as.character(i)
          Pol
   }))

plot(Europe)

If you can, I recommend you look at naturalearthdata. They have high quality shapefiles that are kept up-to-date and are constantly checked for errors (since they are open source if you find an error do email them). Country boundaries are under the Cultural buttons. You will see they are also a bit more lightweight and you can choose a resolution appropriate for your needs.

Answer 2

How about a little intermediate step? I adopted the following code mainly from R-sig-Geo and I think it should do the trick. You'll need both 'maptools' and 'PBSmapping' packages, so make sure you've got them installed. Here's my code:

# Required packages
library(raster)
library(maptools)
library(PBSmapping)

# Download world map
WorldMap <- getData('countries')
# Convert SpatialPolygons to class 'PolySet'
WorldMap.ps <- SpatialPolygons2PolySet(WorldMap)
# Clip 'PolySet' by given extent
WorldMap.ps.clipped <- clipPolys(WorldMap.ps, xlim = c(-20, 40), ylim = c(30, 72))
# Convert clipped 'PolySet' back to SpatialPolygons
EuropeMap <- PolySet2SpatialPolygons(WorldMap.ps.clipped, close_polys=TRUE)

I just tested it and it worked without any problems. It took some computation time to transform SpatialPolygons to PolySet, though.

Cheers, Florian