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

[Roger Bivand]

On Thu, 21 Jun 2012, O'Hanlon, Simon J wrote:

> Dear list,

> For those interested I managed to hack out a solution thanks to some old

No, this was answered in a thread here starting at:

https://stat.ethz.ch/pipermail/r-sig-geo/2012-June/015337.html

just a week ago. If you cite threads, always give a link. Note that 
PBSmapping uses gpclib code, although the PBS authors claim to have 
permission to release it under GPL. If you need to stay on the safe side, 
rgeos solves the problem (also for the planar case) without encumbrances.

Roger


> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no