Okinawa Trough Geophysical and Topographic Modeling by GDAL Utilities and GRASS GIS

This paper presents using GDAL utilities and GRASS GIS for topographic analysis of the raster grids based on GEBCO DEM as NetCDF file at 15 arc-second intervals. The focus study area encompasses the area around Okinawa Trough, Ryukyu trench-arc system, southern Japan, East China Sea and the Philippine Sea, west Pacific Ocean. Several GDAL utilities were applied for data processing: gdaldem, gdalwarp, gdalinfo, gdal_translate. The data were imported to GRASS GIS via r.in.gdal. Data visualization highlighted high resolution and accuracy of GEBCO grid, enabling topographic modelling at the advanced level. The algorithm of DEM processing, implemented in GDAL utility gdaldem, was used for generating multi-purpose topographic models: aspect, hillshade, roughness and topographic indices, such as Topographic Position Index (TPI), Terrain Ruggedness Index (TRI). Thematic maps (topography, geoid, marine free-air gravity) were visualized using GRASS GIS modules for raster (d.rast, r.colors, r.contour) and vector (d.vect, v.in.region, d.legend) data processing. The results demonstrated smoother bathymetry in the East China Sea and rugged relief in the Philippine Sea which corresponds to their different geological and geophysical settings. The presented methodology of the topographic analysis based on DEM demonstrated technical aspects of GDAL and GRASS as scripting approach of advanced cartography.


INTRODUCTION
This paper presents application of scripting cartographic methods for topographic data analysis and modelling by GDAL, Geospatial Data Abstraction Library (GDAL/OGR contributors, 2020) and open-source GRASS GIS (Neteler, 2001;Neteler and Mitasova, 2008). Automated data processing in geological mapping is presented in the existing literature (Gruber et al., 2017;Maggiori et al., 2017;Krčmar et al., 2020;Bongiovanni et al., 2021;Thiele et al., 2021;Arriagada et al., 2021;Balogun et al., 2021). Using scripts and console based methods in cartography presented a new step in further developing of cartographic methods and solutions for machine learning and automatization in GIS. Data processing by GDAL followed by geo-visualization in GRASS GIS, presented in this paper, aimed at topographic modelling for regional physical geographical analysis in the area of the Okinawa Trough and Ryukyu trencharc system. Such methodology has not been qualitatively explored in the existing papers, despite some thematic overlaps (Ohta et al., 2013;Minami and Ohara, 2020;Arai et al., 2016a;Minami et al., 2020;Li et al., 20201;Xu and Chen, 2021). The aim of this paper is to perform a topographic analysis of the selected study area by applying various types of topographic modelling: Topographic Index Ruggedness (TRI), Topographic Position Index (TPI), roughness, hillshade, slope orientation aspect, as well as visualizing geoid and gravimetric grids using GDAL utilities (gdaldem, gdal_translate, gdalwarp, gdal_contour) and GRASS modules ( Figure 1).
The study area includes region of southern Japan: Okinawa Trough and Ryukyu trencharc system, located between the East China Sea (ECS) and the Philippine Sea, Pacific Ocean (Ludwig and Murauchi, 1973), 120°-134°Е, 20°-33°N. The area is notable for the complex topography and variability in the depths: the ECS is shallower comparing to the Philippine Sea. The Okinawa Trough and Ryukyu trench-arc system, stretching parallel, are notable for the rugged bathymetry and local spatial variations. Various factors affect submarine topography, which is reflected by complex structure and local variations in the marine environment (Lemenkova, 2019) highlighting correlation in a complexity of geophysical and geological factors: dynamics of the plate tectonics, Electronic copy available at: https://ssrn.com/abstract=3888343 magmatism (Shinjo and Kato, 2000;Kiminami et al., 2017;Guo et al., 2017), seismicity (Klingelhoefer et al., 2009;Konstantinou et al., 2013), volcanism and sedimentation (Lee et al., 2004;Huh et al., 2006;Minami and Ohara 2016;Yang et al., 2021), subduction channel thickness, landforms (Oiwane et al., 2011;Lemenkova 2020aLemenkova , 2020b, climatic factors (Hayakawa et al., 2015).
Topography in the region of Okinawa Trough and Ryukyu Trench is uneven due to the distribution of the three groups of islands on the forearc basement stretching over the Ryukyu arc (Kizaki, 1986) with Kerama Fault separating southern and central groups of islands. The seafloor of the ECS is well levelled with the depth ranges mostly at 60÷90 m rarely increasing up to 100÷120 m at the outer edge. In the southeastern part of the sea, narrow Okinawa Trough basin extends along the Ryukyu Arc, with depths increasing southwards (800÷2900) m. The Ryukyu Arc stretching parallel to the Ryukyu Trench and Okinawa Trough ( Figure 2) is an arch-shaped gently curved 1000 km long chain of mountain structure. The surface of the Cenozoic structures of the Ryukyu Arc stepwise descends to the seafloor. Geomorphic variations differ regionally: gentler in central and south (40÷50°) at 7 km depth and steeper (70°) at >7 km in the north (Long and Hilst 2006). Seismic investigations of the Okinawa area shown that sediment thickness is ca. 2-3 km on the oceanic plate, crust thickness is 5÷6 km. Upper mantle material underlying the oceanic plate has low seismic velocities, due to the geochemical processes of the mantle peridotites (Klingelhoefer et al., 2012). Trenches located far away in the ocean have dominating 'pelagic' sedimentation type with pyroclastic material coming from the neighbour island arc (Yap, Palau, Mariana, Tonga), while Ryukyu Trench and Okinawa Trough located near the shelf areas (ECS) and large islands -'terrigenous' type with evident sediment filling. In contrast to the ECS, the Philippine Sea, located southwards from the trench, is marked by a small amount of sediments and deeper bathymetry ( Figure 2).
Topography of the seafloor is affected by the rates of the sedimentation, strength of ocean currents and lithological composition of the seafloor rocks (Arai et al., 2016b;Maunde et al., 2021). The sedimentation in the Okinawa region is as follows: in its eastern part, sediment thickness is 0.3-0.5 km, increasing to 1 km in the northern part, decreases to 0.1-0.2 km in the western part and increases near Taiwan and Ryukyu (Litvin, 1987).

