NOAA Marine Geophysical Data and a GEBCO Grid for the Topographical Analysis of Japanese Archipelago by Means of GRASS GIS and GDAL Library 2

: This article analyzes topographical and geological settings in the Japan Ar ‑ chipelago for comparative raster data processing using GRASS GIS. Data in‑ clude bathymetric and geological grids in NetCDF format: GEBCO, EMAG2, GlobSed, marine free‑air gravity anomaly and EGM96. Data were imported to GRASS by gdalwarp utility of GDAL and projected via PROJ library. Method‑ ology includes data processing (projecting and import), mapping and spatial analysis. Visualization was done by shell scripting using a sequence of GRASS modules: ‘d.shade’ for relief mapping, ‘r.slope.aspect’; for modelling based on DEM, ‘r.contour’ for plotting isolines, ‘r.mapcalc’ for classification, ‘r.category’ for associating labels, and auxiliary modules (d.vect, d.rast, d.grid, d.legend). The results of the geophysical visualization show that marine free‑air gravi‑ tational anomalies vary in the Sea of Japan (–30 to above 40 mGal) reflecting density inhomogeneities of the tectonic structure, and correlating with the geo‑ logical structure of the seafloor. Dominating values of geoid model are 3 0–45 m reflecting West Pacific rise, determined by deep density inhomogeneities as ‑ sociated with the mantle convention. Sediment thickness varies across the sea reflecting its geological development with density of 2 km in its deepest part and thinner in central part (Yamato Rise). The aspect map and reclassified map express gradient of the steepest descent. aly Grid EMAG2, GlobSed sediment thickness


