USING R PACKAGES ‘TMAP’, ‘RASTER’ AND ‘GGMAP’ FOR CARTOGRAPHIC VISUALIZATION: AN EXAMPLE OF DEM-BASED TERRAIN MODELLING OF ITALY, APENNINE PENINSULA

: The main purpose of this article is to present the use of R programming language in cartographic visualization demonstrating using machine learning methods in geographic education. Current trends in education technologies are largely influenced by the possibilities of distance-learning, e-learning and self-learning. In view of this, the main tendencies in modern geographic education include active use of open source GIS and publicly available free geospatial datasets that can be used by students for cartographic exercises, data visualization and mapping, both at intermediate and advanced levels. This paper contributes to the development of these methods and is fully based on the datasets and tools available for every student: the R programming language and the free open source datasets. The case study demonstrated in this paper show the examples of both physical geographic mapping (geomorphology) and socio-economic geography (regional mapping) which can be used in the classes and in self-learning. The objective of this research includes geomorphological modelling of the terrain relief in Italy and regional mapping. The data include DEM SRTM90 and datasets on regional borders of Italy embedded in R packages ‘maps’ and ‘mapdata’. Modelling references to the characteristics of slope, aspect, hillshade and elevation, their visualization using R packages: ‘raster’ and ‘tmap’. Regional mapping of Italy was made using main package ‘ggmap’ with the ‘ggplot2’ as a wrapper. The


Introduction
The main tendencies in geographic education have been influenced by the open source GIS and publicly available free geospatial datasets (both raster and vector formats) that can be used by students for cartographic exercises. In view of the increased demand on the distance-based self learning creates the need for the open source programs in educational technologies. The development of the open source geoinformation technologies in cartographic methods recently made significant progress due to such excellent open source GIS as cartographic toolsets as QGIS, Generic Mapping Tools (Lemenkova, 2020a(Lemenkova, , 2020b, ILWIS GIS (Lemenkova, 2013), IDRISI GIS (Da Serra Costa et al. 1996;Lemenkova, 2014), GRASS GIS (Alvioli et al. 2020;Lemenkova, 2020c), SAGA GIS (Vacca et al. 2014;Lemenkova, 2020d), largely used in geoscience research. Besides the software, a variety of open source data exist available from various sources: Landsat TM and Sentinel-2A satellite images by the USGS, widely used in geoscience (Allevato et al. 2019;Lemenkova, 2015), Google Earth aerial imagery, Digital Chart of the World (DCW), raster topographic datasets: elevation DEMs, such as SRTM15, GEBCO, ETOPO1, GLOBE, to mention a few. Therefore, using available sources of such geographic information, students can get free of charge datasets, vector layers, raster grids and satellite images that are constantly updated.
Once the data are captured, another question arises: which software to use? Leaving apart the mentioned above GIS applications, this paper focuses on the presenting non-trivial method of geographic data visualization in geographic education: the use of free and open source R programming language (R Core Team, 2020) for cartographic visualization by students. In particular, its specific libraries (also known as packages) used in this paper are the following ones: ggmap (Kahle & Wickham, 2013), maps, mapdata, sp, raster, tmap. These libraries are specifically tailored for cartographic visualization. In turn, they require the dependent packages, used for auxiliary data manipulation: ncdf4, RColorBrewer (Brewer et al. 2003;Neuwirth, 2014), sp (package for processing spatial data) and sf (Pebesma, 2017).
For a long time, GIS was the only possibility to create maps, while the R language and its libraries was used to obtain the statistical information for modelling, assessment and graphical visualization of the tabular datasets (Lemenkova, 2019a). However, a variety of R packages nowadays presents new opportunities of using machine learning approaches in geographic education. Using R, quantitative and qualitative types of data stored in the tables as 'data.frame' are processed using a variety of statistical methods, and visualized as plots. However, it is also possible to use R not only for statistical analysis but in geographic mapping, which is much less known and popular comparing to GIS spatial applications in geographic and environmental data analysis (e.g. Allen et al. 1999;Suetova et al. 2005;Klaučo et al. 2013aKlaučo et al. , 2013bLemenkova, 2011). Specifically, R can be used for plotting traditional maps, and receiving information from the more complex spatial geomorphological analysis (slope, aspect, elevation, hillshade). The geomorphometry is widely used and applied in landform mapping (Evans, 2012) which makes R especially useful in landscape studies.
R has a relatively straightforward syntax with a variety of packages that can be installed and uploaded additionally, via the RStudio (RStudio Team, 2017) and package installer by R. For instance, using a module 'tm_graticules' can be used for adding grid graticule and its additional elements (ticks, lines, labels) on maps and plotting them in a specific way (e.g. rotated, colored, using enlarged fonts, etc). R package 'tmap' can also be used for creating complex cartographic layouts with legends. Applying various color palettes can be done via the package 'RColorBrewer', adding cartographic elements (legends, multi-level annotations, scale bars, directional compass). The use of 'ncdf4' package is an additional condition for manipulating with various raster formats. The 'tmap' gives an opportunity to create classic maps in R, while 'raster' packages can be used for geomorphological analysis (calculation of slope, aspect, hillshade and elevation) from the SRTM 90 m datasets available by the getData() function, which gives an access to the elevation data.
The research objective of this work concerns the following issues: 1) Introduce application of the cartographic functionality of R available in the packages 'tmap' (Tennekes, 2018), 'raster' (Hijmans, 2017), 'maps' (Becker et al. 2013), 'mapdata'. 2) Perform geomorphometric analysis of the Apennine Peninsula by 'raster' package with a visualized screenshot of R script and final layout output. 3) Present technical detailes of the packages as a visualized codes. 4) Apply 'tmap' and 'ggmap' packages for cartographic visualization. The presented topics should cover a large need of the cartographic issues in geographic education, since both R and data used in this study are open source and can be used in distance-based learning and geographic education.

