Thematic maps in R

A tutorial on how to create thematic maps for South Tyrol in R 🚀

Libraries

For now on the MonalisR package is downloaded via the dev branch. Soon the functions will be merged with the master branch. If completed simply delete the ref parameter.

#library("devtools")
#install_github("mattia6690/MonalisR",ref="dev")
library("MonalisR")
library("tmap")
library("sf")
library("GADMTools")
library("tmaptools")
library("raster")
library("grid")
library("ggplot2")
library("OpenStreetMap")
library("stars")

Background

This is a tutorial on how to potentially use the tmap Package for creating Thematic maps in R.

In this tutorial we will create a thematic map of South Tyrol added by information about the meteorological stations spread across the province.

Data

The datasets used for this tutorial have various sources:
- Meteorological Data provided by the OpenDataPortal South Tyrol and accessed through the MonalisR Package, developed at EURAC research.
- Administrative boundaries downloaded with the GADMTools Packages for Level 3 administrative areas.
- Shapefile of the administrative districs South Tyrol as delineated by the Province and used in e.g. by the weather service.
- An extent of the Mazia Valley as part of the LTSER network in South Tyrol.

South Tyrol Province

In R there are already predefined functions for importing Shapefiles on national scale. With the GADM functionalities we can load the shapefiles of individual nations in diverse levels (down to the Province scale). Here we downloaded the administrative area if Italy and extracted all Shapefiles relative to the province of South Tyrol

ItalySHP   <-GADMTools::gadm_sf_loadCountries("ITA", level=3,basefile = "data/")$sf
BzSHP      <-ItalySHP %>% dplyr::filter(NAME_2=="Bolzano")

plot(st_geometry(BzSHP))

South Tyrolean Districs

Now we can also import our own Shapefiles and display them later. Here we import a shapefile on the Districs of South Tyrol as defined by the Autonomous Province and transform the projection to LatLon coordinate systems. These administrative boundaries can be downloadable at the Geocatalogue SĂĽdtirol

BzDCT   <- sf::st_read("data/Districts_polygon.shp") %>%
  sf::st_transform(4326) %>%
  st_zm(.)
## Reading layer `Districts_polygon' from data source `C:\Users\MRossi\OneDrive - Scientific Network South Tyrol\Projects\R_Meetup_Bz\RBz_Website\BolzanoR_website\content\post\2019-09-11-thematic-maps-tmap\data\Districts_polygon.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 6 fields
## geometry type:  POLYGON
## dimension:      XYZ
## bbox:           xmin: 605877.1 ymin: 5120675 xmax: 765980.7 ymax: 5220330
## z_range:        zmin: 0 zmax: 0
## projected CRS:  ETRS89 / UTM zone 32N
ggplot(BzDCT) +
  geom_sf(aes(fill = NAME_I)) +
  theme_bw()

Meteorological Stations

The Station data for each Station can be importe in using … again … the MonalisR package. In the development version of the package the information on the meteorological stations of South Tyrol can be accessed spatially (as sf Object)

STstations <-MonalisR::getMeteoInfo(format = "spatial") %>%
  dplyr::select(SCODE,NAME_D,NAME_I,NAME_E,ALT,geometry) %>%
  unique

ggplot(STstations) +
  geom_sf(aes(fill = ALT)) +
  theme_bw()

Mazia

This is an example for the Mazia valley located in the Northwestern part of South Tyrol and central Study site of the alpine environment’s LTSER station network.
The following coce imports a shapefile and converts it to the LatLon format

MaziaSHP   <- sf::st_read("data/Mazia.shp") %>%
  sf::st_transform(4326)
## Reading layer `Mazia' from data source `C:\Users\MRossi\OneDrive - Scientific Network South Tyrol\Projects\R_Meetup_Bz\RBz_Website\BolzanoR_website\content\post\2019-09-11-thematic-maps-tmap\data\Mazia.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 621086.8 ymin: 5167146 xmax: 632606.8 ymax: 5184366
## projected CRS:  WGS 84 / UTM zone 32N

In order to plot it properly we add also a basemap to the Mazia extent. This can be done easily by donwloading the Open Street Map Data relative to the extent of the Mazia valley. The extent can be obtained using the bb function in tmaptools, Afterwards the raster is converted to LatLon and plotted as RGB image.
Important: Pay attention to the zoom level. The bigger the shapefile the bigger the zoom.

extent  <- bb(MaziaSHP)
dataras <- tmaptools::read_osm(extent,type="osm",zoom=12)

image(dataras,rgb=c(3,2,1))

Thematic maps

Analogously to the GGplot2 package the syntax of the tmap plots is based on Hadley Wickhams Grammar of graphics. Each layer is imported using the tm_shape parameter and and custiomzied by the following commands. It is possible to add all the layers necessary for thematic maps (e.g. north arrows, scale bars, grids, legends …).

In the following chapters you can find some possiblities on how to create thematic maps of South Tyrol using the tmap package.

Polygons

Polygons are best added using sf objects. Here we added the South Tyrolean districts together with the extent of the Mazia valley and the names of each district both in italian and in German language.

