Pythonで国土地理院の標高データを画像化してみる

前職の勤務中、上長の趣味()で作らされたPythonで国土地理院の数値標高モデルをPNG化するプログラム。
眠らせるのは勿体無いので、次の学習者に託します

言わずもがなですが、下記のコードは『最善のコード』ではありません！
また、私はPythonが専門でもないので、「誰かリファクタリングしてくれたら嬉しいなぁ」なんて思いつつ書いています

※補足(2020/02/14) このコードは、私がまだ一年目のころに書いたコードなので、出来が悪いです。お気を付けください。

やりたいこと

国土地理院の数値標高データを取得し、Pythonを使ってグレースケール画像にしたい
データファイル一つにつき一枚のPNGファイルを生成することにする

落とせるxmlファイルの説明は下記のQiita記事が凄くわかりやすいので参考にしてください

https://qiita.com/tobira-code/items/43a23362f356198adce2

というか、この記事があれば私がワザワザ新しく記事を書く必要もないような・・・

データを準備

さっそく、国土地理院のダウンロードページから、任意のメッシュを選択し取得します

https://fgd.gsi.go.jp/download/menu.php

『基盤地図情報 数値標高モデル』の「ファイル選択へ」というボタンを押すと、ダウンロード画面へ飛びます

データをダウンロードするには会員登録が必要です！
今回は、「10mメッシュ」の「地形図等高線」のデータを使用するので、左側のラジオボタンとチェックボタンをそのように選択し、任意のメッシュをポチッとクリックしてください

落としたZipファイルから、『FG-GML-0000-0-dem10b-yyyymmdd.xml』となっているファイルを解凍して適当なディレクトリに置いてください
コレでデータの準備は完了です。

コードは下記。

import os
import re
import sys
import numpy as np
from PIL import Image
import xml.etree.ElementTree as ET

# xmlファイル入力
DATA = sys.argv[1]
tree = ET.parse(DATA)
root = tree.getroot()
namespace = {
    'xml': 'http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema',
    'gml': 'http://www.opengis.net/gml/3.2'
}
dem = root.find('xml:DEM', namespace)
# メッシュ番号
mesh = dem.find('xml:mesh', namespace).text
# セルの配列数(実際の値は1ずつ追加)
high = dem.find('./xml:coverage/gml:gridDomain/gml:Grid/gml:limits/gml:GridEnvelope/gml:high', namespace).text.split(' ')
highX = int(high[0]) + 1
highY = int(high[1]) + 1

# 画像サイズ設定(データ数 == ピクセル数)
imgSize = highX * highY

# 標高データの配列
dem_text = dem.find('./xml:coverage/gml:rangeSet/gml:DataBlock/gml:tupleList', namespace).text
data = re.findall(',(.*)\n', dem_text)
dataNp = np.empty(imgSize)
for i in range(len(data)):
    if(data[i] == "-9999.00"):
        dataNp[i] = 0
    else:
        dataNp[i] = float(data[i])

# データの開始座標
startPoint = dem.find('./xml:coverage/gml:coverageFunction/gml:GridFunction/gml:startPoint', namespace).text.split(' ')
startPointX = int(startPoint[0])
startPointY = int(startPoint[1])
startPosition = startPointY * highX + startPointX

## 画像生成開始
# データ数が不足している場合(画像の上下に空白が有る場合)
if(len(dataNp) < imgSize):
    add = []
    # 画像下部のデータが不足している場合
    if(startPosition == 0):
        for i in range(imgSize - len(dataNp)):
            add.append(0)
        dataNp.extend(add)
    # 画像上部及び下部のデータが不足している場合
    else:
        for i in range(startPosition):
            add.append(0)
        dataNp[0:0] = add
        add = []
        for i in range(imgSize - len(dataNp) - len(add)):
            add.append(0)
        dataNp.extend(add)

# データの8bit整数化
dataNp = (dataNp / 15).astype(np.uint8)  # 富士山最高所を255に収める場合、dataNpを15で割る
data = dataNp.reshape(highY, highX)

# 数値標高データを画像化
pilImg = Image.fromarray(np.uint8(data))
pilImg = pilImg.resize((int(highX), int(highY)), Image.LANCZOS) # NEAREST

canvas = Image.new('RGB', (highX, highY), (0, 0, 0))
canvas.paste(pilImg, (0, 0))
canvas.save('dem.png', 'PNG', quality=100)

実行は、python .pyファイル .xmlファイル とすればOKです。

補足

各種データの検索

最初のメッシュ番号やセルの配列数・標高データ等の取得は上記の書き方でも良いですし、 こちらのような 感じでループを回して正規表現で探してもいいですね。

ワタシ的にははループを回して探す方が好きだし楽なのでオススメなのですが、それではこの記事と被っちゃうので・・・

後半に出てくる「dataNp / 15」

ゴリゴリのマジックナンバーで恐縮なのですが、なんだか面倒くさいのでそのままにしています。

数値標高モデルでは、水面を除いた各地点の標高がそのまま数値として保存されています(全てのDEMがそうかは分かりませんが)。
富士山頂の地点は3776.0とかでxmlファイルに格納されているはずです。

その値を直接グレースケールのPNG画像にしてしまうと、白黒画像は0~255の値で表現されるので、高度255m以上の土地は全て真っ白になってしまいます。
それを防ぐために、上記コードでは富士山頂の標高3770mを基準として255以下にするために、あえて15で割って調整しているわけです。

なので、別に富士山頂じゃなくて、例えば山頂の標高が333mの函館山を0~255で表現する場合は、1.3で割ってやればいい感じになるはずですね。

ここは、xmlファイルの最高値を調べて自動で割る値を設定させるようにしたら楽そうですね。

最後に

こんな感じで、函館をPNG化したらこんな感じになりました。

市街地の標高が殆ど0で、函館山のいきなり標高上がる感がよく出ています。
ただ、イマイチパッとしないような？
なので、次回は複数のファイルを読み込ませてメッシュ番号ごとに並び替え、デカイ一枚絵にするという記事を書こうと思います。
こうすると、なかなかおもしろい絵が出来ます。
この記事も近いうちに見直して、もっと分かりやすく書き直さねばとは思っています。。。(多分やらない)