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")
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.
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()
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()
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))
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 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
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
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
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))
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