Quantitative Morphometric 3D Terrain Analysis of Japan Using Scripts of GMT and R

Abstract

In this paper, we describe two related scripting methods of cartographic data processing and visualization that provide 2D and 3D mapping of Japan with different algorithm complexity. The first algorithm utilizes Generic Mapping Toolset (GMT), which is known as an advanced console-based program for spatial data processing. The modules of GMT combine the functionality of scripting with the aspects of geoinformatics, which is especially effective for the rapid analysis of large geospatial datasets, multi-format data processing, and mapping in 2D and 3D modes. The second algorithm presents the use of the R programming language for cartographic visualization and spatial analysis. This R method utilizes the packages ‘tmap’, ‘raster’, ‘maps’, and ‘mapdata’ to model the morphometric elements of the Japanese archipelago, such as slope, aspect, hillshade and elevation. The general purpose graphical package ‘ggplot2’ of R was used for mapping the prefectures of Japan. The two scripting approaches demonstrated an established correspondence between the programming languages and cartography determined with the use of scripts for data processing. They outperform several well-known and state-of-the-art GIS methods for mapping due to their high automation of data processing. Cartography has largely reflected recent advances in data science, the rapid development of scripting languages, and transfer in the approaches of data processing. This extends to the shift from the traditional GIS to programming languages. As a response to these new challenges, we demonstrated in this paper the advantages of using scripts in mapping, which consist of repeatability and the flexible applicability of scripts in similar works.

1. Introduction

Cartographic visualization is an important component of many Earth science applications. In numeric land modelling, the discriminative power of 2D and 3D mapping is a key factor in the topographic analysis of land features since it provides the most direct and quickest way to evaluate geospatial information. During the last four decades, a variety of Geographic Information Systems (GIS) have been developed, with commercial ArcGIS software certainly being the most widely used. Land surface modelling in GIS shows that qualitative data analysis, cartographic visualization and the interpretation of the topographic features visualized on maps activate the evaluation of spatial heterogeneity and the variability of the objects and environmental processes on the Earth.

The performance of various GIS is, in many cases, tailored to their specific tasks, among which vector and raster spatial data processing is arguably the most prominent and important functionality of these software. For instance, some are better suited for image processing, such as Erdas Imagine [1], Idrisi GIS [2,3], Integrated Land and Water Information System (ILWIS) GIS [4], and ENvironment for Visualizing Images (ENVI) GIS, while others are best at vector data analysis and visualization, e.g., ArcGIS [5,6,7,8,9]. Most of these GIS are based on a standard interface with a restricted functionality that requires the manual processing of data, although recently machine learning techniques have been applied to spatial data modelling [10] and image analysis [11,12]. However, while these methods drive the geospatial analysis to gather explicit information on the terrain structure and morphometric variations, the core question emerges as the cartographic approach underlying the computational complexity of geoinformation processing. This raises the multi-disciplinary goal of succeeding in improving the cartographic workflow over the performance of the conventional algorithms. For instance, the use of scripts, besides the pre-defined algorithms in GIS, can result in an enlarged cartographic workflow functionality for data being modelled.

GIS-based mapping is a tedious and time-consuming process for cartographic performance as it may involve a large number of separated tasks and operations that are normally made using different commands in the menu toolbar. For this reason, the workflow is split into various steps of data processing. Better still, however, is using the full functionality of the programming and machine learning, applications of scripting languages for plugins, and auxiliary tools that ensure data processing using scripts [13]. Spatial analysis can benefit from the automation of cartographic processes because scripting enables the repeatability of the process. The advanced modelling enables the performance of a more comprehensive analysis of various factors affecting land surfaces and processes.

On the other hand, the importance of machine learning for topographic and geomorphological mapping is well known since it is less error-prone and time-consuming compared to the traditional state-of-the-art GIS. At the very least, this can be a use for scripting for data visualization, such as in Generic Mapping Tools (GMT) [14,15]. Despite certain difficulties in mastering the program, such as its high learning curve, console-based non-visual operation mode, and complex steps of fine-tuning map elements, GMT nevertheless provides a much more powerful functionality of cartographic workflow and increased flexibility in data processing. Furthermore, it yields deep insights into how to process geospatial information in relation to the individual traits of datasets, including the transformation of coordinate systems and georeferencing, extracting attributes and labelling, processing binary formats for gridded datasets, modelling elevation data for morphometric analysis, reading interleaved data formats in one script, and many more.