Introduction
The current study presents GRASS GIS based research which entails the geoprocessing of several different raster grids and data analysis to visualize the regional geological and topographical settings in the Japanese archipelago: DEM based topographic modelling, geoid, marine free-air gravity, sediment thickness of the oceanic seafloor and magnetic anomalies. A special focus of this paper is placed on the methodological description of GRASS modules and snippets of scripting codes, together with visualized output maps demonstrating the cartographic functionality of GRASS GIS.
The main advantage of GRASS GIS is its varied functionalities: to process a wide variety of data which come in different formats (vector, raster netCDF, GRD, GeoTIFF). GIS plays a key role in the processing and visualizing of marine geospatial data, since it enables spatial analysis using advanced data processing and a variety of GIS approaches or statistical libraries with embedded mathematical equations, algorithms and models [1][2][3]. The proliferation of big data in the geosciences has brought a considerable amount of marine geological datasets to the fore and made it accessible for modelling and cartographic visualization.
Automatic data processing has actively been a question in cartography since the development of programming languages for the rapid processing of big data and methods for the automatization and optimization of data processing have constantly been progressing [4][5][6]. Therefore, another important asset of GRASS GIS consists in its applicability with Python and R, advanced high-level programming languages effectively used in geosciences [7,8]. The use of geostatistical models, algorithms of data processing and approaches of data visualization has had a long history in spatial analysis [9][10][11][12][13][14][15][16]. The field of geospatial modeling had made certain advances by the early 1970s, along with rapid development of the IT industry. Advances in computational modelling methods applied in GIS have helped to bring the power of data modeling to classic cartography.
The main research problem consists in highlighting relationships established between the geospatial datasets indicating a correlation between geographic variables showing the complexity of the geologic formation that depends on a set of impact factors. For instance, the impact of the tectonic development occurs at the stage of the relief formation (faults, ridges, deep-se trenches, fracture zones) which finally results in the geomorphic signatures of the modern relief. The seafloor morphology, geomorphic signatures and topographic slope orientation of the Japanese archipelago show remarkably different features in geographically distant regions of the Sea of Japan, Japan Trench and archipelago. Hence, the research problem is to highlight the correlations between these variables, such as Earth magnetic anomaly and marine free-air gravity, combined with sediment thickness, GEBCO relief bathymetry and geoid. Spatial data integrity may reveal correlation of variable geomorphic signatures within and between the datasets overlaid with slope/aspect map and reclassified topographic models. In this way, it can indicate what are the regional differences in geological, geophysical and tectonic settings and how does the topographic slope/ aspect distribution of the modern relief mirrors variations in the geologic settings.
The technical research problem consists of the integration of multi-source thematic datasets to combine different raster data grids into one coherent GRASS GIS project using a scripting approach by importing and transforming GDAL utilities (r.in.gdal, gdalwarp, gdal_translate), and performing numerical spatial analysis by a set of GRASS GIS modules (r.slope.aspect d.rast, d.shade, g.region etc). Aggregation of various high-resolution datasets undertaken in this study aimed to analyze geological, geophysical and topographic settings of the Japanese archipelago region for proper understanding of the deep-sea areas in the oceanic seafloor. There is a potential possibility that accurate orientation of the magnetic lineation zones, gravity anomalies, geologic fault lines, slabs, fracture zones striking in N-S, W-E slope orientation or various combinations of such orientations (e.g. NNW-SEE, NNE-SSW) correlate with geological lineations and their separation lines, and may help with geological prognosis and analysis. Accurate topographic modelling of the slope orientation (by r.slope.aspect and d.shade GRASS GIS modules) enables to compare distribution of the geophysical settings over study area. Different dynamic mechanisms of the seafloor spreading explain the phenomena of the variability of the submarine relief, continental rifting, slab pull and trench formation, during convergence of the Pacific Plate. According to gravimetric surveys and tectonic modelling, tectonic plate subduction plays an important role in trench formation and geophysical settings of the neighboring area. Therefore, a complex analysis of the multi-source datasets (by d.rast GRASS GIS module) clarifies the correlation between the geological and topographic variables.
The main goal of this study was to investigate if combination of diverse available spatial raster datasets including GEBCO terrain topography, Earth Magnetic Anomaly Grid EMAG2, GlobSed sediment thickness global grid, marine free-air gravity anomaly and geoid EGM96, processed by GRASS GIS allows to produce accurate thematic geological maps using scripting approach for complex geospatial analysis of the Japanese archipelago. Specifically, the study intended to produce a series of maps showing contemporary geologic settings in the Japanese Island area and the Sea of Japan at high spatial resolution based on high-resolution raster grids (one to five arc minute resolution) to allow consistent spatially explicit overlay of the topography, geopotential geoid grid, marine free-air gravity and magnetic anomalies, sediment thickness, and topographic modelling of slope and aspect based on raster data obtained from NOAA repository sources. To reach this goal, the study applied a scripting GRASS GIS algorithm to process several raster grids using GRASS code syntax, visualize and integrate various raster data layers in GRASS GIS environment through a sequence of commands. In the end, the study assessed topographic surface of the Japanese archipelago area through the reclassification of the initial topographic grid using algorithm of 'r.mapcalc' data processing and 'if-else' condition in the block of code followed by the 'r.category' operator for data analysis.
Finally, this work was meant to provide answers to the more specific research questions: -How does the functionality of GRASS GIS enable to model and map geographic objects based on the netCDF, GRD and IMG raster numerical grids, aimed to visualize geographically continuous phenomena as xyz datasets (magnetic anomalies, sediment thickness, geopotential grids, gravimetric anomalies)? -To what extent does it contribute to the development of the cartographic automatization with the ultimate aim of rapidly processing of the large datasets common in marine geology?
Answering these questions, the applied approach of GRASS GIS scripting algorithms enabled to create a series of the geological and topographical maps of the Japanese archipelago at the region scale, relating them to the phenomena of geospatial correlation between the topographical, geophysical and geological variables applying automatization techniques as advanced GIS methods.

