It was perhaps inevitable: as soon as I had shown the delegates at an R workshop in Tallinn that maps of Estonia could now be drawn with the
maps package, they wanted to know how to project maps, add site locations, and various other tweaks.
maps package is very easy to use, but a few tricks are needed to get the best out of it and to prevent some surprising behaviour.
I like maps with the land shaded (it can be bit more more difficult to shade the sea, at least with projected maps). Here is a simple map of Scandinavia and the Baltic region.
library(maps) x11(width = 3, height = 6.5) map(xlim = c(4, 35), ylim = c(55, 81), mar = rep(0.1, 4), col = "grey80", fill = TRUE) box()
map() argument sets the margins of the map. I like to keep these very small and roughly the same on each edge. This means I have to tinker with the
width arguments of
x11() to get everything right. I like this map, but there is a glitch: Denmark appears to float above a void which should be occupied by Germany. I think the problem is that
map() prefers to plot whole polygons when using the argument
fill = TRUE.
We can fix this and make other improvements to this map.
loc<-data.frame(longitude=5.3300, latitude=60.3894) x11(width = 3, height = 6.5) #set up map coordinates with bank map map(xlim = c(4, 35), ylim = c(55, 81), mar = rep(0.1, 4), type="n") #plot map map(col = "grey80", fill = TRUE, add = TRUE) #make lake white map(regions = "Lake*", add = TRUE, col= " white", fill = TRUE) box() #add site locations points(loc$long, loc$lat, col=2, pch=16) #add scale bar. Scale correct at latitude of bar map.scale(y = 72, ratio = FALSE, col = 2, relwidth = 0.3) #add north pointer arrows(x0 = 5.5, x1 = 5.5, y0 = 72, y1 = 78) text(x = 5.5, y = 78, "N", pos = 3)
The trick to many of the best things with map is to first call
map() with argument
type = "n" which does not plot anything but sets up a coordinate system that subsequent calls to
map() with argument
add = TRUE can use. In the first call, I ask for the region within the
ylim. In the second, I ask for the whole world, but only the region that fits into the plot is shown. The lake on the Estonian-Russian border looked confusing (micro-nations?), so I’ve filled it white to make it clear. Note, only lakes on national borders are included in the database. If you want to show other lakes (for example Lake Ladoga, east of Saint Petersburg) you need to find an alternative data source.
The scale bar was a bit of a pain as it wanted to have an tick label in the middle which overlaps with the final tick label. This might not be such a problem in a map with a square aspect ratio. I fixed it here by increasing the
relwidth, but was very tempted to edit the function to make it use only the first and last tick label.
map() uses by default a rectangular projection with the aspect ratio set so latitudinal and longitudinal scales are equal in the centre of the map. That is OK for maps covering a small region (for example Estonia) but for larger regions, distortion is inevitable. In the maps above, Svalbard appears to be much larger than it should be. Distortion can be reduced (or increased!) by using one of the many projections that are available with the
mapproj package. I’m going to use a sinusoidal projection to map the same area as above.
library(mapproj) x11(width = 3, height = 6.5) #set up plot coordinates map(xlim = c(4, 35), ylim = c(55, 81), mar = rep(0.1, 4), fill = FALSE, proj = "sinusoidal", type = "n") #add map with same projection as first with proj = "" #set xlim to prevent wrapping map(xlim = c(-130,130), col = "grey80", fill = TRUE, proj = "", add = TRUE) box()
That was not good. The first problem is that with projected maps the
mar argument is ignored and needs to be set with
par() first (this is standard for most other plotting functions). The second bigger problem is that the first call to
map included Russia in the map, so the second call then plots most of the western hemisphere. We can fix this by explicitly telling the code which countries should be used to set the area the map covers.
x11(width = 4, height = 6.5) par(mar = rep(0.1,4)) #force map to cover Norway, Finland and Latvia map(region = c("Norway", "Finland", "Latvia"), proj = "sinusoidal", type = "n") #map everything that does not start with "Lake" map(xlim = c(-130, 130), region = "(?!Lake *)", col = "grey80", fill = TRUE, proj = "", add = TRUE) #add projected site locations points(mapproject(x = loc$long, y = loc$lat), col = 2, pch = 16) #add a latitude-longitude grid map.grid(c(0,40, 55, 80), nx = 5, ny = 5) box()
It takes a bit of tweaking of the list of countries and the size of the plotting window to get the map to cover the required area. I think this looks fairly good.
map.scale() does not work with projected maps, and the code for adding north arrows won’t work with many, so
map.grid() can be useful.
Plotting maps has a certain artistic element, and it generally takes a bit of trial and error to get things just so. Other improvements that can be made include inset maps (hint use
layout()), shapefiles of rivers and lakes, and background raster images. It is also possible to plot maps with the
ggplot2 package, but I haven’t explored this yet.