Study Area
This paper focuses on mapping the Italy with two specific examples of geographic mapping: physical geographic and socio-economic (regional): 1) topographic mapping based on DEM, using 'tmap' and ''raster packages as the main tools. This includes calculation and visualization of slope, aspect, hillshade and elevation using SRTM90 dataset; 2) regional mapping of country borders using 'ggmap' packages as the main tool, 'maps' package (Becker et al. 2013) and a 'mapdata' as its supplement, providing access to larger and higher-resolution databases. Italy has been selected as the country with one of the most representative countries on geomorphology and hilly relief, suitable for terrain analysis.
The landscapes on the Apennine Peninsula evolve under the control of a variety of factors which should be briefly mentioned. Among the most important ones are complex geological setting, fault system, tectonic deformations creating well-exposed faults in the region (Roberts & Michetti). For instance, the uplift occurring during Middle and Upper Pleistocene caused formation of the low relief. The uplift of the Northern Apennine mountain chain shaped the fluvio-lacustrine landforms in Italy (Bartolini, 2003). Landscape formation is also affected by vegetation, climate dynamics and sedimentation (Bertotti et al. 1997;Guido et al. 2020). Surface topography reflects landscape evolution, geomorphic processes, which in turn are influenced by the tectonic events (Bull, 2007), seismicity and earthquakes in faulted mountain fronts affecting river network of the basins and highs. Other factors sculpturing landforms include formation of the orogenic belts, morphotectonics of the Pliocene-Pleistocene regional uplifts affecting topography of the Central-Northern Apennines and Po Basin (Centamore, Nisio, 2003;Zuffetti & Bersezio, 2020) and tectonically active L'Aquila Basin (Cosentino et al. 2017).
As a result, modern topography of the Apennines reflects the interactions between crustal-mantle and surface processes since the Late Miocene which affect the geomorphology of the peninsula in course of geologic, and tectonic-morphologic evolution of the orogenic belt in the Apennines . The formed landscapes record the sequence of geomorphic evolutionary steps, as well as the spatial-temporal relations between the active geological processes. As a result, the relief of Italy is highly uneven and varies regionally within the country, which makes it suitable for topographic visualization, already presented in existing papers on terrain DEM analysis of the selected landscapes in Italy (Tarquini et al. 2007;Ascione et al. 2008;Geurts et al. 2020;Gemelli et al. 2011).

