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

[Hans-Jörg Bibiko]

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