Figure 3
Marine free-air gravity grid (right). GDAL commands in GRASS GIS: 'gdalwarp' for reprojecting grid, 'r.in/gdal' for data import and visualizing metadata by 'r.info' (left). Source: author Important factor affecting the topography in the Okinawa region includes geophysical settings. Thus, regional differences in the magnetic anomalies in the ECS and the Philippine seas near Okinawa Trough and Ryukyu trench-arc system are caused by regional variations in the geochemistry and magnetism of the underlying seafloor rocks and Earth's magnetic field. Linear magnetic anomalies near the Ryukyu Arc are caused by the Antiklinorium of the alpine arc structures and Precambrian anticlinal structures on the shelf of the ECS. Philippine Sea is notable for a typical oceanic linear magnetic field, consistent with the oceanic crust. Regional seismic differences can be illustrated by waves velocities: granite layer on the shelf 5.8-6.3 km/s, under the Ryukyu Arc -5.0-6.2 km/s, in the basaltic layer -6.7-7.2 km/s (Litvin 1987).

METHODOLOGY
GDAL utilities and GRASS GIS commands were applied to the study area for topographic modelling. The methodology includes following steps. First, the topographic raster NetCDF file (rt_relief_UTM52.nc) initially in WGS84 geographic reference coordinate system was warped to UTM proj Zone 56 by GDAL using PROJ library via 'gdalwarp' utility. The file was then imported to the GRASS GIS environment via (Figure 1, left). The metadata of the output file were checked up by GDAL 'gdalinfo' utility ( Figure 1, right): longitude of central and prime meridian, spatial reference details, datum, projection. The same procedure was repeated for other raster grids: geoid (EGM96) and marine free-air gravity which shows difference between the observed acceleration of gravity on the Earth's surface and corresponding modelled value predicted from the Earth's gravity field.  (Figure 4, left). The conceptual approach of GRASS GIS consists in scripting cartographic modelling which means that the maps is being built using a sequence of modules used for each element of map and/or algorithms of data processing. For example, plotting contour lines for gravity map was perform by code using 'r.contour' module: 'r.contour rt_grav_UTM52 out=GravRT_UTM52 step=30', where step=30 assigns an interval for isolines and 'out=GravRT_UTM5' flag stands for defining output file.  Figure 2, 3 and 4.
The next step included topographic modelling using GDAL utility 'gdaldem', which was used for plotting a series of maps for topographic analysis: Figures 5-7 using sequence of five GDAL commands in a script ( Figure 5, left). The mathematical algorithms of the modelling are based on the existing GDAL technical documentation (GDAL/OGR contributors, 2020). The aspect map ( Figure 5, right) was generated from the DEM in NetCDF format (rt_relief_UTM52.nc) using gdaldem utility as follows: 'gdaldem aspect -of GTiff rt_relief_UTM52.nc rt_aspect.tif -trigonometric -zero_for_flat -alg ZevenbergenThorne'. The input raster was geographically North-oriented (UTM Zone 56). The output map (rt_aspect.tif) is stored as a 32-bit float GeoTIFF raster with orientation values 0°÷360° indicating azimuth that slopes are facing. Because the trigonometric angle was used instead of azimuth, the orientations were set up as follows: 0° is an East facing slope, 90° a North facing slope, 180° is a West facing slope and 270° is South-facing slope, respectively. The flag 'zero_for_flat ' assigns 'zero' for absolutely flat areas (slope=0°). The hillshade visualizing the terrain (Figure 6, left) was generated by 'gdaldem' utility using following code: 'gdaldem hillshade -of GTiff rt_relief_UTM52.nc rt_hillshade.tif -az 315 -alt 45 -alg ZevenbergenThorne -b 1 -z 1 -combined'. Here the azimuth was defined as 315° and altitude of the sun (artificial light source). The Zevenbergen and Thorne formula was applied, instead of Horn's formula (-alg ZevenbergenThorne'), to compute aspect, since Zevenbergen and Thorne visualizes better smooth terrains. The 'z' flag was used to set up vertical exaggeration as one, to avoid multiply the elevations (natural view). The '-combined' flag was applied to combine slope and oblique shading. The output is stored as an 8-bit GeoTIFF (rt_hillshade.tif) with a shaded relief effect (Figure 6, left).