Materials and methods
There are a wide variety of packages in R language used for graphical plotting. This research uses the packages 'ggmap', 'maps', 'mapdata', 'sp', 'raster', 'tmap' for data capture, processing and cartographic visualization (Figure 1). Using these packages R demonstrates the functionality of traditional GIS software, briefly described below.
The collected data were pre-processed before they were used as an input data in R. Thus, when the data were loaded, the extent (dimensions) of the tabular data was checked up, as well as the column content: dim(italy). Since the table is rather big-sized, the initial and end part of the tables were inspected using the following functions: head(italy) and tail(italy), respectively. Then the data were inspected for their coordinate projection: crs(alt).
Afterwards, the common geographical projection (lon/lat) was defined using the following code: crs(alt) <-"+proj=longlat +datum=WGS84 +no_defs". If necessary for the students to experiment with projections, the variety of the projections can be set up using this function. For instance, converting geodetic coordinate to UTM Zone 33 on the northern hemisphere (EPSG:32633 WGS 84/UTM zone 33N) can be done using following code: crs(alt) <-"+proj=utm +zone=33". Afterwards, the spatial extent of the data was defined using the following code: e <-as(extent (6,19,36,48), 'SpatialPolygons'). Here the West-East-South-North system is defined as the borders fo the study area ( Figure 2). The digital elevation models (DEM) obtained from the Shuttle Radar Topography Mission (SRTM) create perfect dataset for geomorphological studies due to the open availability and high resolution of the grid. There are various examples of using the SRTM data in geospatial research (Drăguţ & Eisank, 2012;Hirt, 2018) with discussed reliability (Frey & Paul, 2012). The SRTM data sources present a reliable and accurate topographic dataset to gather reliable information about the relief terrain, taking into account the elevation values. The spatial SRTM DEM used in this research have 90m*90m resolution.