Study Area
The study area is located in the Japanese archipelago and Sea of Japan within the square of coordinates 128°E-150°E, 30°N-46°N, Fig. 1. The Sea of Japan is characterized by a narrow shelf along the coasts of Korea, and Japan, with mean depths of 120-150 m [17]. The shelf break marks changes in continental slopes located at depths varying from 20 to 550 m (light blue, Fig. 1). The western coast of Japan is complicated by numerous underwater cliffs and small islands [18]. The continental slope in the northern part of the Sea of Japan is gentle, while along the coast of Korea is formed by a steep staircase-like structure with depths exceeding 3000 m. In the east of Korea there is a marginal plateau with depths of 1000-1200 m [19]. A submarine ridge extends in northward direction from the southern part of Honshu Island, reaching Yamato Rise in central part of Sea of Japan, with depths ranging from 300 to 1000 m. The seafloor of the Sea of Japan is flattened with occasional seamounts and slightly inclined to the north, which is contrasting with other seas of the Pacific Ocean: the Philippine Sea with complex and diverse seafloor relief, geological faults and fracture zones [20] and such prominent geomorphic features as the Philippine and Mariana trenches, containing the deepest point on the Earth [21]. The depths in the northern part of the basin exceed 3600 m, and in the southern part decrease up to 2600-2800 m.
The basement surface of the Mesozoic and Cenozoic structures of Japan descends from the coast to the seafloor of the basin as a steep slope with complex block division. The maximum depths of the basement seafloor in the Sea of Japan exceeds 5 km, and in its central part (Yamato Rise) 2-3 km. The geologic basement of the Japan Trench lies ca. 0.5-1.5 km deeper than its modern seafloor. The Japanese islands with the Japan Trench form a typical arc-trench system with developed continental structures. From the oceanic side, the Japanese archipelago is bordered by a narrow shelf with depths of less than 100 m, which smoothly continues to a slope of the deep-sea Japan Trench [22] (Figs. 1, 2). The profile of the slope is convex, with steepness gradually increasing downwards. The Japan Trench has a flat seafloor, with a maximum depth of 8,046 m according to the GMRT data [23] which is shallower when compared to the deepest trenches of the Pacific Ocean, e.g. Tonga and Kermadec trenches with depths exceeding 10,000 m [24]. The convergent plate margin of the seaward edge of the continental plate is deformed by subduction of the oceanic Pacific Plate. The northern Japan Trench is demarcated from the rigid continental plate as reported by a multichannel seismic survey [25].
The tectonic structure of the seafloor of the Sea of Japan, its continental shelf, Japan Trench and northwest segment of the Pacific Ocean are summarized in previous works [26][27][28][29][30][31]. The seafloor structure of the Japanese archipelago and its geological development have been detailed and reported using the datasets from regular cruise expeditions published as series of volumes "Initial Reports of the Deep Sea Drilling Project" [32][33][34].

Materials and Methods
The methodology of this work is based on the GRASS GIS (Geographic Resources Analysis Support System), an open source GIS, continuously developed by the GRASS Team since 1982 [35]. In contrast to traditional GIS, e.g. ArcGIS [36][37][38][39], GRASS GIS is shell-scripting based GIS with advanced functionality in coding and applicability to Python programming language and GDAL library [40]. GDAL was used for re-projecting raster grids into cartographic projection (Equal Area Cylindrical projection was selected for this research), warping, subsetting and converting data. The comparative analysis of several grids of datasets aims at highlighting correlation between the geographical variables visualized on the maps.
Data selection was performed based on the data quality: reliability (NOAA) and high-resolution. Spatial resolution of the input data is one of the key concepts in GIS since the Earth's complexity restricts studies in full detail. Resolution is critical in determining a dataset's suitability for a given research problem, as it affects the minuteness at which objects can be mapped and interpreted during spatial analysis. Technically, the resolution determines the volume of the processed geospatial grids which affects computer memory and data storage volume.
The mapping is done using the raster grids visualized by GRASS GIS through a sequence of modules using three main methodological blocks: 1) data preprocessing, 2) data visualization (mapping), 3) spatial analysis.
The coding methodology for each of these blocks is explained below stepwise.

