[R-sig-Geo] World Map with Inset Map

[Roger Bivand]

On Wed, 13 Jun 2012, Hans-Jörg Bibiko wrote:

> Hi,
>
> maybe I've overlooked something obvious but I'm trying to generate _one_ plot of a world map (actually projected) with an inset of a specific zoomed area based on SpatialDataFrames.
>
> I followed this example which makes usage of the library "maps":
> http://wiki.cbr.washington.edu/qerm/sites/qerm/images/7/78/MakingAnInsetMapShorter.r
>
> Thus I wrote this (as a minimal and naïve example):
>
> library(maptools)
> data(wrld_simpl)
>
> # plot world map
> plot(wrld_simpl, col = "khaki", bg = "azure2")
>
> # define inset map location and size
> par(plt = c(0.57, 0.87, 0.4, 0.7), new = TRUE)
> # init inset map view port
> plot.window(xlim = c(130, 180), ylim = c(40, 70))
>
> # draw background
> rect(130, 40, 180, 70, col = "azure2")
>
> # plot zoomed inset map
> plot(wrld_simpl, xlim = c(130, 180), ylim = c(40, 70),
> 	col = "khaki", bg = "azure2", add = TRUE)
>
>
> The problem is that the inset map won't be cropped, i.e. you see both 
> maps overlapped.

The plot methods make no promise to crop, in addition, when plot(..., 
add=TRUE), the xlim= and ylim= arguments are ignored, because they would 
have been set in the plot.window() call. If you crop yourself, you can get 
there:

library(rgeos)
bnds <- cbind(x=c(130, 130, 180, 180, 130), y=c(40, 70, 70, 40, 40))
SP <- SpatialPolygons(list(Polygons(list(Polygon(bnds)), "1")),
  proj4string=CRS(proj4string(wrld_simpl)))
gI <- gIntersects(wrld_simpl, SP, byid=TRUE)
out <- vector(mode="list", length=length(which(gI)))
ii <- 1
for (i in seq(along=gI)) if (gI[i]) {out[[ii]] <- 
gIntersection(wrld_simpl[i,], SP); row.names(out[[ii]]) <- 
row.names(wrld_simpl)[i]; ii <- ii+1}
out1 <- do.call("rbind", out)
plot(out1, col = "khaki", bg = "azure2", add = TRUE)

but do remember to reset par() to defaults with:

oldpar <- par(plt = c(0.57, 0.87, 0.4, 0.7), new = TRUE)

before the insert and:

par(oldpar)

afterwards!

Roger

>
> What I'm doing wrong? Or is this with SpatialDataFrames not possible?
>
> I'd be appreciated for any hint!
>
> Kind regards,
> --Hans
>
> PS = sessionInfo():
>
> R version 2.14.2 (2012-02-29)
> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] maptools_0.8-14 lattice_0.20-0  sp_0.9-96       foreign_0.8-49
>
> loaded via a namespace (and not attached):
> [1] grid_2.14.2  tools_2.14.2
>
> _______________________________________________
> 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