In this paper, we extend this idea of using scripts to the morphometric mapping of Japan by incorporating several libraries of R and modules of GMT for cartographic data processing. We propose a scripting framework that, by facilitating the mapping process using powerful methods of programming, can plot 2D and 3D maps more effectively than standard techniques of GIS. We devise a systemic way to use and process geoinformation derived from raster grids for visualising morphometric elements and parameters of terrain and use them to develop a series of scripts that could be extended to other regions to improve cartographic performance in 2D and 3D modelling.

The accurate visualization of complex terrain models can lead to the increased precision and efficacy of maps. Cartographic layout is usually seen as the final stage of spatial analysis summarising and presenting the results in a graphical form. That said, scripting by R land GMT largely facilitates visualization since similar parts of scripts may be reused. The contribution of the paper is twofold. First, we demonstrated that script-based visualization offers the increased precision of the morphological analysis through automation of data processing. Second, a marked improvement in the development of cartographic methods is observed, since using R and GMT in a mapping workflow paved the way for the machine learning methods in geospatial data analysis.

1.1. Background and Motivation

Topography is an important characteristic of the land surface, which is reflected in the landscape and associated with numerous environmental factors. Specifically, the geomorphic and topographic features control and affect vegetation coverage [16,17], habitats and landforms [18,19], hydrology, soil distribution, and local micro-climate settings that are dependent on the topographic exposure, slope steepness, and curvature of the relief. Spatial characteristics that can be retrieved from the topographic and land surface maps have many applications, such as environmental management [20], flood monitoring [21], hydrological and fluvial modelling in riverine ecosystems [22,23], analysis of coastal processes [24], soil management practices, and analysis of crop yields [25,26,27]. Moreover, land surface maps present a background for the analysis of the correlation between the geomorphology and geology [28] and geophysical and geographic factors and processes [29].

Topographic mapping based on the digital elevation model (DEM) is considered one of the most important issues of cartography and is often used as background information for spatial analysis in the geosciences. An inherent research step in the computations of a terrain analysis based on DEM is the use of modelling methods where GIS is traditionally applied. A variety of existing GIS software can be used for the spatial data visualization and topographic analysis considered a background for both socio-geographic and physical-geographic mapping. Since the onset of the development of GIS, technical methods of cartographic visualization and approaches to terrain analysis have been constantly improved.

One of the key components in GIS software and cartographic data processing is map projection. Mapping lands evokes a distortion in their angle, area, and shape depending on the type of map projection, which may be of the conic, azimuthal, cylindrical or miscellaneous types [30]. As a result, minor distortions arise when transforming the coordinates into various types of projections. Therefore, various projections can be better adjusted to map specific study areas depending on their locations and spatial extent.

For instance, the Albers conic equal-area projection is mostly used to plot areas with large longitudinal extents, e.g., Canada and USA. The polyconic projections have true scales on the parallels represented as non-concentric circular arcs. This makes this projection class neither equal-area nor conformal with the least distortion on the parallels having their centres along a central meridian. Likewise, the Lambert conic conformal projection is suitable for regions with a W-E extent with a true scale on the two standard parallels. The equidistant conic projection keeps a balance between the conformal and equal-area types with minimised distortion over the study area and a true scale along all the meridians and standard parallels.

Large regions with a global extent are better mapped using the Lambert azimuthal equal-area or stereographic equal angle projections, while northern regions, such as Scandinavian areas, can be effectively plotted using the polar stereographic projection, where map boundaries are represented by lines of constant longitude and latitude. In this study, we used the cylindrical Mercator projection for plotting a topographic map that is conformal in type with an inserted small global map in a perspective projection.

1.2. Related Work

Various publications focused on Japan were published that, in particular, aimed to understand how the geomorphology of the Japanese Alps is linked to the geologic and surface structures [31], rock glacial processes and distribution of vegetation [32], or reflected in the adjacent bathymetry [33]. Multivariate methods of spatial analysis extract local or regional terrain features from spatial data using coordinates that locate them. Thus, a variety of analytic techniques involves methods of quantitative analysis for detecting and describing spatiotemporal information for a deeper understanding of land surface processes. These can be used as the advanced tools enabling the performance of a terrain analysis on the highly heterogeneous morphological setting of Japan, as reflected in relevant works [34,35,36,37,38].