Data Preprocessing
Data preprocessing included the retrieval of information from the different raster grids aimed at highlighting coherency between geological, topographical and geophysical grids. Due to the differences between the original datasets in terms of content and consistency, it included the analysis of their metadata, spatial characteristics and structure: projection, datum, extent, range, spatial resolution.
The multi-source data origin resulted in their various resolutions: 30 arc-second GEBCO, 2 arc-minute EMAG2, 5 arc-minute GlobSed, 15 arc-minute grid EGM96, 1 arc-minute marine free-air gravity. Thus, data preprocessing in GRASS included the following four steps: 1. The first step involves cartographic projecting of the initial raster NetCDF in WGS 84 datum surface imported to GRASS GIS and warping it to an Equal Area Cylindrical projection by 'gralwarp' utility of GDAL using PROJ Cartographic Projections and Coordinate Transformations Library [47] via GMT [48]. The projecting was done using following code (here is the example for the NetCDF grid of topographic layer (Figs. 1, 2), repeated for other thematic grids (Figs. 3-6), respectively): gdalwarp -t_srs '+proj=cea lat_ts=38 lon_0=138' jt_relief.nc jt_relief_EAC.nc. 2. The second step involves the re-projecting of the raster grids read into the GRASS GIS project using 'r.in.gdal' module by following code: 'r.in.gdal jt_relief_EAC.nc out=jt_relief_EAC title="Topography: GEBCO grid"'. 3. The third step involves applying date and time function of GRASS GIS: r.timestamp map=jt_relief_EAC date='04 May 2020'. 4. The list of the raster grids and information (projection, datum, min/max range) was received by modules 'g.list rast' which visualizes the available raster files in the current project and 'r.info'.
Vector data import from the GMT was done by 'v.in.ogr' command: v.in.ogr trench.gmt out=trench. The compatibility of the GRASS GIS with native GMT formats is an advantage of the GRASS GIS, as GMT is a functional cartographic toolset often used in geophysical mapping [49][50][51].

