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. 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
World Map with Inset Map
3 messages · Hans-Jörg Bibiko, 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
On 14 Jun 2012, at 09:09, Roger Bivand wrote:
On Wed, 13 Jun 2012, Hans-J?rg Bibiko wrote:
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. 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)
Dear Roger, thank you so much! It works like a charm :)
but do remember to reset par() to defaults with:
That goes without saying. It was only a minimal example :) Regards, --Hans