By choosing a proper GIS, data handling can yield a large number of extended features. However, providing a standardised solution for selecting a proper GIS is a challenge since, in most cases, it has a similar functionality of both raster and vector data analysis. Nevertheless, whichever GIS is chosen, it requires active work in an interface menu with the manipulation of input data and modelling to present a cartographic visualization. Another point is the issue of open source availability, which has restricted access for commercial GIS. The increasing variety of GIS necessitates that alternative tools with more flexible and varied approaches to cartographic visualization and geospatial modelling should be explored. Scripting approaches can provide a substitute functional tool for cartographic visualization and geospatial modelling. Examples of scripting toolsets in cartography include the Geographic Resources Analysis Support System (GRASS) GIS [39,40] and GMT used in topographic and geomorphic mapping [41,42].

1.3. Contribution

This paper presents an application of the GMT cartographic scripting toolset [43] and the R language [44] for the spatial analysis and terrain mapping of Japan. We demonstrate the main advantages of both methods, which consist of their straightforward and logical language syntax, a scripting approach, and open-source access. We use the powerful functionality of several GMT modules for the fine-tuning of the cartographic visualization, as well as R packages to process both tabular and geospatial data directly from a console. We considered the existing cases of using R libraries for geospatial data processing and visualisation. For instance, special packages for geographic data processing include gstat [45], RStoolbox [46], terra, and raster [47].

This study makes a technical contribution to the development of cartographic methods by adapting the R language for mapping instead of the traditional GIS. We combined the GMT scripting toolset and R language for spatial data visualization and morphometric analysis in a cartographic framework, which included following general steps: (1) processing the data of the General Bathymetric Chart of the Oceans (GEBCO), Earth Topography Global Relief Model (ETOPO1) and ETOPO2 by GMT for 2D and 3D mapping by scripts; (2) the importing of Shuttle Radar Topography Mission (SRTM90) DEM by R; (3) visualising a DEM using the ‘raster’ package; (4) computing the topographic surface parameters of slope, aspect, and hillshade using R scripts; and (5) cartographic visualization by the ‘tmap’ package of R with additional cartographic elements (legend, histogram, grid, ticks, annotations, compass directions, and a bar scale).

The structure of this document is as follows. First, we discuss the key issues relevant to the morphometric modelling based on DEM in a cartographic analysis of the geospatial data. Specifically, slope, aspect, hillshade and elevation are reviewed as essential elements of the terrain analysis. Second, we present a case of the Japanese archipelago as an example for mapping and briefly outline the major geographic features and morphology of the Japanese Alps. Third, we point out that our combination framework can be easily extended to other regions with heterogeneous terrain morphology with minor modifications of the scripts. To this end, we provide technical notes on scripts, followed by a methodological description of the performed morphometric analysis using R and GMT. The full scripts are listed in the Appendix as a technical cartographic reference for similar studies. Furthermore, the methodological part demonstrates several screenshots of the performed scripting process in the RStudio environment, where maps of slope, aspect, hillshade and elevation are plotted, and the GMT for 2D and 3D modelling.

The results present the computed and visualized maps of the Japanese archipelago made by a GMT scripting toolset and in RStudio. The latter ones include the morphometric elements (slope, aspect, hillshade, and elevation) with statistical elements of data distribution as histograms and a map of prefectures of Japan visualized by the ‘ggplot2’ graphical package of R. A discussion of the machine-based approach to cartographic modelling and visualization follows the Results section. Finally, we provide conclusions regarding the advantages of using GMT and R in cartography compared to the state-of-the-art GIS techniques based on a standard graphical user’s interface and ways of modelling geospatial data by scripting languages from the console.

2. Materials and Methods

The key methods included the GMT with diverse modules and packages of Rm, including ‘tmap’ for thematic mapping and ‘ggplot2’, ‘ggmap’ [48], ‘maps’ and ‘mapdata’ for the data capture and thematic regional mapping of prefectures of Japan. The package ‘mapdata’ contains the base location of the binary files of prefectures boundaries of Japan used by the map drawing functions. The ‘raster’ package was used for morphometric analysis, which also included the depending packages, such as ‘sp’, a package which operates with classes and methods for spatial data, and ‘sf’, a package operating with ‘Simple Features’ as object classes in an R syntax.

2.1. Study Region

The study region is Japan, shown in Figure 1.

The geology of Japan is characterised by an early stage of mountain range formation comprising young and active island arcs [49]. As a result of complex tectonic movements, uplift, and denudation processes, the Japanese Alps consist of a mountain range cut by river and glacial valleys, running through the main island of Honshu. During the Quaternary period, the major axis of the Japan Alps experienced more than 2000 m of uplift, which, together with the denudation processes, resulted in the exposure of the Takidani Granodiorite. The regions is divided into three major parts: the Northern, Central, and Southern Alps. The northern region located in the Nagano, Toyama and Gifu prefectures includes the Hida Mountains with high seismicity along one of major Quaternary faults, the Atotsugawa fault, with repetitive earthquakes [50].

The geologic structure of the Hida terrane shows a complex geologic structure with gneisses formed from Permo-Carboniferous clastic sediments during a single metamorphism at ca. 250 Ma [51]. The Northern Japanese Alps are characterised by lateral and terminal moraines and outwash terraces of glacial origin. The Central Alps include the Kiso Mountains located in the Nagano prefecture [52] and are characterised by a granite structure. The highest peak of the Japanese Alps is recorded in Mt. Kita (3193 m) in the Southern Alps, or Akaishi Mountains.

The Japanese Alps had long been exploited, and their geologic structure, geomorphologic setting [53], and tectonic evolution [54] have been reconstructed. In addition to the environmental value and strong effects on the climate and vegetation setting of Japan, the Japanese Alps contribute to the economy of the country, being a source of natural mineral resources such as timber and minerals, as well as habitat for diverse species and vegetation, including medical herbs. Moreover, rice paddy fields are cultivated on the slopes of the mountains. Finally, the Alps are considered a potential source of geothermal energy [55]. All these factors make the Japanese Alps a key land surface object for nature and the society of Japan.

2.2. Datasets Preprocessing

The geospatial data processing was performed in RStudio environment, Figure 2.

The workflow included the following general scheme: (1) loading packages, (2) obtaining data, (3) inspecting data.frame, (4) setting up the coordinate system, (5) calculating terrain characteristics (slope, aspect, hillshade) by the ‘raster’ package, (6) visualizing maps on a screen in RStudio, (7) plotting cartographic aesthetics in the ‘tmap’ package, and (8) mapping the prefectures of Japan in the ‘ggplot2’ package. The data for the terrain analysis include digital elevation models (DEM) from the Shuttle Radar Topography Mission (SRTM). The SRTM is widely used in research [56] due to its acceptable resolution (30 m × 30 m) and open-source availability: the National Aeronautics and Space Administration (NASA) Shuttle Radar Topography Mission (SRTM) provided DEMs for over 80% of the land surface on the Earth. The slope, aspect, hillshade and elevation maps were plotted using the RStudio environment [57] using the ‘raster’ package as implemented by [58]. Specifically, the following workflow was used for data pre-processing. The data were captured using the ‘getData()’ function of the ‘raster’ package of R from the available geospatial datasets of the University of California, Davis campus, CA, U.S.

The data were then reprojected to the Lambert conformal conic projection (LCC) by the following function: $crs\left(alt\right)<-"+proj=lcc+lat\_1=-30+lat\_2=40+lon\_0=140+datum=WGS84".$ The coordinate system is an important aspect of geographic data which is implemented in the R environment by a PROJ library [59]. It should be specially mentioned that the SRTM shows elevation data not as a bare-earth model but as a surface, which includes dense canopy forests and built-up areas in the estimation of the terrain. Although such nuances might be worrying in hydrological modelling, the SRTM is generally acceptable for country-level mapping, as in our case for mapping Japan using SRTM DEM.

2.3. Methods

2.3.1. 3D Modelling by GMT Scripts

Scripting by Generic Mapping Tools (GMT) followed the existing cartographic experience [60] with the use of a 3D package for the terrain modelling of Japan, as seen in Figure 3.

It should be noted that 3D modelling is an important approach to data processing, with many applications in engineering and natural sciences [61,62,63,64]. Most 3D modelling approaches in mapping address the problem of the visualization of the land surface in the forms of multiple views of the terrain by 3D data acquisition, e.g., light detection and ranging (LiDAR) [65,66]. The high cost of special hardware for obtaining LiDAR data and substantial manual processing reduce the operational flexibility of such approaches and limit their availability as practical applications for topographic mapping.

In this part, we present 3D terrain modelling based on the GMT open soft toolset. Such approach makes 3D mapping a more practical routine and enables many applications in Earth science with advanced cartographic data visualization of the terrain. Here, the 3D shapes of the terrain are presented in perspective projections for data interpolated by splines with heights obtained from the 3D coordinates of the raster grids of the ETOPO1 and ETOPO2 datasets. The differences in the local ruggedness of the terrain as represented by ETOPO1 and ETOPO2 are visualised and compared in 3D perspective plots with a rotated azimuth view. The algorithms included the ‘grdview’ module of GMT, as shown in scripts in the Appendix A.

2.3.2. Mapping the Prefectures of Japan

To illustrate the location of various prefectures of Japan, we used the "ggplot2" packages of R for plotting the regions of prefectures using arguments in the dataset for visualization. The data were captured by the ‘maps’ package of R [67] and ‘mapdata’, which provide the map databases, including the data on Japan: its prefectures, areas, etc. The advantage of this approach consists of the integrated use of ‘mapdata’ with the ‘ggplot2’ package of R. The ‘ggplot2’ package is a common graphical package designed for general purpose scientific visualization; however, it is applicable for cartographic purposes as well. Here, the main elements, i.e., the polygons of the prefectures in Japan, were plotted using the $‘geom\_polygo{n}^{\prime }$ function, which operates with data frames containing the coordinates of polygons and values associated with each of them:

geom_polygon(data=japan, aes(x=long, y=lat, fill=region, group=group)

The annotations of the axes were added using the ‘xlab’ and ‘ylab’ functions. The color scales were defined in the RColorBrewer and extended from the default, fixed number of colors to the number of prefectures in Japan, which was inspected by the ‘length’ function, as follows:

length(unique(japan$region))

Afterwards, the number of colors in the color palette was defined using the ‘colourCount’ function:

colourCount = length(unique(japan$region))

Following that, the color palette was expanded to the required number (47 prefectures of Japan) by the use of the ‘colorRampPalette’ function, as follows:

colorRampPalette(brewer.pal(name="Spectral", n = 8))(47).

Now, when the color palette was adjusted to the dataset, the map of prefectures of Japan was colored by assigning each individual color to each of the 47 prefectures, as follows:

scale_fill_manual(values = getPalette(colourCount)).

Other, additional cartographic elements included adding the titles, subtitles and captions by the ‘labs’ function, which enables the modification of the labels, annotations and captions. The rest of the aesthetics were added using $theme\left(\right)$ function, as shown in Figure 4, which provides an illustration of the script on the left and the resulting output map on the right part of the menu. These include, among others, plotting the legend, its orientation and annotations, defining the text size and font characteristics, selecting the color of background, and choosing the types of grid ticks and the frequency of breaks on a cartographic grid. The final map is shown in Figure 4.

The map of the prefectures of Japan (Figure 4) presented the visualization of the 47 prefectures made using the ‘ggplot2’ package of R, as well as the ‘mapdata’ and ‘maps’ packages used for data capture and processing. Each prefecture was colored as single-colored polygons modified from the default palette of R (‘Spectral’) using the palette extension by the ‘colorRampPalette’ function.

2.3.3. Mapping Morphological Features

Modeling slope, aspect, and hillshade is a technique for visualising terrain determined by DEM as a numerical data source. Here, the slope and the aspect were plotted first as major features. Following that, we visualised the hillshade, which is a derivative from the slope and an aspect of the elevation. The parameters of hillshade were adjusted and visually changed using various illuminating positions of the light source as degrees of simulated sun angle rotation. The four derivatives of the SRTM90 DEM grid were modelled and visualized using R as follows: slope, aspect, hillshade, and elevation.

The slope gradient, aspect and hillshade were defined by the algorithm of the ‘raster’ package of R at any point of the input raster grid (SRTM90) using local neighbourhood analysis. These variables are automatically defined by the machine using values of altitude (elevation) and its derivatives at or around each cell point on a raster grid representing the land surface. The ‘alt’ has a formal class ‘RasterLayer’ of package ‘raster’ with 12 slots.

Slope

The terrain characteristics of slope were calculated and modelled using the ‘raster’ package of R using the following sequence of commands. First, the color palette of R was created for visualization: $cols<-rainbow\left(255\right)$. Second, the slope and aspect were computed by following: $slope=terrain(alt,opt=‘‘slope")$. Third, the slope was visualised: $plot(slope,col=cols,main=‘Slop{e}^{\prime },xlab=‘‘lon",ylab=‘‘lat")$. Here, ‘alt’ represents the abbreviation of ‘altitude’.

Aspect

The aspect was modelled using the function ‘terrain’ with the option ‘aspect’ as follows: $aspect=terrain(alt,opt=‘‘aspect")$, followed by the definition of colors: $cols<-terrain.colors\left(255\right)$. Afterwards, using the same function $plot\left(\right)$ as for the slope modeling, the map of the aspect was plotted using the following snippet of code: $plot(aspect,col=cols,main=‘Exposur{e}^{\prime },xlab=‘‘lon",ylab=‘‘lat")$. The full script of R is provided in the Appendix A.5.

Hillshade

The topographic hillshade represents a remarkable phenomenon in cartography that occurs when shadow caused by the elevation cue three-dimensional shape perception. Various techniques exist to produce a monochrome 3D view of a terrain relief with the relative position of the artificial illumination representing the shadows from the natural sun. The effects from this shading are well reflected in a highly rugged terrain such as Japan, where high-elevation land surface uses strong luminance contrast, while low heights use low contrast. In attempts to visualise reliefs using the shading effects of a monochrome palette, efforts have been made since the onset of the GIS development and are considered in this work. The computation of hillshade was carried out using the following algorithm: $hill=hillShade(slope,aspect,angle=40,direction=270)$. The color palette was applied using the code snippet: $cols<-topo.colors\left(255\right)$. Plotting the map was carried out using the plot function as described above.

Cartographic processing

Mapping the morphological parameters in the ‘tmap’ package [68] included the processing of the computed raster layers by a specifically designed cartographic package where more elements could be added to the layout and more control over their aesthetics was available. The package used a function $tmap\_style$, where the general style of the layout was defined. Then, $tm\_shape$ was used to control the main data: the name of the raster, title, and subtitle. The $tm\_raster$ function was applied to process the color scale of the numeric variables (values of slopes in degrees), labels, and other cartographic details in the legend; see Figure 5.

The $tm\_scale\_bar$ and $tm\_compass$ functions show the parameters of the annotations of the auxiliary elements on the map. Thus, the $tm\_layout$ function controls the variety of the cartographic aesthetics necessary for proper visualization, such as the fonts and ticks or color of panel labels. The elements of the map affect the perception of the cartographic layout, which is crucial for overall visual appreciation by the reader. Therefore, map elements were adjusted using relevant functions of the package. The final map was visualized in an RStudio environment by calling the raster for inspection (here, ‘map1’) and saved using the $tmap\_save$ function: $tmap\_save(map1,‘‘Slope\_Japan.jpg",height=7)$. The same procedure was repeated with all the other three files (slope, aspect, and hillshade). The script visualized in RStudio is presented in Figures where the left parts demonstrate the script and the right part shows the output maps.

3. Results

The results of the geospatial modelling of Japan are structured into the 2D and 3D plots made using GMT, the regional mapping of prefectures performed by R, and morphological mapping (slope, aspect, and hillshade sections). Various R packages were used in the workflow, as described above and presented in the Appendices. The slope gradient, being a derivative of the altitude, presents a scale-dependent variable which changes according to the reduced or increased DEM resolution. Since the input data of this research are a 30*m SRTM DEM, the results of the slope modelling refer to the given resolution.

Marked physiographic variations in the Japanese Alps control the type and distribution of morphometric parameters in the following two regards. First, the heights change between the Hokkaido, Kyushu and Honshu Islands, with the highest elevation points in the central part of the Honshu (Mt. Fuji). Therefore, the variation of topographic ruggedness and slope steepness is primarily controlled by the geographic location, with the largest difference between the extreme highest and lowest points in the Central Alps. Second, the largest earthquake recorded in Japan, the Tohoku earthquake recorded ca. 371 km NE of Tokyo in 2011 (Miyagi prefecture, see the map in Figure 4), affected the topography of the country and increased slope instability and the risk of landslides.

The comparison of the exposure (Figure 5, left) and elevation maps (Figure 5, right) shows the trends in the North-South-West-East directions with regard to the topographic altitude of the land surface. The functional options related to the visualization by RStudio are discussed in previous sections. The presented maps are based on the ETOPO1 and ETOPO2 data with 1 and 2 arc second resolution, respectively, used for 3D modelling by GMT, as shown on the surface plot of the land relief of Japan with varied rotation. The 3D mesh model with isolines drawn on top of the surface is based on ETOPO2, and the grey-shaded topographic ‘waterfall’ plot based on the ETOPO2 grid with a view rotation of 115/30° as shown in Figure 6.

A 3D map with a view rotation of 65/30°is presented in Figure 7. The SRTM90 DEM, which shows the morphometric models, represents the slope, aspect, and hillshade relief in the land surface of Japan. The number of pixels (over 16,000 on a raster grid) was the greatest in the ‘gentle’ slope level (yellow color) compared to the others: 12,000 pixels for the ‘moderate’ slope (orange color), over 10,000 for the ‘strong’ slope (light red color), and 8000 for the ‘very strong’ slope (magenta color); three bins are covered by the class ‘extreme’ slope (purple color), and the rest (blue color) are represented by the less than 3000 pixels group, that is, the steepest mountain slopes.

The slope map shows the steepness of the mountain sides in the Japanese Alps in the Northern, Central and Southern Alps. The highs and hills in the raster grids are visualised in different sub-regions of Japan. The slope directions (Figure 8) revealed the following variations in the data: ‘gentle’, ‘moderate’, ‘strong’, ‘very strong’, ‘extreme’, and ‘steep’ slopes of the mountainous regions of Japan.

The aspect map (Figure 9) demonstrates that slope orientation, according to the compass direction (W-E-S-N), differs in various parts of the mountains chains of the Japanese Alps.

Thus, in the central part of the area (around 36° N), many slopes have a primarily west orientation (as shown by a red color), which is correlated to the geographic distribution of the Central Alps. However, in other parts, the southern (yellow color) and eastern (orange color) orientations demonstrate the prevailing values. Similarly, the northern orientation (green color) of the slopes contributes the least in the examined dataset.

Furthermore, the relationship of the slope aspect and elevation values indicating the altitude of the mountains might be the points of correlation. This should be emphasised since the impact of the geomorphic patterns on the morphometric characteristics helps detect the trends, showing the relief in the Japanese Alps. The hillshade map (Figure 10) shows a complex model based on the previously created maps of slope and aspect.

The levels of illumination were set at angle = 40°and azimuthal direction = 270°. The elevation map was used as a basis for the relief overlaid by the hill shading. The results of the regional mapping of prefectures performed by R are shown in Figure 11.

4. Discussion

The results of the SRTM90 DEM modelling identified major morphometric characteristics, such as slope, aspect and hillshade and demonstrated the uneven elevation in the topographic maps made using R scripting. Using a variety of R packages implies both the processing of data, exploration of data frames, manipulation of its structure and spatial data analysis (e.g., using the ‘raster’ package), and the cartographic aesthetic visualization of the layouts supported by specially designed packages such as ‘tmap’, ‘maps’ and ‘mapdata’.

The analysis of topographic attributes and DEM is of fundamental importance for reconstructing the genesis and development of landforms and, more generally, the geological setting of a specific area. The new ways of measuring, sensing, and analysing relief morphology carried out in Japan suggest that the range of the relief is reinforced by the high-elevation landforms, as well as external factors such as geomorphological processes and the climate, which contribute to eroding and modifying them. In fact, a strong correlation between geomorphology, topography and climate is commonly known, as also reported in Japan [69]. As proof of this, they pointed out that the climatic changes since the Late Glacial period have been responsible for the passage from the glacial to the permafrost environment in the current alpine zone, leading to modifications in the geomorphological processes and in the relief of the northern Japanese Alps.

Other important findings in geomorphic studies in Japan are well summarized by [70], who noted extensive sedimentation in mountain piedmonts and coastal fluvial plains and abundant sediment in steep watersheds. They furthermore pointed out the occurrence of hydro-geomorphological events in the areas of earthquakes and volcanic eruptions, as well as the post-glacial development of hillslope and flood processes along alluvial fans strongly controlled by the climate. In addition to the environmental factors, the relief directly influences the social-economic patterns through the possibility of road and building constructions and the potential triggering of landslides, which depend on the rock properties and slope steepness.

The geomorphic processes are largely driven by the gravity of the Earth and controlled by the slope steepness of the relief. As a result, the intensity of the surface processes, such as landslides, is affected by the displacement gradients connected to the geometric curvature and ruggedness of the relief. In such a way, practical applications of land surface maps in 2D and 3D representations enable the evaluation of the slope steepness quantitatively. In turn, the information retrieved from the calculated rates at which landslides are dismantling mountain slopes can be used for hazard and risk assessment for practical purposes of engineering geology, as well as for estimating potential slope instability [71]. Furthermore, other applications of land surface maps include the evaluation of environmental risks caused by climatic factors, such as precipitation intensity and the repeatability of rainfall or downpours [72,73] and visualising water surface dynamics in the context of flood and drought applications [74]. Finally, risk assessment and management in mountainous regions should consider other natural events, such as earthquakes or typhoons [75,76,77].

Due to the importance of topographic data visualization and the need for updated maps based on high-resolution data, there are many research reports on methods of topographic and bathymetric mapping and DEM applications for geospatial mapping. The latter includes, for instance, land surface classification based on DEM, reliefs, and geomorphological modelling using various GIS and geomorphometric computations of slope and aspect, as well as issues concerning the visualising methods of hillshade and DEM [78,79,80]. Since morphometric studies should always be supported by a spatial topographic representation, the questions of technical tools of mapping always remain actual for DEM-based studies. Our research contributed to this topic and showed that the integrated use of GMT and R in morphometric studies is an effective approach both for thematic mapping and for spatial analysis.

Graphics and maps, when created well, can provide eye-catching and detailed information on the land surface. This is why graphical approaches and methodologies of cartographic visualization are of high importance in Earth sciences. Despite a specific language syntax and approaches to data analysis that require mastering the tool, the GMT scripting toolset and R language are shown to be very promising tools for geospatial visualization, morphological modelling, and mapping, supporting the research in various aspects due to their different functionality.

Cartographic data visualization, regardless of what GIS software is used, is an integral part of the complex workflow of geospatial research. To mention some steps in a simplified process of geospatial data analysis using cartographic tasks, this includes data capture, preprocessing, projecting, analysis of content of the datasets, modelling variables, visualising and plotting the maps with controlled layout, and others. Applying the methods of scripts extended to the machine learning approaches to a cartographic workflow facilitates the process of spatial analysis. Specifically, it helps to increase the accuracy and precision of the final results and maps and significantly decreases the time of data processing due to the automation of geospatial data processing. Moreover, we have shown that the use of scripts and programming methods are more efficient and therefore more suitable for cartographic workflow due to repeatability. Our methods could be easily extended to also use other datasets and geographic extent of data. In fact, the computational procedure adopted here to investigate the setting of reliefs and mapping land features in 2D or 3D models could be applied to any target scales, be it global downscaled modelling or upscaling to regional (prefecture-level) and local (city-level) terrain models.

5. Conclusions

In this paper, we have designed adaptive script-based algorithms for plotting morphologic maps using GMT and R that are able to visualise maps accurately and effectively. Our two algorithms provide a tradeoff between the computational approach of R for geospatial data processing and the cartographic performance of GMT. For future relevant works, we suggest that one chooses the appropriate algorithm based on their mapping goals and available dataset. Both programs have free access as open-source tools and extensive functionality for geographic data analysis and visualization. The presented methods can be used and adopted to other regions and areas with changed relevant coordinates and modified attributes for 3D modelling (e.g., elevation range). We have shown that the use of both these tools for mapping is a very effective approach for spatial analysis that can be used instead of or besides GIS, as a complimentary script-based method.

We furthermore identified the existing challenges such as the need to use, install, and load a variety of packages in R to fit in various tasks, such as: geomorphometry, classic mapping, spatial analysis, data capture, and data conversion. In contrast, GMT enables us to perform all the steps of a cartographic workflow in the same session by calling necessary modules. On the one hand, using scripts can be an easier approach compared to GIS, but on the other hand, it can be a challenge, since it requires finding a suitable package, installing and activating it, and using it by applying its specific functions and syntax using GMT or R.

The integration of our method based on the GMT and R algorithms for processing various geospatial data is straightforward since many data formats can smoothly be imported and processed both by GMT and R. Furthermore, R has certain embedded datasets that we used for plotting the prefectures of Japan and for visualising physical features of relief. Working on grey-level images for 3D modelling by GMT enabled us to detail the land surface features and demonstrated the high accuracy of the ETOPO grids. With respect to our work, which exploited script-based mapping, we have designed and presented a series of thematic maps on Japan that can be applied for other regions and countries.

As we have shown in this study, the use of GMT and R for 2D and 3D mapping in geographic analysis is highly effective and recommended, especially for perspective visualization of the terrain. However, we should also notice that this method is not yet as popular as traditional GIS. This can be explained by the non-trivial approach of using scripts and programming techniques in mapping and cartographic data processing, as well as certain skills required for coding and mastering the syntax of both the tools. Nevertheless, the improved workflow of mapping and results suggests that both GMT and R are effective for geosciences and recommended for mapping purposes.