[R-sig-Geo] Subset polygon(s) from a country from wrld_simpl

[Roger Bivand]

On Fri, 6 Apr 2012, Lars Dalby wrote:

> Dear list
>
> I am trying to make a subset of wrld_simpl to use as mask for a raster
> layer.
> I have a list of countries and I can do the subset quite easily. However I
> would like to have one (or several) of the polygons from Russia in the list
> also - Kaliningrad on the coast of the Baltic. This is where I run into
> trouble.
> How do I cut these out and subsequently subset wrld_simpl accordingly?
>
> require('sp');require('maptools')
>
> cntr <- c('Denmark', 'Germany', 'Norway', 'Ireland', 'United Kingdom',
> 'France', 'Italy', 'Sweden', 'Finland', 'Spain', 'Portugal', 'Latvia',
> 'Estonia', 'Slovenia', 'Belgium', 'Netherlands', 'Austria', 'Poland',
> 'Switzerland', 'Morocco', 'Algeria', 'Tunisia', 'Slovakia', 'Lithuania',
> 'Croatia', 'Czech Republic')
>
> data(wrld_simpl)
> NW <- wrld_simpl[which(wrld_simpl at data$NAME %in% cntr),]
>
> #I tried to use the overlay function, but that didn't quite do the trick:
>
> rus <- wrld_simpl[which(wrld_simpl at data$NAME == 'Russia'),]

length(slot(rus, "polygons"))

shows that each country is a single Polygons object (by the way, 
wrld_simpl at data$NAME should be wrld_simpl$NAME). If you need to get below 
that, you have to dismantle the Polygons object:

>
> x <- 21.335449; y <- 54.788017
> P <- data.frame('x' = x, 'y' = y)
> coordinates(P) <- ~x+y

proj4string(P) <- CRS(proj4string(rus))
over(rus, P)
lrus <- slot(slot(rus, "polygons")[[1]], "Polygons")
for (i in seq(along=lrus)) {
   lrus[[i]] <- Polygons(list(lrus[[i]]), as.character(i))
}
SPrus <- SpatialPolygons(lrus, proj4string=CRS(proj4string(rus)))
pip <- which(over(SPrus, P) == 1)
kaliningrad <- SPrus[pip,]
row.names(kaliningrad) <- "Kaliningrad"
NW1 <- spRbind(NW, kaliningrad)
plot(NW1)

Should Kaliningrad become an ISO country, it would be promoted at some 
stage. As of now, I think that South Sudan is not distinguished in 
maptools, but I'm not sure about rworldmap; it isn't in cshapes which 
covers changes 1946-2008. The data is intended to cover whole countries at 
a fairly coarse vector scale; maybe GADM has the administrative boundaries 
to let you zoom in.

Roger

> ov <- overlay(P, rus)
> ov # 1
>
> # Now if I move the point further east and repeat:
>
> x <- 90.335449; y <- 54.788017
> P <- data.frame('x' = x, 'y' = y)
> coordinates(P) <- ~x+y
> ov <- overlay(P, rus)
> ov # still 1- so I guess this refers to the single object in the polygons
> slot of rus, which is not quite what I'm after.
>
> #sessionInfo()
> R version 2.14.1 (2011-12-22)
> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
>
> locale:
> [1] C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] dismo_0.7-17    maptools_0.8-14 lattice_0.20-6  foreign_0.8-49
> [5] raster_1.9-70   sp_0.9-97       ncdf_1.6.6
>
> loaded via a namespace (and not attached):
> [1] grid_2.14.1  tools_2.14.1
>
>
> Any suggestions on how to solve this would be highly appreciated.
>
> Best regards and Happy Easter!
>
> Lars
>
> Lars Dalby
> Msc. PhD student
> Ecoinformatics & Biodiversity
> Department of Bioscience, Aarhus University
> Denmark
>
> --
> View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Subset-polygon-s-from-a-country-from-wrld-simpl-tp7442490p7442490.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> 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