More maps in R

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.

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

The mar to 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 height and 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 xlim and 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.