rmapsgisr-maptoolsrmaps

Using a projection with xlim/ylim in R with the "maps" package results in a bigger map


I'm using R with the maps/mapproject/maptools packages to plot some maps and noticed a behavior which seems counter-intuitive to me and actually limits what I can do.

Drawing a map of Europe (with the limits taken from the ETRS89 / ETRS-LCC, so without Iceland and also clipped in the East) without specifying any projection:

library(maps)
map("world",xlim=c(-10.67,34.5),ylim=c(31.55,71.05), interior = T)

Europe map without projection specified

The result is as expected, the limits are being used and the resulting map follows them.

The projection used by default by maps is, as per the help:

The default is to use a rectangular projection with the aspect ratio
chosen so that longitude and latitude scales are equivalent at the 
center of the picture.

This is not a good projection for my needs, I will use an LCC projection with the parallels as indicated in the above spatialreference.org link:

library(maps)
map("world",xlim=c(-10.67,34.5),ylim=c(31.55,71.05), interior = T, projection="lambert", parameters = c(45,65))
box()

Europe map with projection specified

The result is unexpected in that it includes a much bigger area (going very far north and including Russia), essentially making the map unusable.

What's stranger is that when using a grid the original limits are clearly considered:

library(maps)
library(mapproj)
map("world",xlim=c(-10.67,34.5),ylim=c(31.55,71.05), interior = T, projection="lambert", parameters = c(45,65))
map.grid(cex=0.1 , col="grey30")
box()

Europe map with LCC projection and grid

What I would like to have (and what I assumed would be the result of the above code) was a rectangular crop that included the limits I specified (adjusted due to the projection used so having more area than the rectangular one above would be expected). Additionally there is a white space around the map and the border which is present whenever one uses a projection with map().

The question is: is there a way to have this result when using map/mapproj/maptools? I tried to artificially change xlim/ylim without good results since it seems to work in big intervals (i.e. changing them doesn't produce an effect until suddenly half of Europe disappears with the next decrement).


Solution

  • The bounds given on spatialreference.org are given as "bottom left, top right coordinate", which coincidentally comes very close to what you wrote. Correctly these coordinates would be:

    xlim=c(-10.6700,31.5500), ylim=c(34.5000,71.0500)
    

    But that does not seem to be the problem here.

    Playing around with the second xlim parameter it appears that map() tries to include neighbouring lines, as long as some of those line's coordinates are within the limits. The country of Russia fits in there (which causes a huge jump around 29-31 on the 2nd xlim parameter), along with some African and eastern European countries. Turning on and off boundary somewhat confirms this, but this hides a lot of country's borders.

    A workaround I found was to explicitly exclude neighbouring countries from being drawn at first. Then draw the map a second time using add=T. Use col="white" on the first drawing to not draw in the same place twice, it thickens the lines.

    library(maps)
    library(mapproj)
    map("world",regions="(?!Russia|Morocco|Algeria|Tunisia|Turkey|Ukraine)",col="white",xlim=c(-10.6600,31.5500), ylim=c(34.5000,71.0500), interior=T, projection="lambert", parameters=c(10.44,52.775))
    map("world",xlim=c(-10.6600,31.5500), ylim=c(34.5000,71.0500), interior=T, projection="lambert", parameters=c(10.44,52.775),add=T)
    map.grid(cex=0.1, col="grey30")
    box()
    

    There seem to be some lines missing still due to xlim and ylim, removing that will again include Iceland and other countries and islands. Or maybe there's a better approach altogether.

    resulting plot