. Introduction
Current paper introduces the use of Generic Mapping Tools (GMT) for the cartographic workflow aimed at the geological data visualization. Among the wide range of GIS software used for mapping and geospatial data modelling, GMT stands apart from the traditional tools such as ArcGIS (ESRI Team, 2010), QGIS (QGIS, 2019) and MapINFO. The particularity of the GMT comparing to the standard GIS consists in its fundamentally different approach towards geodata processing. The GMT, originally developed by Wessel and Smith (1998) is not a classic GIS but a toolset of the shell/bash scripts that should be written by a cartographer for mapping rather than operating with standard GUI and panel tools. Each line of the code in GMT script presents a command used for specific purposes: visualizing specific elements, modelling grid, generating colour palette, drawing lines on the map, adding logo, scale bar, setting up parameters (fonts, offset, map layout), adding directional rose and so on. The sequence of the several GMT commands together combine a script that can be run from the command line of the GMT console. Running a GMT script results in an output map. The resulted map is generated as a PostScript image file (ps) that can be further converted into standard graphical format (e.g., tiff, jpeg).
Marine geological data analysis include two major approaches: GIS based data modelling and statistical analysis. The first approach, the statistical analysis of the marine geological observations can be found in various works (Brown et al., 2012; Gallo et al., 2015; Harris and Whiteway, 2011; Lemenkova, 2019c,e; Meibom and Anderson, 2003; Ravisankar et al., 2015; Strick et al., 2018), in which the authors test various statistical approaches, software and libraries by introducing the latest algorithm technologies specifically for oceanographic data analysis. The second approach, a GIS mapping for the geospatial data mapping, consists of the multiple reports with mostly usage of ArcGIS. A brief summary of the up to date applications of GIS in the marine studies can be read in the following works: Ichino et al. (2015); Metz et al. (2014); Rovere et al. (2019); Sen et al. (2016); Yesson et al. (2011). Examples of the photogrammetric approaches of image analysis include, for example, Li-DAR and processing digital elevation model (DEM) (Altuntas, 2019), GNSS and GPS data geodetic modelling (Banasik and Bujakowski, 2017; Tomaszewski et al., 2017). The traditional GIS, for example, ArcGIS, is used in the thematic mapping: environmental analysis (Klaučo et al., 2013, 2017); geomorphological assessment (Lemenkova et al., 2012); marine cartography (Suetova et al., 2005). Advanced data analysis with programming and scripting based approaches may be used for the urban mapping, for example, CityGML XML with 3D modelling (Janečka, 2019), marine geological mapping, for example, R (Lemenkova, 2019d) or Python (Lemenkova, 2019c), combination of GMT and ArcGIS (Gauger et al., 2007).
Current paper presents the advantages of GMT to increase its popularity. The fundamental difference of the GMT from traditional GIS lies in its scripting approach which enables semi-automated processing of the data in large volumes (e.g., big data in geosciences). In total, the open source GMT collection consists of approximately 80 command-line modules for the processing of geographic and Cartesian data sets. Several grid and vector layers are embedded in the GMT, such as the world coastline (Wessel and Smith, 1996) that can be visualized directly using coordinates and choosing projection. A variety of cartographic projections and transformations (over 30) are presented in the GMT enabling cartographer to select coordinates, suitable map orientation and additional elements (e.g., standard parallels, meridian). The effectiveness of the GMT consists in flexibility and variety of the modules that enables to perform high-quality mapping. The flexibility of the usage of the GMT modules for data modelling proves it to be an excellent tool set for the geomatics and cartographic data visualization.
The GMT compatibility with Unix commands enables a cartographer to insert echo, rm Unix commands directly in the GMT code and to integrate GMT data modelling with further data processing to Matlab. Unix echo is a command available in shell Unix scripts. Its meaning is as follows. It puts the string texts passed as arguments to the visualization on the computer screen or printed output. In other words, it enables to put the text to the screen or a computer file as part of a sequential pipeline of Unix commands. For example, using echo, one can write directly a piece of text on any location of the map, as illustrated in Fig. 14. The text on graphs A and B in Fig. 14 in red, blue and brown colours (’Pacific Ocean, Pacific Plate, ’Median stacked profiles’, . . . ) were written using the echo Unix command. The same can be seen in Fig. 15, subplots B–E and G–J, where formulae (placed in lower left corner on the light yellow background) were also visualized using echo command. The placing of echo command can be seen in the fragment in Listing 10 (GMT code for trend modelling). The rm Unix command is an abbreviation meaning remove files. It removes the auxiliary files after the execution of the script to avoid the overstock of the computer memory. It is commonly used in the programming environments aimed to delete unnecessary files automatically, after the batch processing is finished. Besides, GMT enables both 2D and 3D data modelling and high-quality data visualization, which is explained in the following sub-chapters stepwise with code snippets.
Figure 2
Topographic map of the study area: Sea of Okhotsk, Kamchatka Peninsula, Greater Kuril Chain and Kuril– Kamchatka Trench. Numerical data source: ETOPO 1 Global Relief Model 1 minute raster grid.
Figure 4
Composite overlay of the 3D-topographical mesh model on top of the 2D grid contour plot. Contour bathymetric map source: ETOPO 5 min grid resolution. Azimuth rotation: 165/30.
Figure 5
Composite overlay of the 3D-topographical mesh model on top of the 2D grid contour plot. Contour bathymetric map source: ETOPO 5 min grid resolution. Azimuth rotation: 135/30.
. General Situation
Geographic Location
The study are is focused on the Kuril–Kamchatka Trench, Kamchatka Peninsula and the Sea of Okhotsk located in the Far East, North Pacific Ocean (Figure 1). The Sea of Okhotsk is a semi-enclosed sea, connected in terms of bathymetry with the open Pacific Ocean through the Bussol and Krusenstern straits (Maiorova and Adrianov, 2018). The Sea of Okhotsk is not isolated, has maximal depths of 3374 m. It has deep-sea straits to the NW Pacific Ocean, the deepest of which are the Bussol with depths reaching 2300 m (Avdeiko et al., 2007; Brandt et al., 2018) (Figure 2). The Sea of Okhotsk is connected with the Pacific Ocean hydrologically 14through the Bussol Strait. The upper water masses of the Kuril–Kamchatka Trench region are strongly influenced by the currents in opposite directions: the Oyashio coming from the Arctic southwards into the Pacific Ocean and the Kuroshio coming from the South China and flowing in the northeast direction (Belkin et al., 2009; Tyler, 2002).
Geology
The region of the Far East region is notable for particularly high seismicity, repeated earthquakes and unstable geologic settings. The particular seismic features of the area are submarine earthquakes generated within the island-arc slope of the Kuril-–Kamchatka Trench (Figure 3) seismic belt that can reach M = 9.0 and trigger tsunami waves. Comparing to other natural hazards, such as earthquakes, floods, and typhoons, tsunami is ranked fourth in terms of damage to human life and infrastructure losses.
Among the most recent such events, there is 2006–2007 Great Kuril Earthquake Sequence that involved coupled under thrusting and extensional faulting on a large scale in this area (Lay et al., 2009). Most of the earthquakes within the Kuril– Kamchatka Trench are located in an area between the trench axis and the edge of the Kuril islands and land shelf, where the Pacific plate is subducting beneath the island-arc zone of Eurasian continental lithosphere (Lobkovsky and Sorokhtin, 1979).
Tectonics
During the Late Eocene, the Kuril–Kamchatka subduction zone and volcanic front was located in the Sredinny Mountains of the Kamchatka Peninsula, and migrated ca. 150 km eastwards thereafter, towards the eastern shore of the Kamchatka Peninsula. This migration resulted in the formation of the Sredinny Mountains on the Kamchatka Peninsula, a series of volcanoes on the Central Kamchatka Depression as well as Eastern Volcanic Plateau (Barr and Spagnolo, 2013).
The seismic situation is reflected by the bedrock in the Kuril– Kamchatka Trench, mostly dominated by the Quaternary and Miocene–Plocene volcanic complexes (Persits et al., 1997). As a consequence, the tsunamis in the Far East region in general and in Kuril–Kamchatka area in particular, are frequent, extensive, and devastating compared to the other regions of the Pacific Ocean (Hatori, 1971; Soloviev, 1968, 1972).
Marine Biology
The uniqueness of the Sea of Okhotsk consists in its richness of the marine deep-sea biota. In particular, the high level of biodiversity is shown in the eutrophic area of the Kuril– Kamchatka Trench, which is explained by the turbulence of the water masses. The abyssal echiurans fauna in the Kuril– Kamchatka Trench and the deep-sea basin of the Sea of Japan show an intermediate level of species biodiversity and species richness (Maiorova and Adrianov, 2018). However, in other geographic places of the Kuril– Kamchatka basin, a very rich marine fauna is recorded (Fischer et al., 2015; Zenkevich, 1963). In general, the geographic location of the trench and abyssal plain affects its biodiversity (Schmidt et al., 2019); despite a vertical separation of the seafloor of the trench from the euphotic zone, it underlies a richly productive boreal region, because it is located near the coasts of the Kamchatka Peninsula and Greater Kuril Chain. Steep slopes of the trench serve as depositional centres of the organic carbon accreted on the seafloor via the lateral transport from the continental margins forming deepsea biological hotspots (Danovaro et al., 2003; Itoh et al., 2011).
. Methodology
Methodology includes application of the GMT scripting toolset for the raster thematic maps visualization and automated digitizing of the profiles across the Kuril–Kamchatka Trench crossing the trench in a perpendicular direction. A sequence of the GMT codes is explained below in the relevant subsections. The GMT coding was used to visualize raster and vector data, perform geomorphological modelling, data analysis and comparison of the two segments of the Kuril–Kamchatka Trench, as described below. Several cartographic projections were tested for maps visualization in GMT: the Oblique Mercator projection (shown in Figure 1) was used to show Kuril–Kamchatka Trench area along a great circle other than the Equator or a Transverse Mercator. Figures 7 to 14 are presented in Mercator projection. Two maps are plotted in the conic projection: equal-area Albers projection (Figure 3) and equidistant conic projection (Figure 6). Figure 4 and 5 are presenting 3D modelling visualizing geomorphology of the study area.
Research questions and purpose
The main purpose of the research presented here is to provide a technique for the spatial modelling and mapping of the oceanic trench by remote sensing methods: visualizing raster data grids using GMT without having a fieldwork in the study area. To achieve this, the role of the geological and tectonic setting that may affect the geomorphology of the trench were investigated and reviewed in the Section 2 (General Situation).
The cartographic means and GMT modules taken to investigate Kuril– Kamchatka geospatial setting is described below. Therefore, the objective of this research was to derive new GMT based techniques applied for cartographic mapping, visualizing and geomorphic modelling of the geospatial data with a case study of the Kuril– Kamchatka Trench region, Far East, North Pacific Ocean.
The general questions of this research are: first, how different the topographic bathymetric settings are in northern and southern parts of the trench, and how to perform modelling technically by means of GMT scripting using the available data sets. With this aim, Figure 14 and Figure 15 visualize the variations in the gradient slope of the trench. Second question was to understand, if the topography varies across the length of the trench. With this aim, the topographic contours were mapped to visualize the whole study area of the Sea of Okhotsk and North-West Pacific Ocean (Figure 3) and modelled through the raster grid specifically for the Kuril–Kamchatka Trench (Figure 10). The geophysical settings of the study area are visualized through modelled geoid (Figure 6 and Figure 7) and gravimetry (Figure 8 and Figure 9). The enlarged fragment of the gravity along the Kuril–Kamchatka Trench is presented in Figure 11.
Data
In this research, high-resolution geodata collection was performed from the website USGS and available embedded layers of the GMT from the SOEST, University of Hawaii, and NOAA Laboratory for Satellite Altimetry (Wessel and Smith, 1996; Smith and Sandwell, 1997b):
The main data include ETOPO1 1 arc-minute global relief model of Earth’s surface from NOAA, used for visualization of the bathymetry taken from the official website: ETOPO1 the original grid raster file (earth_relief_01m.grd) in NetCDF format was produced for modelling and mapping in GMT, that is, Figure 2, Figure 3, Figure 14 (modelled profiles are based on the ETOPO1), Figure 15 as a visualization of the modelled profiles based on ETOPO1.
The ETOPO5 5 arc-minute global relief model of the Earth’s surface (earth_relief_05m.grd) was used for modelling Figure 3 (contour map), Figure 4, Figure 5 (3D maps).
Gridding bathymetric contours from the xyz data were based on the data taken from the public official website: Global seafloor topography from satellite altimetry. From there, the necessary square was selected by the coordinates following West-East-South-North usual system (in this case, the square was the following: 144◦W–162◦W; 40◦N–51◦N). The derived xyz data were used as a table and then processed and visualized on Figure 12, Figure 13.
The raster data on gravity and geoid were used from the Scripps Institution of Oceanography and visualized on Figure 7, Figure 8, Figure 9.
Coastal borders on all the maps were based on the NOAA: Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG).
3D-topographical mesh modelling
The 3D-modelling of the trench highlighting the bathymetric patterns of the study area is presented in Figure 4 and Figure 5 plotted using the main module grdview, which reads a 2D raster grid on Kuril–Kamchatka area and produces a 3D perspective plot.
Numerical data source for the grid: Smith and Sandwell (1997a). Visualizing ETOPO data set with 5 minute resolution was performed through the sequence of main GMT modules: img2grd and grdimage. Because processing of the 1-min ETOPO presented a too detailed and large 3D map, it was questionable if applying higher resolution data continuously over the whole region of the Kamchatka area is beneficial, since the gain in information was too small for the price. Therefore, the 5-min resolution for the 3D modelling was chosen for the selected segment presented of Figure 4 and Figure 5 (rotated). The makecpt module was used for generating colour palette. Further techniques include drawing a mesh painted by a coloured surface. In this case a -Qsm command in the presented code below stands for added mesh lines plotted as contours on top of the surface, and added a plane at the z-level equal to -7500 (that is a bathymetric depth in this case). It is worth noting that two maps are combined together on one layout, presented with a set up distance between them: a 2D oblique contour and a 3D map. The chosen colour palette is artificial ’rainbow’. As can be seen from the code below (Listing 1), all commands generate a part of the map that is read into the ’ps’ file defined initially. At the next steps, the elements of the 2D map are added stepwise. The 3D grid was added on the map through the following code (Listing 1). Here the -J -R stand for the projection derived from the basic 2D map, R stands for the region, -p165/30 shows the azimuth angle of the 3D grid, -Qsm -N-7500 shows the baseline of the level of the 3D. Figure 4 is showing that the mesh 3D grid is the rotation azimuth angle 165◦ versus 135◦ in the Figure 5.
Geoid image modelling
The geoid data were taken from the Satellite Geodesy research group at Scripps Institution of Oceanography, University of California San Diego, Global Geoid Online Data Set. The geodetic situation in the study area was visualized by the GMT modules (Figure 6) using the following GMT code snippet (Listing 2): stepwise. The output image showing modelled geoid is presented in Figure 6. The scale colours used for the data visualization are artificial ’rainbow’ to highlight the differences in the elevation of the shape that the ocean surface would take under the influence of the gravity and rotation of the Earth. The geodetic situation in the study area with shaded overlapped relief of the adjacent study area was visualized by the sequence GMT modules (Figure 7) using the following code (Listing 3).
Gravitational regional modelling
Gravity model reflects the deformation of the Earth’s surface and is connected with the tectonic settings of the study area, such buoyancy stress of the tectonic slabs, vertical gravity signal in the forearc (Boutelier and Oncken, 2011; Ramberg, 1967). The key GMT modules grdimage and grdcontour were used to generate gravity image using available numerical data (Sandwell et al., 2014) with shading using geoid data. Marine free-air gravity modelling is a useful approach for computing non-isostatic bathymetry. In this case, the additional effects of the sedimentation on the seafloor, and lithospheric subsidence related to age are removed from the bathymetry. As a result, a non-isostatic topographic map presents an improved approximation to the deformed shape of the Earth’s seafloor and tectonic plate (Zhang et al., 2014). The selected colour palette is ’GMT_haxby’ which is a GMT embedded Bill Haxby’s colour scheme specially designed for geoid and gravity. The data on the marine free-air gravity anomalies (FAGA) were derived and estimated from the available satellite altimetry missions from the CryoSat-2 and Jason-1 radar altimeter data. Modelling the marine free-air gravity anomaly along the Kuril– Kamchatka Trench (Figure 8) was performed using GMT modules by the following code snippet (Listing 4) using available data set (Smith and Sandwell, 1997b). Figure 9 shows a marine free-air gravity modelling along the study area of the Kuril– Kamchatka Trench on the junction between the Pacific and Okhotsk tectonic plates. Finally, the base map including grid, title, coastline, direction rose and scales were added via the GMT pscoast module.
Surface modelling
In addition to the main study area including the trench area (140◦E/170◦E/40◦N/60◦N), the area of the trench was miniaturized for surface modelling compared to that in the major map series with the enlarged scale of the square with coordinates 144◦E/162◦E/40◦N/51◦N, as can be noticed on the respective maps (Figure 10 and Figure 11). Surface modelling by the curvature splines from the ASCII data was performed and analysed with a given data of topography and gravity. The geoid data were taken from the Satellite Geodesy research group at Scripps Institution of Oceanography, University of California San Diego, Global Geoid Online Data Set: http://topex.ucsd.edu/cgi-bin/get_data.cgi Two datasets, on topography and gravity, were modelled using these data, as described below in the respective sub-sections.
Topographic modelling
To perform surface modelling of the topography in the study area by GMT surface module (Figure 10), the continuous curvature splines of 0.25 was set up to derive explicitly high resolution data in the output grid. The resolution of the initial ASCII XYZ (Sandwell and Smith, 2009) input raw data set was 1 minute. This topographic surface model has its bottom slightly below -10000 m, and maximal values not significantly higher than 2000 km on the Kamchatka Peninsula. The GMT modules blockmean, surface, grdgradient, grdimage and psbasemap were used based on the relevant methodology (Wessel and Smith, 2018) to model the topographic surface along the trench (Figure 10). The colour palette for topographic surface modelling was selected as ’GMT sealand’, which is a Smith bathymetry/topography with scale R = –6000/ + 3000 (m), H = 0, C = HSV. As can be noticed from the above code (Listing 5), the block-mean filter was used before the data processing to achieve data smoothing and interpolation. At the output data, the azimuth shading was selected to 45◦, that is north-east illumination, as the best contrasting the geomorphic features in the study area.
Gravity modelling
Gravity surface modelling was done with the aim to visualize the geophysical situation and settings in the study area. Figure 11 shows the surface gravity model along the Kuril– Kamchatka Trench computed from the ASCII XYZ table data (ASCII XYZ table data source: (Sandwell and Smith, 2009)). Modelling numerical ASCII XYZ table data was done using surface GMT module designed for raster data processing, and auxiliary modules: blockmean for interpolation and smoothing, grdgradient modules was used for making gradient illumination with azimuth 45◦, grdimage was used to visualize the resulting image. The GMT modules psbasemap, pscoast, pstext and psscale were used to add cartographic elements: grids, coastline, colour legend, directional rose, scale bar and text annotations. The colour palette was selected as ’GMT relief’ which is a Wessel/Martinez colours for topography with elevation range (m) R = –8000/ + 8000, H = 0, C = RGB. The block-mode filter was applied to the initial raw standard input data set to arbitrarily read located x, y, z triples. In this way, the raw data matrix was initially standardized and normalized for the decimating and averaging x, y, z data as a pre-process before running surface. It also ensured that each determined variable now had equal weighting in the statistical analysis. The output data were written as estimates of the position x, y and value z for every non-empty block in a grid region of the selected trench area to the standard grid output. The most crucial GMT code snippets of code used to generate surface by surface module (for data modelling) and grdimage module (for visualization) are presented in Listing 6.
Comparison of the 2 modelling approaches: nearneighbor and XYZ2grd module
Modelling surface in GMT can be performed using several various approaches nearneighbor, surface, greenspline, triangulate, as well as XYZ2grd that reformats existing data to a grid structure. In this research two GMT modules were tested: XYZ2grd and nearneighbor (Figure 13).
The difference between the two approaches consists in the algorithms used in the modelling and the relevant commands. The topographic data were downloaded from the available dataset: http://topex.ucsd.edu/cgi-bin/get_data.cgi In this particular case, the coordinates of the study area square are: -144◦–162◦E; 40◦–51◦N.
Nearest Neighbour
The second modelling approach was done using nearneighbor GMT module (Figure 13b) using methodology of Nearest Neighbour algorithm for grid modelling. The algorithm of the Nearest Neighbour was used to model the surface from the xyz data using the commonly known approach that consists in the finding the point in a data set that is the most similar in its attributes to a neighbour point. The closeness is typically expressed in terms of a dissimilarity function: the less similar the points with coordinates, the larger the function values. The GMT code used for the modelling is presented in Listing 7.
In this case, the module reads (x, y, z) triples coordinates and depths values from the standard input table and uses embedded into program mathematical algorithm ’nearest neighbour’ to assign an average value to each node that have point within a radius centred on the node. The average value is computed as a weighted mean of the nearest point from each sector inside the search radius. The weighting function used is
where
XYZ2grd
The second XYZ2grd approach, that is XYZ2grd module reads one or more xyz tables and creates a binary grid file from the raw file. Modelling using XYZ2grd (Figure 13a) is done using script explained stepwise in Listing 8 where topo_KKT.xyz is the name of the processed file, xyz2grdKKT.nc is the name of the output file defined by the command -G, -R command shows the lon/lat region in min/max coordinates (144◦/162/◦/40◦/51◦ in this case), -I5m stands for the grid spacing and -V for verbosity level (that is, computer’s explanations of possible errors in the code). The explanations of the table examination is as follows: N is the total number of observations, the coordinates for the minimal depth (-9677 m) are 144.0083◦E 39.9976◦N, the coordinates for the maximal depth (2143 m) are 162.0083◦E 50.9968◦N.
After that the contour lines, cartographic elements (grids, annotation of Mercator projection, information on the computing date, GMT logo) were used by auxiliary GMT modules: grdcontour, pscoast and logo.
Automatic digitizing of the cross-section profiles using grdtrack GMT module
Plotting cross-sectional profiles, or transects along the selected segment of the study area is a common technique in geological and geographical investigations where the area is studied using a sequence of the cross-sections and then interpolated for the region with certain approximation level. The techniques of the semi-automated digitizing of the bathymetric data is presented with a case study of AutoTrace vectorizer by (Schenke and Lemenkova, 2008) where the authors illustrate vectorization of the paper maps (e.g., in .bmp or .tiff format) and converting them into a digital format with isolines as vector objects. Examples of the geological cross-section profiles plotting and data analysis can be seen in Chen and King (1998); Hayes (1966). More recent work include: Cortés-Aranda (2018); Feng et al. (2018); Lee et al. (2018); Lemenkova (2018a); Montaggioni et al. (2018); Vázquez et al. (2015).
A series of the topographic cross-section profiles was plotted along the trench using the sequence of the GMT modules grdtrack, trend1d, grdtrack: 26 profiles on the south and 31 on the north. Their characteristics are: 400 km long, spaced 20 km, sampled every 2 km.
Technically, a data set of the transects was digitized automatically using grdtrack module and a sequence of the auxiliary modules. The full details of the GMT bash code sequences and characterization of the module options and elements including 5 UNIX progs (echo, rm, cat) are provided as follows (Listing 9). The methodology is explained for the northern segment, and repeated correspondingly for the southern with changed coordinates for the start and end points.
First, the two points (starting and ending) were selected along the Kuril– Kamchatka Trench using UNIX cat utility (153.5◦ 45.5◦ and 158.5◦ 50.0◦ stand for the lon/lat coordinates for start and end points, correspondingly). Then the cross-track profiles were generated with set up parameters. They were then stacked using the statistical mean. Now as the northern segment was digitized, the procedure was repeated for the southern part with the coordinates 148.0◦E 43.0◦N and 153.5◦E 45.5◦N. The orthogonal cross-section bathymetric profiles were then drawn automatically using the presented script (Listing 9) across the Kuril–Kamchatka Trench with a length of 400 km (Figure 14).
In this study, the cross-section orthogonal profiles were taken through the Kuril–Kamchatka Trench to show the variety of the steepness angles. The cross-track profiles are 400 km long, spaced 20 km, sampled every 2 km to achieve evenness of the sampling (Figure 14). Two segments are digitized and compared: northern and southern parts separated by the Bussol Straight connecting the Sea of Okhotsk with the open Pacific Ocean. The depth values were calculated by each profile, and their mean were demonstrated as red lines (Figure 14A and Figure 14B) and median values were modelled thereafter. The output model is provided at Figure 14C showing the initial bathymetric map with applied artificial colours (’rainbow cpt’) and profiles across the Kuril–Kamchatka Trench. Automatically digitized two series of the profiles are shown in Figure 14A and Figure 14B, respectively. As compared to the traditional digitizing methods available in ArcGIS and other GIS software, automated digitizing process enabled by the GMT increased precision and speed of the plotting, as well as highlighted submarine geomorphic features of the trench, such as oceanward slopes, terraces and compare slope steepness on the continental side.
Trend modelling using trend1d GMT module
Final part of the research includes modelling gradient curves of the Kuril–Kamchatka Trench in northern and southern segments using trend1d module of GMT. Using various mathematical algorithms and statistical methods for data analysis in geoinformatics is presented in multiple works, for example, regression analysis (Bazan-Krzywoszanska and Bereta, 2018), cluster analysis (Lemenkova, 2019a), hierarchical dendrograms (Lemenkova, 2019b) and Pearson correlation analysis Hlotov et al. (2019); Lemenkova (2018b). Combination of the mathematical approaches with geospatial modelling increases the precision and reliability of data modelling by visualizing general trends in behaviour of the studied phenomena, the direction of the curves when on the graphs enables to highlight the relationship between the variables to identify possible dependencies. Since directed fieldwork is not always possible due to the remote location of the study object (as in case of the hadal trenches) or is expensive, the modelling of the available data sets is the optimal solution of geodata processing. In view of this, using mathematical approximations, advanced tools of data analysis for the visualization is indisputably important.
This research step used tables generated at the previous research part (automatic digitizing of the cross-section profiles using grdtrack GMT module). Now the tables contained the x-y coordinates and z values for the bathymetric data set. The tables were processed to visualize the overall trend of the slope curvature of the two segments of the trench. This was assessed using statistical data processing embedded in GMT. Modelling trend regression for the cross-section profiles (Figure 15), computing trends for the profiles by four types of mathematical models:
m2t = a + b · t,
m3t = a + b · t + c · t2,
m5t = a + b · t + c · t2 + d · cos 2π · t + e · sin 2π · t,
(t) = y(t) – m5t (plotting residuals).
The resulting graph of the trend regression models for the cross-section profiles of Kuril–Kamchatka Trench was plotted using the GMT codes sequence. The principle of this work is included in the following Listing 10.
. Results
The processed data included using ETOPO data (1 and 5 min resolution grids cut via the grdcut module), CryoSat and Jason-1, together with GMT-embedded data sets available at modules pscoast, psbasemap (Figure 1 and Figure 2). Main geomorphic elements that compose the landscapes along the Kuril Kamchatka Trench and the Sea of Okhotsk were visualized and mapped. The particular distribution of the spatial variations over the study area in correlation with the geophysical data is presented on the marine gravity anomaly field and modelled geoid of the study area. A resulting series of the thematic bathymetric and geological maps covering seafloor and coastal areas consists of easily interpretable maps presented both in natural topographic colours and artificial colours visualizing geomorphology of the trench, the geological and geophysical settings.
Along the south-eastwards direction of the study area (Figure 8), the values of the marine gravity anomaly decrease with the biggest relative marine gravity anomaly is 450.0 mGal; along the set of profiles southwards, the density decreases reaching -340.0 mGal, and the biggest relative pressure of gravity anomaly is 300.0 mGal (Figure 9). Thus, the absolute values of the highest relative pressure gravity anomaly (450.0 mGal) are detected on the north-west part off from the trench. The detected gravity values range from -220 to +560 mGal relative to the model of the study area. A visualization of the grid data set of the Kuril–Kamchatka Trench (Figure 9) derived from CryoSat-2 and Jason-1 (Sandwell et al., 2013) are compared with the results from the numeric modelling from the ASCII XYZ format. The findings from the surface gravity modelling section show a systematic decrease in gravimetric values notable in the centre of the trench, highlighted by dark navy colour (Figure 11).
With regards to surface modelling from the grid data sets, two different methods are used to model vertical marine gravity gradient: using grdimage GMT module (Figure 8 and Figure 9), and generating surface from the ASCII XYZ table data with results compared in the output visualization (Figure 10 and Figure 11). Coloured geoid image modelled for the study area shown high values of over 80 mGal in the southern part of the trench (Figure 7) while lower values are detected in the south-eastern part of the study area (168◦E 42◦N). Likewise, different strategies were investigated to estimate the geoid visualization from the observed raster grid data and ASCII data sets.
Automatic digitizing cross-section profiles topographic modelling (Figure 14) revealed that southern part of the trench’s geomorphology is characterized by more steep slope gradients while northern part has more gentle and soft shape form of the slopes. The results of the comparative analysis of the two distinct parts of the trench located north- and southwards from the Bussol Strait show that southern part is deeper reaching -8,200 m depth while northern part has -7,800 m at maximal records. At a southern segment, the seafloor depths increase with increasing latitude gradually while moving southwards. The profiles located in the northern segment are shallower, which illustrates tectonic and geological local variations, as well as different sedimentation of the northern and southern part of the trench. For example, the trench-parallel, linear system of the islands in the southern part of the Greater Kuril Chain provides a natural barrier to the extensive downslope sediment transport from the Okhotsk Sea to the trench. The submarine terraces are presented along the southern part of the trench axis as compared to the more straight northern part (Figure 15). As a consequence, such differences affect the landforms of the profiles in the respective segments. The seafloor geomorphology on the north is characterized by the gentle slope patterns and depths not exceeding 7,400 m. Northern part of the trench is the island-ward sloping surface of the Kuril islands. Here, the profiles are mainly flat and take a parallel form at the shallower depths that extend to about 7,200 m in the seafloor. The northern profiles show shallow depths and gentle slope shape following the enlargement of the trench valley. Stacks of local small hills and bulges across the trench indicate sediment deposition along the trench slope in its northern part.
. Conclusions
The research novelty lies in the comparative thematic mapping of the Kuril– Kamchatka Trench based on the several raster grids visualization, and automatically digitized profiles for geomorphological modelling of the trench. The developed methodology of the GMT-based mapping applied for the geomorphic modelling tested, presented and explained functionality of the several GMT modules. Automated digitizing of the orthogonal profiles crossing trenches in the perpendicular direction showed variations in the bathymetry in its northern and southern parts. The shape of the landforms and steepness gradient of the Kuril– Kamchatka Trench were visualized and compared. The traditional manual vectorizing is usually a laborious work prone to mistakes. Using GMT for the automatization of the digitizing, as demonstrated in this research, improved the quality, speed and precision of the data processing.
The vertical gravity gradients calculated from the raster grids show a decrease in the gravimetric values notable in the centre of the trench. The results on the geoid models showed high values of over 80 mGal in the southern part of the trench while lower values in its south-eastern part (168◦E 42◦N). The complex relationship between the geodetic and geophysical properties of the Kuril– Kamchatka Trench involves diverse factors: geoid height, vertical deflection, vertical gravity gradient (that is, the curvature of the ocean surface), gravity gradient, ocean dynamic topography and gravity anomaly. Taken together, they affect the geomorphic structure of the trench. Mapping gravity anomaly contributes to the precise bathymetric studies, since the correlation between the gravity and depth may facilitate to more precisely identify some submarine landforms: small seamounts, linear ridges, local deformations on the seafloor and so on. Global gravity grids reveal submarine volcanoes on the seafloor greater than 1000 m (Smith and Sandwell, 1997b). Since bathymetry of the seafloor reflects the complex movements of the tectonic plates, when they move apart from the spreading centre. Therefore, a complex investigation of the study area is important for the thorough understanding of the topographic patterns. Presented approach contributes to the development of the methodology of the GMT based marine mapping.
Comparison of the two surface modelling approaches for the generation vector contour map from the ASCII XYZ data, XYZ2grd and nearneighbor, shows that XYZ2grd presented more precise bathymetry lines, better detailed contouring. That is, even small geomorphic objects, such as small islands, are shown with higher level of smoothness. On the contrary, the nearneighbor approach showed generalized lines, less smooth and more straight, and rigid contours. Based on this, the results derived from the XYZ2grd based modelling in large-scale digital topographic models proposes better results in the context of cartographic aspect of contextual generalization. That is, the given identical input data, having two alternative methodologies to follow, the solution of XYZ2grd approach would be better.
The presented GMT based methodology demonstrated in this paper can be tried with other grid raster data sets covering other regions of the hadal trenches. Working on such hard-to-reach region as hadal trench, it is hard to collect field data for the assessment. Therefore, selecting reliable tools for the precise modelling of data is crucial for the correct understanding of such areas. It is commonly known that hand-made manual cartographic digitizing is a tedious and laborious work prone to errors. Automated digitizing of the profiles made by the machine with set up parameters provides reliable results and enables to perform geographic analysis of the remotely located region.
At the same time, the proposed sequential method of the raster data sets processing by means of the GMT is not fixed to the specifically marine areas and should also work for the mountainous areas. Data accuracy improvements may be achieved by applying other data sets (e.g. 15-sec GEBCO or SRTM). However, this should be compensated by the computer memory usage for the higher resolution data sets. Varying from location within the tectonic plate, the geomorphic conditions of the hadal trench are strongly affected by the geological conditions of the study area.
Developing an algorithm of the semi-automated digitizing of the geomorphic profiles aimed to visualize and model gradients of the slopes of the hadal trench and for comparing various different segments between two points provides more quantitative insights into the structure and geomorphic settings of the submarine geomorphology that still remains a question for the marine geology. The Kuril–Kamchatka has a high level of seismicity owing to a variety of the complex of geologic and geomorphic processes: volcanic hot spots, tectonic plates subduction, submarine erosion all affecting the trench geomorphology. The geologic constitution of Kamchatka region naturally determines the morpho-structural features of the Kuril–Kamchatka Trench and relief, as well as the processes of the formation of the submarine features (terraces, slopes, seamounts) closely related to the lithosphere activity. The geological evolution of the Kamchatka region is the result of the subduction of the Pacific and Okhotsk tectonic plates, separated by the Greater Kuril Chain, beneath the Eurasian Plate at the Kuril– Kamchatka Trench. Therefore, close correlation between the geology and tectonics controls the hydrological processes of the Sea of Okhotsk and geomorphology of the Kuril– Kamchatka area. In view of this, precise visualization of the geospatial data help to better understand the structure of the oceanic trenches, the least located natural phenomena on the Earth. It furthermore enables to analyse possible correlations between various processes and factors: tectonics, geology, sedimentation, oceanography affecting the landforms formation.
The importance of the effective cartographic approach of such geologically vulnerable zones as Kamchatka region can hardly be underestimated. Advanced GMT scripting toolset illustrated in the presented work included a variety of modules and techniques: data filtering, interpolation, gridding, projecting, Nearest Neighbour surface modelling, creating layouts, perspective 3D-mesh plotting and plain 2D modelling, visualizing contour maps, surface raster grid modelling with various colour palettes, artificially illuminated surfaces. In particular, a sequences of the GMT modules combined as scripts were used for contour mapping of the digital bathymetric data, raster grids visualization with sun-illumination and shading and xyz data modelling. Various colour palettes were tested as functionality of the GMT with the colour palette names explained. Besides cartographic instrumentation, GMT includes statistical data analysis and modelling trends by mathematic approximation. Therefore, as demonstrated, the GMT proves to be a powerful yet flexible cartographic toolset that enables a range of approaches in geospatial data modelling and geomatics. Based on the presented methodological approaches and results, GMT is strongly recommended for mapping and geoin-formation modelling.
