국토 지리원 DEM을 Linux (Ubuntu) 명령 줄에서만 GIS에서 사용하기 쉬운 형식으로 변환하고 싶습니다.

GIS의 정보를 싣기 위한 베이스로서 수치 고도 모델(DEM)에서 읽은 고도도와 지형 음영도 같은 것이 있으면, GIS의 정보가 비쳐 보입니다.
일본의 경우는 국토지리원(GSI)이 국토기반정보로서 5m, 10m메쉬의 DEM을 공개하고 있습니다.
GSI DEM은 제공하는 JPGIS (GML) 형식의 형태로 제공됩니다.
Metadata와 같은 텍스트 편집기에서 열면 파일 상단에 DEM 픽셀 크기 및 위도 경도와 같은 데이터에 대한 정보가 들어 있습니다.

개인적으로는 geotiff 형식으로 되어 있으면 매우 편리합니다만, geotiff 형식으로의 포맷 변환에 관해서 확실히 조사하면, Windows 환경에서 변환해 주는 소프트등이 히트합니다.
그러나 리눅스 사용으로는 일일이 외부 소프트웨어를 시작하여 파일 형식을 변환하고 있어도 역시 귀찮습니다.

여기에서는 GSI DEM을 다운로드해 온 곳에서 geotiff 형식으로의 변환+알파의 형식 변환을 터미널상에서만 완결하는 방법을 나타냈습니다.

이번은 bash 환경의 예입니다.
자세한 설명은 스크립트 아래에 설명되어 있습니다.

gsidem2geotiff.sh

#!/bin/bash

#私の環境で足りなかったものを追加．
#足りないものがあれば各自追加．
#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はDLしてきたFG*zipのまま．
python fgddemImporter/fgddem.py FG-GML-*zip

#geotiffファイルをmosaic (merge.tif)．海領域の値を-9999から0に．
#-9999 (-srcnodata指定) の値を0 (-dtsnodata指定) にする．
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

#geotiffの中身を表示
gdalinfo merge.tif

#不要だと思いますがgeotiff->ESRI ascii形式に変換．-9999をNaNとする．
#-projwinで切り出しています．西から時計回りに端の座標 (メートル) を指定．
gdal_translate -projwin ${west} ${north} ${east} ${south} -a_nodata -9999 -of AAIGrid merge.tif dem.asc

#EOF

중요한 것은 국토 지리원에서 다운로드 한 DEM 파일을 zip 파일에서 geotiff로 변환, geotiff 모자이크, ESRI ascii 형식 (자르기)까지의 변환이 터미널에서 (이번에는 Ubuntu18.04LTS, bash 환경)에서 완결 할 수있는 일입니다.

실천예 -GSI DEM을 merge하고 일부 영역만 잘라내고 싶다-

이번에는 사쿠라지마를 예로 들어 봅시다.
절차의 습관입니다.
1. GSI HP에서 DEM 파일 다운로드
2. JPGIS (GML) 형식을 Geotiff 형식으로 변환
3. UTM 좌표로 변환
4. UTM 좌표 지정으로 잘라 내기 및 ESRI ascii 형식으로 변환

이번에는 아래와 같은 영역의 데이터를 사용합니다.


위의 스크립트에 잘라내기 영역을 지정하는 구체적인 값을 넣은 것입니다.
이 절삭 영역은 UTM 좌표에서 (west, north, east, south) = (650500 3515000 665300 3490000)로 설정되었습니다.

결과 (QGIS로 그리기)

우선 원하는대로 히로이 영역에서 GSI DEM을 가져와 원하는 영역만 잘라서 UTM 좌표로 변환할 수 있었습니다.
여기에서는 gdalwarp의 커멘드로 mosaic과 UTM 변환으로 원 쿠션 놓고 있습니다만, EQA 좌표의 mosaic한 geotiff 파일이 불필요하면 옵션을 연결하는 것으로 mosaicUTM 좌표 변환을 일발로 할 수 있습니다.

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

먼저 EQA 좌표 지정으로 잘라내기 → UTM 좌표 변환해 보자

방금전은 EQA→UTM 좌표 변환→UTM 좌표 지정 영역 잘라내기를 했습니다만, 굳이 EQA 영역에서 잘라내기→UTM 좌표 변환을 하면 조금 외형이 바뀝니다.

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

결과 (QGIS로 그리기)

위도 경도 지정 자르기 → EQA에서 UTM 좌표 변환이므로 좌표로서는 올바르게 변환되었지만, 외형으로서 반 시계 방향으로 회전한 영역이 되었습니다.
이것을 더 이용하면, 약간 사용하기 어려울지도 모릅니다.

마지막으로

쉘 스크립트로 파이썬 스크립트를 읽는 것은 사악한 느낌입니다.
파이썬만으로도 잘 할 수 있기 때문에 ...