Calculate terrain characteristics: slope and aspect
A complex geomorphic analysis may include computation of indices such as mountain front of water divide sinuosity, asymmetry, drainage basin elongation, relief ratio, hypsometry, normalized steepness, and concavity, integrated with geomorphological analysis (Giano et al., 2018). This study presents the most important geomorphic features that are often included in geographic education curriculum and also embedded in R algorithms: computation and mapping of slope steepness, aspect orientation, visualization hillshade and DEM elevation. The slope stability analysis is based on calculating the slope steepness often applied in studies on geologic hazard risks, such as assessment of landslides in mountains, rock avalanches (Antonielli et al. 2020;Lemenkova et al. 2012).
Calculating slope, aspect and hillshade is based using existing slope and aspect generation mathematical algorithms (Ritter, 1987) of geomorphological models. Other computational examples of the terrain include, for instance, landscape metrics, used in environmental and sustainability studies assessing the vulnerability and ecological significance of the landscapes (Klaučo et al. 2014(Klaučo et al. , 2017. In R, the calculation of slope and aspect maps was done using the following codes: slope = terrain(alt, opt = "slope") and aspect = terrain(alt, opt = "aspect"), respectively ( Figure 2).
After the slope and aspect were computed using the scripts presented in Figure 2, the preliminary visualization without cartographic elements was done using the plot(slope) and plot(aspect) functions, respectively. Afterwards, the response of R for calculated slope and aspect unit analysis was translated to the 'tmap' package as the spatial input data where the additional cartographic elements were added on the maps (graticule, grids, annotations, title, subtitle, legend, etc) using multiple specially designed functions, as follows: tm_scale_bar, tm_compass, tm_raster, tm_layout, tm_graticules. The full script used for plotting the maps of slope and aspect is visualized in Figure 2. The overlay of these cartographic elements, slope and aspect maps, their legend and histograms of data distribution was performed after the slope and aspect maps were prepared based on the SRTM DEM data of Italy. The final maps are presented in Figure 3 (left: slope, right: aspect) and Figure 4 (left: hillshade, right: elevation).

Visualizing topographic maps: elevation and hillshade
The elevation DEM data of SRTM was used to create a hillshade map. The raster package provides an easy access to the SRTM 90 m resolution elevation data with the getData() function. For example, it is possible to download the elevation data for the whole country using the getData() function. The hillshade map shows the topographical shape of highs, hills and mountains using levels of gray on a map. The role of this kind of maps is to display relative slopes, though not in absolute heights. However, despite its relative estimation of the topography, the value of the shaded relief maps consists in its approach: they give the student an immediate appreciation for the surface topography, because a hillshade map helps to visually estimate the depths (Horn, 1981).
The hillshade map was created based on the terrain characteristics using the 'raster' package and then visualized as classic map with elements using 'tmap' Figure 2.

Figure 3. Map of slope (left) and aspect (right) plotted in R using scripts presented in
package. Thus, initially, the hillshade was calculated based on the input data from the previous research steps: slope and aspect using algorithms of hillshade computing (Jones, 1998). Both of them were calculated using 'raster package' with the terrain function. Computing of a hillshade was done by the hillShade() function using the following code: hill = hillShade(slope, aspect, angle = 40, direction = 270). Now the hillshade was computed using slope and aspect of the terrain with the option argument set to "slope" or "aspect", respectively. The previously created slope and aspect objects were used as input data and two new arguments -angle as 40° (the elevation angle of the light source) and direction as 270° (he direction (azimuth) angle of the light source (sun)) -were set up. The output hillshade map was plotted using plot(hill) function. Once the hillshade was computed, it was transferred to the 'tmap' package for further cartographic plotting. Here the visualizing was done using stacking of data layers similar to the traditional overlay in GIS, with the layer of the hillshade object (hill) colored using grey hues.

Figure 4. Map of hillshade in monochrome colors (left) and elevation map
in terrain colors (right), plotted in R.

Mapping regional division
The advantages of the 'ggmap' package consists in the convenient mapping and plotting the regional map of Italy ( Figure 5 and Figure 6). The maps package contains the outlines of Italy as country that have been with R as embedded data. The 'mapdata' package was used, since it has a higher-resolution vector outlines. The 'maps' package was used for its plotting functionality in addition to the 'ggplot2 package which operates on data frames (Wickham, 2009). Therefore the data from the 'maps' package were forwarded into a data.frame format of the 'ggplot' package ( Figure 5, left).
The country data (Italy) and SRTM data (DEM) were also collected and loaded to RStudeio using common R syntax (Bivand et al. 2008). The country data on Italy origins from the NUTS III (Tertiary Administrative Units of the European Community) database of the United Nations Environment Program (UNEP) GRID-Geneva datasets made in 1989 (Becker & Wilks, 1995).  ', 'ggmap', 'maps' and 'mapdata'. Script used for visualization (left) and resulting output map (right).
The function 'tm_shape' of library 'tmap' was applied for the creation of cartographic layout and then saved in jpg format using R. The 'tmap' library defines the shape objects plotted on a map (Tennekes, 2018) and visualized in Figure 6. The base graphics of R was head directly to ggplot2 package. The 'gg-plot2' package was used as a wrapper for the raster map visualization, enabling the finer operation with spatial objects. Using some functions of ggplot2 enables easier interaction with the data in R 'maps' package. Here the main map was plotted using geom_polygon function (geom_polygon(data = italy, aes(x=long, y = lat, group = group)), and aesthetics was added using additional elements of ggplot2 syntax: fill = "pink", color = "blue", linetype = 1, size = 0.2).
The annotations of the XY axes were made using the options xlab and ylab: xlab("Longitude") + ylab("Latitude"). The titles and the subtitles were also added on a mp using the following arguments of the 'ggplot2' package: labs(ti-tle="Italy", subtitle = "Mapping: R", caption = "Packages: ggmap, ggplot2, mapdata, maps"). The legend has been added on a map using the code: guides(fill = guide_legend(reverse=TRUE)). The disadvantage of the 'ggplot2', however that it does not support different cartographic projections. Once a map layout is prepared, the map was converted to the 'jpg' format and visualized in Figure 6.

Results and discussion
The visualized DEM derivatives included geomorphometric maps of slope, aspect and hillshade relief maps which shows the trends in the orientation (North-South-West-East) and slope steepness of the raster grids in different sub-regions of the Apennines. The slope directions of the rural areas of Italy revealed variations in the data showing 'gentle', 'moderate', 'strong', 'very strong', 'extreme' and 'steep' slopes of the mountainous regions of Italy. The number of pixels (over 75,000 on a raster grid) was the greatest in the 'moderate' slope level comparing to the others: 35,000 pixels for the 'gentle', 22,000 for 'strong', 18,000 for 'very strong', 11,000 for 'extreme' and 9,000 for 'extreme' slopes. As the influence of the geomorphological patterns was clear on the slope and aspect maps correlating with a general elevation map layer, detecting geomorphometric trends shows effective results by the presented and used R packages.
The trends in aspect data distribution in the Apennines visualized using R package 'tmap' apparently reflect the structure of the geomorphological patterns of the study area. When comparing the four major classes of the slopes orientation (West-North-East-South) with the subdivided (North-West, North-East, South-West, South-East) in Figure 3, the integrated dataset shows four classes: 1 st class (East, colored red), the 2 nd class (West, colored orange) 3 rd class (South, colored yellow) and the 4 th class (North, colored blue). The data assessment demonstrated the subdivision of the slopes by DEM. The visualized elevation map of the Apennines identified heights in meters grouped in five major classes ( Figure 4) and colored by terrain.colors palette of R. The classes include heights from 0 to 4,000 meters. Additionally, a hillshade map was visualized as a raster calculation based on slope and aspect layers (Figure 4).
The applied parameters of R libraries were sufficient for the terrain mapping and geomorphological analysis, as presented in this work, and therefore, R can be strongly recommended in similar research studies. All maps presented in this article are made and plotted in R. Although there are special cartographic toolsets and software has functionality for the complex cartographic workflow using scripting approaches e.g. GRASS GIS or GMT (Lemenkova, 2020e, 2020f), the direct use of R language presents the principally next steps in the development of the cartographic methods, since using programming language is completely based on the machine learning algorithms for data evaluation, which could be effectively used for geoinformatics and mapping.

Conclusions
Automatization in topographic modelling is largely presented by the use of scripting languages (Lemenkova, 2019b(Lemenkova, , 2019c(Lemenkova, , 2019d or specially designed software that uses approaches of machine learning in geoinformatics (Gauger et al. 2007;Iwahashi & Pike, 2007;Schenke & Lemenkova, 2008;Alvioli et al. 2016;Kuhn et al. 2007). Both practical and theoretical conceptual approaches are the core issues of the GIS-based geomorphic mapping. Therefore, the advances in R programming and development of its new libraries present new ways in geographic education by rapidly evolving scripting technologies used in cartographic visualization. The extended possibilities of R can successfully support and top up the existing GIS methods. Besides, free and open source R language is available for every student, which makes it especially useful in geographic education due to its availability and functionality.
However, nowadays, the applications of R are mostly presented in cases of the statistical analysis (Lemenkova, 2018, while the use of its cartographic packages for geographic visualization have not yet been adequately addressed. This paper aims to fill in this lacuna and to present the illustration of how R libraries designed for cartographic needs can be used in classes of geography by the students for geomorphological analysis with a case study of Italy, known for the diverse relief patterns (mountains, slopes, hills, plains, highs). Mapping approaches by R presented in this research demonstrated maps of the terrain analysis based on DEM of the Apennines.
The value of using R in geosciences consists in practical reproducibility of scripts in geography and geosciences. Therefore, the use of R in geographic education is a highly promising approach. However, there are some possible drawbacks and disadvantages comparing to the traditional GIS, which should be mentioned as well. Thus, the syntax of R may not always easily accessible by the beginner students and require preliminary skills in coding. Second, the presented packages require installation. Third, some concepts of R may not be always straightforward for the students and require additional classes and special teaching on script approaches. Fourth, the plotting of maps always require some familiarity with general GIS concepts, such as projections, plotting, operating with spatial data, etc.