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

Question

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)

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()

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()

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).

Answer 1

The basic behaviour of "map()" is to extract all the lines within xlim and ylim, but also lines that are partially in the domain are extracted from the database. When there is no projection, xlim and ylim are also used as limitations for the plot, but when the data is projected that would obviously not work. (You can not simply project the xlim and ylim, as this would not be a rectangular domain). The maps code is not sophisticated enough to find a nice solution and just sets the plot limits to the range of the data (which includes any unwanted lines). As (recent) maintainer of 'maps', this is actually something I may try to fix. Unfortunately, mapproj doesn't include inverse projections, which are required to do it correctly like you would expect. Limiting the map at the bottom (in stead of extending at the top) is much easier.

As a workaround, you could also try the following (using the output value of map):

mymap <- map("world",xlim=c(-10.6600,31.5500), ylim=c(34.5000,71.0500),  projection="lambert", parameters=c(10.44,52.775), plot=FALSE)

Now you can look at the range of mymap$x and mymap$y:

range(mymap$x, na.rm=TRUE)

or just plot(mymap) to see the range and decide which interval you need.

Then, finally, you can plot a map with exactly the limits you want, for instance

plot(mymap, type="l", asp=1,xaxt="n",yaxt="n",xlab="",ylab="", xlim=c(-0.2,0.15), ylim=c(-0.75,-0.35))

Not ideal, but hopefully good enough. The "asp=1" is to keep the aspect ratio on the map correct, even if you resize the plot window.

Answer 2

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.