GSI_DEMをUbuntuコマンドライン上だけでgeotiff変換→UTM変換→ascii変換

6302 ワード
シェルスクリプト Python GDAL DEM geotiff Python テキストリンク

国土地理院DEMをLinux (Ubuntu) のコマンドライン上のみでGISで使いやすい形式に変換したい

GISの情報ののせるためのベースとして数値標高モデル (DEM) から読み込んだ標高図や地形陰影図みたいなものがあると，GISの情報が映えて見えます．
日本の場合は国土地理院 (GSI) が国土基盤情報として5m，10mメッシュのDEMを公開しています．
GSI DEMは提供しているJPGIS (GML) 形式という形で提供されています．
Metadataのようなものでテキストエディターで開くと，ファイル上部にDEMのピクセルサイズや緯度経度などデータに関する情報が記載されています．

個人的にはgeotiff形式になっていると大変便利なのですが，geotiff形式へのフォーマット変換に関してぱっと調べると，Windows環境で変換してくれるソフトなどがヒットします．
が，Linux使いとしてはいちいち外部ソフトを立ち上げてファイル形式を変換していては，やはり面倒くさいです．

ここでは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のmosaic，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座標変換なので座標としては正しく変換されたのだが，見た目として半時計回りに回転した領域となりました．
これを更に利用するとなると，やや使いにくいかもしれません．

最後に

シェルスクリプトでpythonスクリプトを読み込むのは邪道な気がします．
pythonのみでもうまくできるはずなので...