Lots of students want maps showing their field sites for their thesis. Several have come to me with code they found on the internet but couldn’t get to work. The problem is that the plotting package
ggplot2 has evolved since that code was written and some things no longer work (for example
opts instead of
ggplot2 now supports map projections and Norway is more complicated than Canada in that it has remote territories.
This post updates the code from code to make a map with an inset showing the location against a regional map.
I’m going to plot the seedClim locations.
Load the required packages and some data.
library(ggplot2) library(maps) library(mapdata) library(grid) dat <- read.table(header = TRUE, text = " siteID latitude longitude Temperature Precipitation Alrust 60.8203 8.70466 2 1 Arhelleren 60.6652 6.33738 3 3 Fauske 61.0355 9.07876 3 1 Gudmedalen 60.8328 7.17561 1 3 Hogsete 60.876 7.17666 2 2 Lavisdalen 60.8231 7.27596 1 2 Ovstedal 60.6901 5.96487 3 4 Rambera 61.0866 6.63028 2 3 Skjellingahaugen 60.9335 6.41504 1 4 Ulvhaugen 61.0243 8.12343 1 1 Veskre 60.5445 6.51468 2 4 Vikesland 60.8803 7.16982 3 2 ") dat$Temperature <- factor(dat$Temperature) dat$Precipitation <- factor(dat$Precipitation) precipLab <- c("Very dry", "Dry", "Wet", "Very wet") tempLab <- c("Alpine", "Intermediate", "Lowland")
Find the latitude and longitude ranges required for the high resolution map. Here I’ve added a buffer round the locations.
xlim <- range(dat$longitude) + c(-1, 0.5) ylim <- range(dat$latitude) + c(-0.5, 1)
Get the map data
# Low resolution map of Norway for the inset map norwaymap <- map_data("world", "Norway") # High resolution map of Norway for the main norwaymapHires <- map_data("worldHires", "Norway")
Unfortunately, even the high resolution map is not high-resolution enough – lots of important islands off the Hordaland coast are missing, for example Askøy near Bergen. I need to find a higher resolution map.
Make the maps
The inset map showing all of Norway
a <- ggplot() + geom_map(data = norwaymap, aes(x = long, y = lat, map_id = region), map = norwaymap, colour = NA, fill = "grey60") + geom_rect(data = data.frame(), aes(xmin = xlim, xmax = xlim, ymin = ylim, ymax = ylim), colour = "red", fill = NA) + coord_map(xlim = c(3, 33), ylim = c(57, 72)) + labs(x = NULL, y = NULL) #print(a)# print this map until you are happy with it.
geom_map draws the map.
geom_rect shows the limits of the high-resolution map.
coord_map projects the map. Any of the projections in the
mapproj package can be used – there are probably better ones than the default Mercator projection. I had to set the x and y limits to exclude Svalbard, Jan Mayen and Bouvetøya (a subantarctic island).
Zoomed in map, with the site locations
b <- ggplot(dat, aes(x = longitude, y = latitude, colour = Precipitation, shape = Temperature, fill = Precipitation)) + geom_map(aes(x = long, y = lat, map_id = region), data = norwaymapHires, map = norwaymapHires, colour = NA, fill = "grey60", inherit.aes = FALSE) + geom_point(size = 3) + coord_map(xlim = xlim, ylim = ylim) + labs(x = NULL, y = NULL) + scale_shape_manual(breaks = 1:3, labels = tempLab, values = c(24, 22, 25)) + guides(shape = guide_legend(override.aes = list(fill = "black")))+ scale_colour_brewer(breaks = 1:4, labels = precipLab, palette = "RdBu") + scale_fill_brewer(breaks = 1:4, labels = precipLab, palette = "RdBu") #print(b)
I’m using various
scale_?_??? commands to set the shape and colour of the points. I need to use the
guides function as otherwise the temperature legend has empty rather than filled shapes.
maptheme <- theme( axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), panel.border = element_rect(fill = NA, colour = "black"), panel.background = element_blank() )
This theme can be added to the plots. It sets a dark rectangle round the map and removes the grey background and axis annotation. An alternative would be
ggthemes::theme_map() but it moves the legend.
Printing the maps
Print the two maps together, one an inset of the other using the
grid package for flexible alignment of
grid.newpage() vp_b <- viewport(width = 1, height = 1, x = 0.5, y = 0.5) # the larger map vp_a <- viewport(width = 0.4, height = 0.4, x = 0.6, y = 0.8) # the inset in upper left print(b + maptheme, vp = vp_b) print(a + maptheme, vp = vp_a)