[R-sig-Geo] Crop a SpatialPolygonsDataFrame by a bounding box/anotehr polygon

[O'Hanlon, Simon J]

Dear list,
For those interested I managed to hack out a solution thanks to some old posts from this thread and others. Basically PBSmapping::clipPolys() does what I need, but then you have to fiddle around converting SpatialPolygons to PolySets and back again. If anyone knows a way I can solve my original problem of how to clip SpatialPolygons by another polygon without coercing to PolySets I'd be most grateful. Below is my code which now works (however it is not 'bulletproof' - you have to subset the original shapefile to only the country polygons which will appear in the final plot otherwise you will get mismatch between the resulting polygon IDs and the rownames of the original dataset when you convert back from a PolySet to a SpatialPolygonsDataFrame object. I'm sure this can be solved rather easily too, but this rough solution works for me for now).

New code below:

Thanks,

Simon
---

#	Load relevant packages
library(ggplot2)
library(maptools)
library(rgeos)
library(plyr)
library(PBSmapping)


#	Download and unzip the Natural Earth Country Polygon data
oldwd <- getwd()
tmp <- tempdir()
setwd(tmp)
url <- "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/50m-admin-0-countries.zip"
dest <- paste(tmp,"\\tmp.zip",sep="")
download.file(url,dest)	#File is 1.3Mb
unzip(dest)


#	Read in the data to a SpatialPolygonsDataFrame and subset to West African countries I wish to appear in plot
wld <- readShapePoly("ne_50m_admin_0_countries")
wa <- wld[c(21,22,43,81,82,83,84,144,151,158,160,190,195,213),]


#	Convert to PolySet, clip polygons and convert back to SpatialPolygons
wa.ps <- SpatialPolygons2PolySet(wa)
wa.c <- clipPolys(wa.ps,xlim=c(-18,4),ylim=c(4,18))
wa.sp <- PolySet2SpatialPolygons(wa.c, close_polys=TRUE)


#	To attach the original attribute data to the new polygons the rownames of the data must match the polygon IDs
rn <- sapply(slot(wa.sp, "polygons"), function(x) slot(x, "ID"))
rownames(wa at data) <- rn
wa.new <- SpatialPolygonsDataFrame(wa.sp,data=wa at data)


#	Work out the centroids of the new polygons 
wa.df <-as.data.frame(wa.new)
wld.lab <- data.frame(country = (wa.df$NAME), coordinates(wa.new))


#	Plot map and country labels
plot(wa.new)
text(wld.lab[,2:3], labels=wld.lab$country, col="red")


#	Cleanup
unlink(tmp,recursive=T)
setwd(oldwd)

-----Original Message-----
From: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of O'Hanlon, Simon J
Sent: 21 June 2012 11:29
To: r-sig-geo at r-project.org
Subject: [R-sig-Geo] Crop a SpatialPolygonsDataFrame by a bounding box/anotehr polygon

Dear all,
I would like to crop a particular contiguous area of country polygons to a rectangular bounding box, and return new polygons, from which I can work out the centroids. I would like to do this so that I can use the centroids of the cropped country polygons to work out nice label placements for a map. I am having trouble cropping a SpatialPolygonsDataFrame object by another polygon using gIntersection (I suspect this is not for this purpose?).

Can anyone help edit the following to make it work? The gIntersection command fails. I think I need to find an alternative, but I am not sure what to do yet. If anyone can suggest a command I'd be most grateful.

Many thanks in advance,

Simon
---
#	Load relevant pacakges
library(ggplot2)
library(maptools)
library(rgeos)
library(plyr)

#	Download and unzip the Natural Earth Country Polygon data
oldwd <- getwd()
tmp <- tempdir()
setwd(tmp)
url <- "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/50m-admin-0-countries.zip"
dest <- paste(tmp,"\\tmp.zip",sep="")
download.file(url,dest)	#File is 1.3Mb
unzip(dest)

#	Read in the data to a SpatialPolygonsDataFrame
wld <- readShapePoly("ne_50m_admin_0_countries")

#	Create bounding box polygon which we would like to use to crop polygons
bb <- readWKT("POLYGON((-18 4, 4 4, 4 18, -18 18, -18 4))")
proj4string(bb) <- proj4string(wld)

#	Use gIntersection to crop the polygons to the borders of the bounding box - THIS FAILS
wa <- gIntersection(wld,bb)

#	Work out the centroids of the new polygons (but we can't do that until we have the new polygons!)
wa.df <-as.data.frame(wa)
wld.lab <- data.frame(country = (wa.df$NAME), coordinates(wa))

#	Cleanup
unlink(tmp,recursive=T)
setwd(oldwd)

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo