[PYTHON] GSI_DEM to geotiff conversion → UTM conversion → ascii conversion only on Ubuntu command line

I want to convert the Geospatial Information Authority of Japan DEM to a format that is easy to use with GIS only on the command line of Linux (Ubuntu).

If there is something like an elevation map or terrain shading map read from a digital elevation model (DEM) as a base for placing GIS information, the GIS information will look good. In the case of Japan, the Geospatial Information Authority of Japan (GSI) publishes 5m and 10m mesh DEMs as national land infrastructure information. GSI DEM is provided in the form of the provided JPGIS (GML) format. When you open it with a text editor with something like Metadata, you will find information about the data such as the pixel size and latitude / longitude of the DEM at the top of the file.

Personally, it is very convenient if it is in geotiff format, but if you take a quick look at the format conversion to geotiff format, software that converts in the Windows environment will hit. However, as a Linux user, it is still troublesome to start up external software and convert the file format.

Here, I showed how to complete the conversion to geotiff format + alpha format conversion only on the terminal from the place where GSI DEM was downloaded.

This time is an example of bash environment. A detailed explanation is given below the script.

gsidem2geotiff.sh

#!/bin/bash

#Added what was missing in my environment.
#If there is something missing, add it yourself.
#sudo apt-get install git osmium-tool gdal-bin python python-pip python-numpy python-gdal python-matplotlib python-beautifulsoup python-lxml
#pip install gdal 
#git clone https://github.com/minorua/fgddemImporter.git

#unzip GSIDEM zip file
unzip PackDLMap.zip

#GSI DEM->geotiff #GSI DEM has been DL FG*As a zip.
python fgddemImporter/fgddem.py FG-GML-*zip

#geotiff file mosaic(merge.tif).. The value of the sea area-From 9999 to 0.
#-9999 (-srcnodata specification)Value of 0(-dtsnodata specification)To.
gdalwarp -of "ENVI" -srcnodata -9999 -dstnodata 0 FG*tif merge.tif

#transform EQA2UTM
#EPSG code| EQA:4326, UTMzone52:32652, UTMzone53:32653, UTMzone54:32654 
gdalwarp -s_srs EPSG:"4326" -t_srs EPSG:"32652" merge.tif merge_utm.tif

#Show the contents of geotiff
gdalinfo merge.tif

#I don't think it's necessary, but geotiff->Converted to ESRI ascii format.-Let 9999 be NaN.
#-It is cut out with projwin. Edge coordinates clockwise from west(Meter)Is specified.
gdal_translate -projwin ${west} ${north} ${east} ${south} -a_nodata -9999 -of AAIGrid merge.tif dem.asc

#EOF

The important thing is ** Conversion of DEM file downloaded from Geographical Survey Institute from zip file to geotiff, conversion of geotiff to mosaic, ESRI ascii format (cut out) is on the terminal (this time Ubuntu18.04LTS, bash environment) It can be completed with **.

Practical example-I want to merge GSI DEM and cut out only a part of the area-

Let's take Sakurajima as an example this time. It is a review of the procedure.

1. 1. Download DEM file from GSI HP
2. Convert from JPGIS (GML) format to Geotiff format
3. 3. Convert to UTM coordinates Four. Cut out by specifying UTM coordinates and convert to ESRI ascii format

This time, we will use the data in the area below.

I just put a concrete value to specify the cutout area in the above script. The cutout area this time is UTM coordinates (west, north, east, south) = (650500 3515000 665300 3490000).

Result (drawn in QGIS)

For the time being, as requested, I was able to cut out only the area where I wanted to bring the GSI DEM in a wide area and convert it to UTM coordinates. Here, the gdalwarp command is used to convert mosaic and UTM, but if you don't need the mosaic geotiff file of EQA coordinates, you can convert the mosaic UTM coordinates in one shot by connecting the options.

gdalwarp -s_srs EPSG:"4326" -t_srs EPSG:"32652" -of "ENVI" -srcnodata -9999 -dstnodata 0 FG*tif merge_utm.tif

Cut out by specifying EQA coordinates first → Try UTM coordinate conversion

Earlier, EQA-> UTM coordinate conversion-> UTM coordinate specification area was cut out, but if you dare to cut out in the EQA area-> UTM coordinate conversion, the appearance will change a little.

gsidem4sakurajima.sh

#!/bin/bash

#GSIDEM2geotiff
python fgddemImporter/fgddem.py FG-GML-*zip

#mosaic geotiff files -> merge.tif
gdalwarp -of "ENVI" -srcnodata -9999 -dstnodata 0 FG*tif merge.tif

#show geotiff information
gdalinfo merge.tif

#Extract Sakurajima region
gdal_translate -projwin 130.5870 31.6389 130.7414 31.5371 -a_nodata -9999 merge.tif merge_ex.tif

#transform EQA2UTM
#EPSG code| EQA:4326, UTMzone52:32652, UTMzone53:32653, UTMzone54:32654 
gdalwarp -s_srs EPSG:"4326" -t_srs EPSG:"32652" merge_ex.tif merge_utm.tif

#show merge_utm.tif information
gdalinfo merge_utm.tif

#Convert geotiff2ascii
gdal_translate -of AAIGrid merge_utm.tif dem.asc

Result (drawn in QGIS)

Latitude / longitude specified cropping → EQA to UTM coordinate conversion, so the coordinates were converted correctly, but the area was rotated counterclockwise as it looked. If you want to use this further, it may be a little difficult to use.

Finally

It feels awkward to load a python script with a shell script. You should be able to do it well with python alone ...