Creating a map with inset using tmap
Plotting maps with tmap
The tmap R package greatly simplifies creating good looking maps. Like in ggplot, but also similar to most GIS tools, it allows you to build up a map layer by layer. If tmap is new to you, good starting points are the Geocomputation with R online book and the tmap get started page. And you can of course google for examples.
One of the features I like, is the ability to include inset maps, like in Figure 1. To be able to do this, you need the viewport function of the grid package. This is not a function I use a lot, so every time I want to make a map with an inset map, I have to look up the details. One of those details is how to properly align the legend and the inset map. As a note to self, but perhaps also useful to others, below the code I used to create this figure.
Get the libraries
First step is to load the required libraries. Besides the aforementioned tmap and grid libraries, the following packages are used: the tmaptools, sf, sp and rgrass7 libraries. The rgrass7 library is required because I want to use two layers I keep in a GRASS GIS database. If you want to import for example a geotif or vector layer in a geopackage, you can use the raster package and sf packages instead.
Main map
The readRAST and readVECT functions are used to read in the raster and vector data from a GRASS GIS database. To import raster data, rgrass7 needs to be instructed to use the older sp package, while for the vector layer, the sf package can be use. This is done by the functions use_sp and use_sf respectively.
Next, create a vector layer of the study area based on the bounding box of the raster layer, using the bb_poly function from the tmaptools package. This is later used to show the location of the study area in the inset map.
To be able to properly align the different map elements, you need to know the aspect ratio of the input layers. To calculate this, extract the lower and upper x- and y-coordinates using the st_bbox function from the sf package. With those values, calculate the width and height, and subsequently the aspect ratio.
Now you have all the data, you can create the main map. I won’t go into all the details, the functions and parameters are clearly explained in the help pages. Just make sure to define the inner and outer margins (here both set to 0), the position of the legend, and the legend height and width. Important, you’ll notice that I used negative numbers for the legend width and height. The minus sign is used to tell tmap to use the exact width and length as provide. If you give a positive number, tmap will decrease the actual legend width and height automatically based on the legend content and font sizes.
mainmap <- tm_shape(rd) +
tm_raster(style="cont", title="Distance (m)") +
tm_shape(n2) +
tm_fill(col="white", alpha=0.5, lwd=0.5) +
tm_borders(col="lightgrey", lwd=0.01) +
tm_layout(legend.outside=F, legend.text.size = 0.8,
legend.title.size=0.9,
legend.frame=TRUE, legend.position=c(0.985, 0.985),
legend.just = c("right", "top"), legend.width=-0.25,
legend.height=-0.4, outer.margins = c(0,0,0,0),
inner.margins = c(0,0,0,0))Inset map
The raster layer covers a region in the central Oder river basin in Poland. The idea is for you to add an inset map of Poland, including the region covered by the raster layer. You can get the boundary data from the GADM database of Global Administrative Areas. The data can be downloaded in various formats, including the rds data format. The rds file contains a sf vector object. By combining the readRDS and url functions, you can download and import the required files in one go (you’ll need to get the download urls first of course).
### Get boundary data
link1 <- "https://biogeo.ucdavis.edu/data/gadm3.6/Rsf/gadm36_POL_0_sf.rds"
link2 <- "https://biogeo.ucdavis.edu/data/gadm3.6/Rsf/gadm36_POL_1_sf.rds"
PolCntr <- readRDS(file = url(link1, "rb"))
PolProv <- readRDS(file = url(link2, "rb"))
### Project to EPSG:3035 (the same crs as the raster data)
PolCntr <- st_transform(PolCntr, crs=3035)
PolProv <- st_transform(PolProv, crs=3035)The two layers will be used as input for the inset map. To properly align this inset map with the legend, we need to know the aspect ratio of the vector data. Again, use the st_bbox function to get the minimum and maximum x- and y-values, and use these to calculate the aspect ratio.
Now you’re set to create the inset map. Set the outer.margins to 0. If you want to set the inner margins to a certain value (as I have done here), make sure they are the same on all sites. Otherwise, you’ll need to adapt the aspect ratio.
Create the final map
For the final map, use the viewport function to determine where the inset map will be printed, and at what size. For the latter, define the width and height. The ratio of these two need to correspond to the aspect ratio of the inset map. So first set the width to match the width of the legend of the main map. Next, compute the corresponding height by multiplying the width by the aspect ration you calculated earlier.
Now you just need to ensure the inset map aligns neatly with the legend by setting the x and y and the justification of the viewport relative to the x and y location.
And, as the very last step, save the final map as a png image. To avoid white spaces on the top and bottom, or left and right, make sure the aspect ratio of the output image is the same as that of the map.
All the code
Quite a bit of code for such a simple map. But on the other hand, once you have created the map, it is very easy to update or adapt it if for example new data comes available.
# libraries
library(tmap)
library(tmaptools)
library(raster)
library(rgrass7)
library(grid)
library(sf)
tmap_mode("plot")
# Main map
#-------------------------------------------------------------------------------
# import raster data
use_sp()
rd <- readRAST("natura2000dist")
sg <- bb_poly(rd, projection = 3035)
# import natura2000 data
n2 <- read_sf("natura2000.gpkg")
# Aspect ratio raster data
xy <- st_bbox(sg)
asp <- (xy$ymax - xy$ymin)/(xy$xmax - xy$xmin)
# Create map
mainmap <- tm_shape(rd) +
tm_raster(style="cont", title="Distance (m)") +
tm_shape(n2) +
tm_fill(col="white", alpha=0.5, lwd=0.5) +
tm_borders(col="lightgrey", lwd=0.01) +
tm_layout(legend.outside=F, legend.text.size = 0.8,
legend.title.size=0.9,
legend.frame=TRUE, legend.position=c(0.985, 0.985),
legend.just = c("right", "top"), legend.width=-0.25,
legend.height=-0.4, outer.margins = c(0,0,0,0),
inner.margins = c(0,0,0,0))
# Inset map
#-------------------------------------------------------------------------------
# Get boundary data and project to EPSG:3035
link1 <- "https://biogeo.ucdavis.edu/data/gadm3.6/Rsf/gadm36_POL_0_sf.rds"
link2 <- "https://biogeo.ucdavis.edu/data/gadm3.6/Rsf/gadm36_POL_1_sf.rds"
PolCntr <- readRDS(file = url(link1, "rb"))
PolProv <- readRDS(file = url(link2, "rb"))
PolCntr <- st_transform(PolCntr, crs=3035)
PolProv <- st_transform(PolProv, crs=3035)
# Get bounding box and aspect ratio of Poland map
xy <- st_bbox(PolCntr)
asp2 <- (xy$xmax - xy$xmin)/(xy$ymax - xy$ymin)
# Create inset map
insetmap = tm_shape(PolCntr) + tm_fill(col="lightgrey") +
tm_shape(PolProv) + tm_borders(lwd = 1, col="grey") +
tm_shape(sg) + tm_borders(lw=2, col="red") +
tm_layout(inner.margins = c(0.04,0.04,0.04,0.04), outer.margins=c(0,0,0,0))
# Create viewport
w <- 0.25
h <- asp2 * w
vp <- viewport(x=0.985, y=0.585, width = w, height=h, just=c("right", "top"))
# Final map
#-------------------------------------------------------------------------------
tmap_save(mainmap,filename="dist2nat2000.png",
dpi=100, insets_tm=insetmap, insets_vp=vp,
height=asp*91, width=91, units="mm")