Pythonで複数の数値標高データを一枚絵にする

前書いた記事の続き？のようなもの
『Pythonで国土地理院の標高データを画像化してみる』
この記事では、複数のxmlファイルをポポンと放り込んで、適切に並べて一枚絵にしたいと思います

準備

今回は、任意の都府県のメッシュデータを、まるごと落とします。

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

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

『基盤地図情報 数値標高モデル』の「ファイル選択へ」というボタンを押すと、ダウンロード画面へ飛びます。
今回は、大阪府のデータを使用することにします。
左側にある10mメッシュの10(地形図の等高線)にチェックを入れておいてください

この記事に書いてあるコードでは、国土地理院で提供されているデータをそのままグレースケール化しています
つまり、10mメッシュのデータの場合、xmlファイル一つ毎に1125x750のグレースケール画像が生成され、それを一つ一つ並べて一枚の大きな画像を生成することになります。
大阪府の場合は、横は6x1125px,縦は10x750pxの馬鹿デカイ画像が生成されることになるので、選択する都府県が小さいに越したことはありません。
もし一枚絵の画素数を小さくしたい場合は、後述の箇所を適当に弄ってください。

さて、落としてきたZIPファイルは一度解答すると、以下のようなデータ構造になっているかと思います。

解凍したディレクトリ
|-- FG-GML-0000-00-DEM10B.zip
|-- FG-GML-0000-00-DEM10B.zip
|-- etc...

それぞれのZIPファイルを更に解凍しても良いのですが、少し面倒くさいので今回はZIPファイルのまま扱います

コード全体は、下記の感じ

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

WIDTH = 1125
HEIGHT = 750

# 準備
def set_value(dir, rate):
    global WIDTH
    global HEIGHT
    # ファイル名からメッシュをリスト化し並び替える
    meshList = []
    for f in glob.glob(os.path.join(dir, "*.zip")):
        base = os.path.basename(f)
        meshList.append(base[7:11] + base[12:14])
    meshList.sort()

    # 処理予定データ数
    denominator = len(meshList)
    # 処理済みデータ数
    numerator = 0

    # メッシュリストから最東西南北端のメッシュを調べる
    top = right = meshList[-1]
    bottom = left = meshList[0]
    for meshNumber in meshList:
        str(meshNumber)
        if top[0:2] < meshNumber[0:2] or (top[0:2] <= meshNumber[0:2] and top[4] <= meshNumber[4] and top[4] <= meshNumber[4]):
            top = meshNumber
        if bottom[0:2] > meshNumber[0:2] or (bottom[0:2] >= meshNumber[0:2] and bottom[4] >= meshNumber[4]):
            bottom = meshNumber
        if right[2:4] < meshNumber[2:4] or (right[2:4] <= meshNumber[2:4] and right[5] <= meshNumber[5]):
            right = meshNumber
        if left[2:4] > meshNumber[2:4] or (left[2:4] >= meshNumber[2:4] and left[5] >= meshNumber[5]):
            left = meshNumber
    # 最端のメッシュから、キャンバスサイズを出す
    vartical = (int(top[0:2])-int(bottom[0:2])) * 8 + int(top[4]) - int(bottom[4]) + 1
    horizon = (int(right[2:4])-int(left[2:4])) * 8 + int(right[5]) - int(left[5]) + 1
    maxArea = vartical * horizon
    point = top[0:2] + left[2:4]  + top[4] + left[5]
    # canvas生成
    canvas = Image.new('RGB', (int(WIDTH / rate) * horizon, int(HEIGHT / rate) * vartical), (0, 0, 0))

    # ファイル読み取りと画像化
    for zipfile in glob.glob(os.path.join(dir, "*.zip")):
        paste_image(zipfile, point, canvas, rate)
        numerator += 1
        print('%d / %d' % (numerator, denominator))
    canvas.save('dem.png', 'PNG', quality=100)


# 画像化メソッド
def paste_image(z, p, c, r):
    global WIDTH
    global HEIGHT
    point = p
    canvas = c
    rate = r
    with zipfile.ZipFile(z, "r") as zf:
        for info in zf.infolist():
            inner = info
        with zf.open(inner) as zfile:
            root = ET.fromstring(zfile.read())
            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

            # 画像サイズ設定(データ数 == ピクセル数)
            dataSize = 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(dataSize)
            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) < dataSize):
                add = []
                # 画像下部のデータが不足している場合
                if(startPosition == 0):
                    for i in range(dataSize - 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(dataSize - 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 / rate), int(highY / rate)), Image.LANCZOS) # NEAREST

            # 貼り付けるメッシュの座標を計算
            targetX = (int(point[0:2]) - int(mesh[0:2])) * 8 + int(point[4]) - int(mesh[4])
            targetY = (int(mesh[2:4]) - int(point[2:4])) * 8 + int(mesh[5]) - int(point[5])

            canvas.paste(pilImg, (targetY * int(WIDTH / rate), targetX * int(HEIGHT / rate)))

# ZIPファイルの入ったディレクトリを入力
DIR = sys.argv[1]
RATE = int(sys.argv[2])

set_value(DIR, RATE)

pythonの実行時に、ZIPファイルの詰まったディレクトリを指定すれば動くはずです

補足

set_valueメソッドで、canvasを生成する箇所

canvas = Image.new('RGB', (1125 * horizon, 750 * vartical), (0, 0, 0))

ココの1125, 750はそれぞれ数値標高データのセル数を手動で書いてあるだけです
コード上に出てくるhighX, highYと値自体は同じなので、二度手間だしhighX,Yは毎画像生成の度に値を探してたりもう色々クソなコードになっています。
(また今度書き直します、多分)
誰かリファクタリングしてくれたら嬉しいなぁ、なんて思いつつ・・・
それと、

dataNp = (dataNp / 15).astype(np.uint8)

の15は何ぞ？の答えは、前回の記事に書いてますのでご参照願います。
https://qiita.com/artistan/items/99266407702d4152e9d5

デカイ画像はリサイズしよう

最初の方でも書きましたが、このままのコードではデッカイ画像が生成されて不便です。
なので、set_valueメソッドの下記箇所を

canvas = Image.new('RGB', (1125 * horizon, 750 * vartical), (0, 0, 0))

から

canvas = Image.new('RGB', (1125/10 * horizon, 750/10 * vartical), (0, 0, 0))

にして、コメントアウトしているpilImageをresizeする箇所を

pilImg = pilImg.resize((int(highX)/ 10, int(highY) / 10), Image.LANCZOS) # NEAREST

とかにすれば、10分の1に縮小された画像が生成されて便利よねー、という感じ

最後に

上記コードで生成し、適当に縮小した大阪府がこちらになります

大阪はそんなに山の高さがないので面白くないですが、もっと起伏にとんだ場所なら中々見ごたえがあると思います
グレースケールだと山の形が葉脈みたいになっていて、なんだか新鮮ですね
分かりにくい説明しか出来ませんでしたが、読んでいただきありがとうございました