tmraw<-tm_shape(BzDCT)+
    tm_grid(col = "grey90",ticks = T)+
    tm_fill("grey95")+
    tm_borders("grey")+
    tm_compass(type = "8star", position = c("right", "top")) +
    tm_scale_bar(breaks = c(0, 20, 40), size = .8,position = "left")+
    tm_text("NAME_I",size=.5,ymod=-.5)+
    tm_xlab("Longitude")+
    tm_ylab("Latitude")+
    tm_layout(legend.bg.color="white",legend.frame = T)+
  tm_shape(BzDCT)+
    tm_text("NAME_D",size=.5,ymod=-1)+
  tm_shape(MaziaSHP)+
    tm_borders("green",lty=2)+
    tm_add_legend(labels="Mazia/Matsch Valley Extent",col="white",border.col="green")

tmraw

Polygons and Points

In order to add some point we can just add futher layers to the map. Here we added all the meteorological stations within the extent of South Tyrol to the tmraw map we generated previously. In the legend we evidenced the altitudes of each station continuously and added the desired breaks manually

tmstat<-tmraw +
  tm_shape(STstations)+
  tm_symbols("ALT",size=1,
             breaks = seq(250,2250,500),
             title.col ="Meteorological Stations (Altitude)",
             palette="Blues")

tmstat

Polygons and Raster

With the tm_rgb function it is possible to combine both a raster and Polygons on one map. Here we combined the downloaded OpenStreetMap image of the Mazia valley with the actual extent of the map and added all the necessary layers needed for a thematic map. This can be further customized by adding tmap elements.

tmmatsch<-tm_shape(dataras) +
    tm_rgb() +
    tm_grid(ticks = T,alpha = 0)+
    tm_compass(type = "8star", position = c("right", "bottom")) +
    tm_scale_bar(breaks = c(0, 2.5, 5), size = .8,position = "right")+
    tm_layout(legend.bg.color = "white",legend.frame = T)+
    tm_xlab("Longitude")+
  tm_shape(MaziaSHP)+
    tm_borders(col="red")+
  tm_add_legend(labels="Mazia/Matsch Valley Extent",col = "white",border.col="red")

tmmatsch

Combine Maps

We can combine two different datasets in one map. Therefore we have to do them separately and combine them in one plot later on. Oftentimes this approach is useful for delimiting the actual study site within a broader scale.
In this example we will create a map of the Vinschgau area including the Mazia valley together with the meteorological station that can be found there. Additionally we will create an inset of South Tyrol indicating where the Vinschgau valley is located within the study site.

The map of the Study site is done simply by adding the layer. We added no further details since it is only thought to roughly demostrate where we are

vinschgau<-BzDCT %>% dplyr::filter(NAME_D=="Vinschgau")

tmview<-tm_shape(BzDCT)+
  tm_fill("grey95")+
  tm_shape(vinschgau)+
  tm_borders("orange")+
  tm_shape(MaziaSHP)+
  tm_borders("green",lty=2)

tmview

Now we first select the meteorological stations in the Vinschgau region and later create a dedicated plot for the region.

STstations.vi<-st_intersection(STstations,vinschgau)

tmvi<-tm_shape(vinschgau)+
  tm_grid(col = "grey90",ticks = T)+
  tm_fill("grey95")+
  tm_borders("orange")+
  tm_compass(type = "8star", position = c("left", "bottom")) +
  tm_scale_bar(breaks = c(0, 5, 10), size = .8,position = "left")+
  tm_xlab("Longitude")+
  tm_ylab("Latitude")+
  tm_layout(legend.bg.color="white",legend.frame = T)+
  tm_add_legend(labels="Venosta/Vinschgau Valley Extent",col="white",border.col="orange")+
  tm_shape(MaziaSHP)+
  tm_borders("green",lty=2)+
  tm_add_legend(labels="Mazia/Matsch Valley Extent",col="white",border.col="green")+
  tm_shape(STstations.vi)+
  tm_symbols("ALT",size=1,
             breaks = seq(250,2250,500),
             title.col ="Meteorological Stations (Altitude)",
             palette="Blues")+
  tm_text("SCODE",ymod=1)

tmvi

And now we combine both images and export them. I recommend to combine AND save them directly. Otherwise you may have problems to completely understand the position of the viewports

tmap_save(tmvi,filename = "Vinschgau.jpg",width = 2000,height = 2000,
          insets_tm = tmview,
          insets_vp = grid::viewport(0.83, 0.86, width = 0.35, height = 0.3))

Further Information

There are a lot more posibilities to use tmap. This is just a bief overview on the potential of the package. Find more information consider the following links:
- Tmap tutorial page.
- Geocompr book chapter
- Rpubs presentation
- A shiny Tmap app
- The Publication in IJSS
- A Presentation at OpenGeo MĂĽnster

Avatar
Mattia Rossi
Institute for Earth Observation

My Research Interests comprise Spatial and EO Data Processing, Statistical Computation with R.

Related

comments powered by Disqus