RESULTS AND DISCUSSION
Roughness of topography ( Figure 6, right) shows a key geomorphometric variable widely used in geosciences (Hobson 1972;Frankel and Dolan, 2007;Grohmann et al., 2009). Topographic roughness indicates variability of elevation values, enabling to identify landforms and geomorphic processes. The roughness was plotted using embedded algorithms of GDAL using variability elevation DEM values from the initial NetCDF grid. The algorithm consists in computing the area ratio through assessment of similarities between the surface and flat area of square cells. Thus, roughness expresses standard deviation of the elevation values from a best-fit plane. Roughness was computed by following GDAL command: 'gdaldem roughness rt_relief_UTM52.nc rt_roughness.tif'. The Terrain Ruggedness Index (TRI) was defined by 'gdaldem' utility as the mean difference between a central pixel and its surrounding cells (Figure 7, left). The TRI was used to quantify topographic heterogeneity of the Okinawa Trough region. The calculating of the TRI was done by GDAL for each point on the grid: the algorithm computed average values across all grid cells as an elevation terrain ruggedness of the surface area. The computation was done using code: 'gdaldem TRI -compute_edges -of GTiff rt_relief_UTM52.nc rt_TRI.tif'. Here the '-compute_edges' flag was used to compute cells at raster edges and 'near-nodata' values. The output map is presented on Figure 7 (left). As can be seen, the most rugged areas (TRI 500÷725) are coloured dark brown and correspond to the mountainous areas and deep-sea Ryukyu Trench. The TRI of 275÷325, yellow-coloured areas, correspond to the moderate ruggedness, and the lowest values (TRI 0÷50) are notable for the near-flat well-levelled areas of the seafloor (well seen for the ECS seafloor, while more rugged for the Philippine Sea seafloor with TRI 50÷275).
The Topographic Position Index (TPI) is computed using DEM elevation values. The concept is based on the comparison of the elevation of each cell to a mean elevation of a neighbour cell. The TPI was generated using following GRASS GIS code: 'gdaldem TPI -compute_edges -of GTiff rt_relief_UTM52.nc rt_TPI.tif'. The output map shows TPI for the Okinawa area (Figure 7, right). Positive TPI values (Figure 7, right, green to brown colours) show locations higher than the average of the nearby area. Vice versa, negative TPI values (Figure 7, right, blue to black colours) represent locations lower than the surroundings. The TPI values near zero are flat areas (the slope is almost zero) or areas of constant slope. As can be seen, the zero TPI values (light blue colours) are characteristics for the most of the seafloor, negative TPIfor the deep-sea trenches and troughs (Ryukyu Trench, Okinawa Trough) and positive TPI are noticeable for the terrestrial areas. High TPIs can be noticed on the Taiwan Island (brown colours) and around the Ryukyu Arc.

CONCLUSIONS
GDAL and GRASS GIS scripting techniques in cartography demonstrated in this research presents a new step in cartographic development with an accent on machine learning and application of advanced scripting and coding in GIS. Applications of GIS in geoscience are diverse and well known: ArcGIS and QGIS (Klaučo et al., 2013(Klaučo et al., , 2017Jamieson and Stewart, 2021;Tsunogai et al., 2002;Ruan et al., 2020;Fang et al., 2020;Subarno et al., 2016), GMT. However, the use of big data and advanced programming tools in geoscience (Lemenkov and Lemenkova, 2021) requires use of console-based data processing rather than GUI.   Geology, 231(1-4)