Cartographic Mapping
Topographic dataset was processed based on the isolines generalization highlighting the particularities of the relief. Data at different levels of representation containing specific information (resolution and smoothness of magnetic and gravity anomalies, topography, geoid, sediment thickness) were integrated in a GRASS GIS project.
The resulting cartographic mapping and data visualization included following eight steps: 1. Displaying a raster map by initialized GRASS GIS monitor: d.mon wx0 and creating spatial subsets of the region of Japanese archipelago using the 'g.region module': g.region raster=jt_relief_EAC -p 2. Applying color palette for the maps via 'r.colors' GRASS module. Following color tables were selected for the maps: 'srtm_plus' (Fig. 1), 'haxby' [52] (Figs. 2, 6), roygbiv (Fig. 3), celsius (Fig. 4), rainbow (Fig. 5), aspectcolr (Fig. 7) and defined colors by 'rules' function for Figure 8. Color assignment was done by code: r.colors jt_relief_EAC col=haxby 3. Spatial analysis of slope, shaded relief and plotting isolines were performed using following modules: 'r.slope.aspect', 'd.shade', 'r.contour' and 'd.vect', as follows in the GRASS code: r.slope.aspect elevation=jt_relief_EAC aspect=jt_aspect_EAC d.shade shade=jt_aspect_EAC color=jt_relief_EAC (the 'd.shade' module resulted in overlay of the topographic DEM and aspect models, in which raster cell values are equally spaced, gridded according to a matrix structure of topography and aspect which gives the effect of shaded relief visualization (Fig. 2). r.contour jt_relief_EAC out=TopoJT_EAC step=2000 -overwrite d.vect TopoJT_EAC color='75:23:3' width=0. 4. The next step included plotting a vector border box around the study area (red line on Figure 1) by module 'v.in.region' using following code: v.in.region output=jt_EAC_bbox v.info map=jt_EAC_bbox d.vect jt_EAC_bbox color=red width=3 fill_color="none". 5. Cartographic projected coordinate grids along the map borders were plotted using following GRASS code: d.grid size=5 color='172:219:250' border_color=yellow width=0.1 fontsize=8 text_color=red bgcolor=white. 6. Legend to the left of each map was plotted using following code where self-explaining cartographic attributes, data range and visual settings of the legend are plotted by flags, such as font size, color, border, background color (here is example for the topographic map ( Fig. 1): d.legend raster=jt_relief_EAC range=-9759.701,3700.742 title=Topography,m title_fontsize=8 font=Arial fontsize=7 -t -b bgcolor=white label_step=1000 border_color=gray thin=8. 7. Annotating textual information was done by following code example: d.text text="Honshu Island" color=green size=2.5 font="Trebuchet MS". 8. Adding the title of the maps was piped by Unix redirection combining 'd.title' and 'd.text' GRASS modules using following code example: d.title map=jt_relief_EAC | d.text text="Shaded topography: Japan Islands" col-or="green" size=3.

Spatial Analysis
After the raster grids were visualized (Figs. 1-6), the next step included topographic spatial analysis (Figs. 7, 8). The spatial analysis performed on the topographic dataset can be considered as a derivation process of the initial reference grid (GEBCO) aimed to retrieve information on the slope and aspect of the topography. The spatial analysis was made by following steps: 1. The compass orientation was generated and classified by four major directions (N, E, S, W) using adopted code [53] developed based on the algorithms of data analysis [54][55][56] The resulting map re-classified from Figure 7 is presented on Figure 8.

Results
The topographic map (Fig. 1) and derived shaded relief map (Fig. 2) are based on the GEBCO global bathymetry and topography grid with a spatial resolution interval of 15 arc seconds where bathymetry is produced using a combination of shipboard and satellite altimetry predicted data.
The map of the geoid in the Sea of Japan (Fig. 3) is based on the visualized EGM96 raster grid. It shows the transition zone of geoid undulations from the Pacific Ocean to the continental area. Composition relationships in geophysical datasets are based on the analysis of several sub-entities (marine free-air gravity anomaly values distribution correlated with geoid and topographic isolines over the study area). Thus, the region of the Japan Archipelago is located in the vast West Pacific rise of the geoid (the long-wave part of the spectrum of the geoid). Therefore, as can be seen in Figure 2, the dominating values are 30-45 m (orange to red colors on the map), followed by 26-30 m (yellow colors). Lower geoid values (17-25 m, green and 10-17, cyan) are mostly repeating topographic depressions of the Pacific Ocean and the lowest −5 to 10 (blue to purple) corresponding to the Japan Trench topography (Fig. 3). The nature of the geoid is determined by deep density inhomogeneities that may be associated with the mantle convention.
Map of the marine free-air gravity (Fig. 4) shows gravimetric patterns in the basin of the Sea of Japan. There are weak positive anomalies (0-5 mGal, cyan colors, Fig. 4) or negative anomalies (0 to −17 mGal, dark blue colors, Fig. 4) and in some depressions, larger negative values (−17 up to −30 mGal, purple colors, Fig. 4) are observed. The Yamato Rise (in the central part in the Sea of Japan) is distinguished by positive Faye anomalies (>20 mGal, orange to red colors, Fig. 4). Above the Japanese Islands, the marine free-air gravity anomalies reach values above 40 mGal and higher (grey-pink colors, Fig. 4). They sharply decrease to a deep-sea trench, up to strongly negative values (below -30 mGal, white color, Fig. 4).
The marine free-air gravitational anomalies field varies significantly in the Sea of Japan (from −30 mGal to above 40 mGal) reflecting density inhomogeneities of the tectonic structure, but in general they correlate with geological structures of the seafloor. On the shelf they are mostly positive (20-40 mGal, orange to red colors). The shelf area is characterized by linearly elongated anomalies parallel to the shelf isolines. In the Japan Trench, they reach values lower than -30 mGal. An accurate and consistent visualization of the marine free-air gravity field in the Japan region is crucial for geophysical studies of western Pacific Ocean and thus also for a better understanding of the geological structures of the Japanese archipelago. The map of the magnetic anomalies of Japan (Fig. 5) shows the magnetic anomalies zone, which is one of the three magnetic anomalies zones of the Pacific Ocean forming a system of ancient magnetic anomalies which reflects the basin development in Mesozoic period (the other two are Phoenix and Hawaiian). Theoretical basis of the magnetic anomalies of the oceans is described in the existing studies [57,58] and explained by regional particularities of the magnetic lineament anomalies [59] or anisotropy of magnetic susceptibility [60]. The geometric relationships analyzed from the comparison of datasets on magnetic anomalies, topography, geoid, and marine free-air gravity are meant to translate spatial relationships between the geographical and geophysical entities. These are especially concerned with the notions of the geometric similarities notable on data value distribution (depicted as colored areas) and isolines (depicted as lines), the intersection as overlay in a GRASS GIS project and spatial connection between the data. Thus, the Japanese system of the magnetic anomalies occupies the western and northern parts of the Northwest Pacific Basin which correlates with the topographic subdivision of the basin based on the elevation distribution. The magnetic anomalies here are oriented as strips in a northeast direction (Fig. 5). The magnetic anomalies of Japan region are separated from the Hawaiian system by the Shatsky Rise which is the 3 rd Earth's largest oceanic basaltic plateau formed during a period of geomagnetic reversals [61,62]. The magnetic anomalies are then 'absorbed' by the Japanese Trench. The amplitudes of the magnetic anomalies exceed 400 gammas in their highest values (Fig. 5, red colors).
The correspondence between the dataset of sediment thickness and topographic relief represent two closely correlated phenomena of marine geological settings: sedimentation and bathymetry. Map of sediment thickness (Fig. 6) shows distribution of the marine sediments as important variables in geological datasets. Dependencies in sedimentation distribution can be clearly observed comparing two maps (Figs. 1, 6) showing higher sedimentation values in shallow shelf areas and river mouth estuaries and almost absent or low values in the open ocean basin.
Visualization of sediments reflects seafloor sedimentation records the Earth's tectonic movements during geological history, and may assist in reconstructing paleoclimate data. Visualized sediment maps can assist in highlighting correlations of oceanological parameters, such as ocean current directions and strength, changes in sea level and sea temperature. Besides, seafloor sedimentation is important for studies of the marine biological productivity and variations in environmental conditions, biological evolution and ecology since it indirectly reflects the complexity of the climatic and geological processes [63]. Therefore, it has always been important to visualize seafloor sediment coverage for geological studies.
Sediment coverage in the Japan Trench region is contrasting with other trenches of the Pacific Ocean, e.g. Middle America Trench having Quaternary trench-filling turbidites, fault scarp in lower Miocene chalk, pelagic clays and diatom ooze [64], or Mariana Trench with sedimentary cover of basalt of the leitic composition, a complex of parallel dykes of diabases, an isotropic gabbro, a banded gabbro-ultrabasic complex [65]. In the deepest part of the basin of the Sea of Japan, the sediment thickness represented by Meso-Cenozoic sediments is about 2 km. In its central segment on the Yamato Rise and other the sediment coverage is thinner and has a sharply variable thickness: it decreases to a minimum on the individual peaks and steep slopes and increases in local depressions and saddles [66]. Some data received from drilling in the Yamato Rise discovered 529 m thick layer of siliceous-siliceous silts, gradually changing to diatomaceous silts of the Upper Miocene, followed by volcanic sand and tuffs. Further reading regarding regional patterns of the sedimentation at the seafloor of the Pacific Ocean and theoretical processes of sediment deposition through modelling and observations are presented recently [67,68] and in earlier works [69].
The dependency relationships were established based on the analysis of the datasets visualized as several raster grids. They were used to compute the attribute values of the marine free-air gravity, topographic grids, sedimentation and geoid data with slope/aspect reclassified map. For geographic databases, the geometry was considered as an attribute for data comparison. The topographic spatial analysis is presented on Figures 7 and 8. Figure 8 shows the results of the reclassification of a continuous field map showing aspect of the topographic slope in the Japanese archipelago (Fig. 7).
The initial aspect angle ranged from 0 to 360 degrees with the origin in the east and a counter-clockwise rotation. The reclassified map (Fig. 8) shows the aspect map with values grouped into four aspect clusters reflecting the terrain compass orientation (North-South-West-East). The topographic parameters describe the geometry of the terrain and seafloor surface in a given point, such as slope and aspect. Estimating and visualizing a slope aspect and steepness is useful for the analysis of such hazardous geological parameters as soils creep, debris avalanching and landslides [70]. Transport processes on the slopes include surface processes (e.g. flow) Fig. 7. Aspect map of the slope topography: Japanese Archipelago region. Mapping done using GRASS GIS, based on GEBCO elevation raster grid and subsurface process (leaching) for middle-scale spatial analysis and large mass movements that are the subject of gravity tectonics, for a small-scale analysis.
A geometrical arrangement of patterned ground on the slopes is represented by geomorphological landforms such as polygons, local depressions, circles, steps etc. Slope analysis is essential for the geomorphological viewpoint, since both terrain and submarine ground surfaces are composed of the landform elements related to slope: curved geometric shapes, related to upslope, downslope, and orientation. Hence, altitude, extent, slope, curvature, ruggedness, and aspect are essential parts of the spatial analysis based on DEM [71]. The slope of the terrain surface allows gravity to induce the flow of materials (debris, rocks) and water. Therefore, it lies at the core of many geological models.
The aspect of the slope is a compass direction of the maximum rate of variation. Continuous classes in the topographic aspect modelling were first visualized (Fig. 7) and then modeled based on the 'r.mapcalc' re-classification algorithm applied to the initial dataset of GEBCO, merged in four classes in the derived grid by means of the attribute re-calculation (Fig. 8). The reclassified maps express the orientation of the gradient of the steepest descent, converted to a compass attributes N-S-W-E. Comparison of the map on Fig. 8 with topographic DEM (Figs. 1, 2) shows that slope aspect in the Japanese archipelago region corresponds to the general topography. Mapping done using GRASS GIS (modules r.slope.aspect, r.mapcalc, r.category), based on GEBCO raster grid

Conclusion
A unique aspect of the submarine geomorphology comparing to the terrestrial studies lies in the fact that seabed is not available for direct observations. Marine geological datasets could only be received by the remote sensing methods of the seafloor survey and data post-processing using GIS. In view of this, the selection of the advanced GIS tool for data processing is clear. The GRASS GIS, as demonstrated in this paper, presents an excellent tool for the processing of the raster datasets enabling to visualize, perform spatial topographic analysis and map thematic grids.
The scripting approach provided by GRASS GIS affords the opportunity to combine a number of thematic raster grids and visualize layers using selected code lines of GRASS GIS syntax. Such a collection of codes can be then saved as shell scripts. Using shell scripts in GIS, rather than plotting maps through a Graphical User Interface (GUI), provides significant editing and analytical possibilities due to the principle of repeatability. Shell scripting presents a new way in cartographic automatization through visualizing spatial big datasets in a rapid and effective way. In this sense, this paper presented a contribution to the development of the cartographic methods with a case study of complex analysis of the marine geological settings of the Japanese archipelago by raster datasets.
The GRASS GIS scripting-based mapping relied on using a specific code syntax used for processing a series of raster datasets, as demonstrated in this research. Such a technical approach may use the advantage of the repeatability of scripts and perform a series of maps with identical spatial extend and coverage reflecting correlations between topographic setting, geophysical anomalies and geologic conditions of the study area. As one of the goals driving this research was the visualization of the series of the thematic maps by scripts written on GRASS GIS syntax, such a paradigm demanded a careful evaluation and standardization of the input multi-source data (projections, resolution, spatial extent, coverage) to reduce the related uncertainty of the output spatial series of the geological, geophysical and topographical maps.