はじめに

国土地理院が公開している数値標高モデル(DEM)の利用方法を記載します。
DEMデータを2次元標高データに変換する方法を記載します。

ファイル名

基盤地図情報（数値標高モデル）データ 5mメッシュ（標高）
のファイル名は次の形式になっています。

フォーマット

FG-GML-pppp-qq-rr-DEM5X-yyyymmdd.xml

ppppは1次メッシュ番号
qqは2次メッシュ番号
rrは3次メッシュ番号
Xは測量方法で、A:レーザー測量、B:写真測量です。
yyyymmddは年月日です。

FG-GML-5338-06-89-DEM5A-20130702.xml

メッシュ番号

1次メッシュは幅が経度1度、高さが緯度40分の四角形になります。
1次メッシュ番号はこの四角形の南西端の緯度経度で表します。
1次メッシュ番号は4桁で表します。上位2桁が緯度の1.5倍になります。下位2桁は経度の下2桁になります。

2次メッシュは1次メッシュを縦横に8等分したものになります。
南西端を(0,0)として(南北,東西)に番号付します。番号の範囲は0~7になります。
南北、東西の順に並べた2桁の数字が2次メッシュ番号になります。

3次メッシュは2次メッシュを縦横に10等分したものになります。
番号付は2次メッシュと同じ考え方になります。

ファイルの内容

主な情報のみ抜粋します

mesh番号

1次,2次,3次メッシュ番号を並べたものです。

  <mesh>53380689</mesh>

GridEnvelope

メッシュを構成するグリッドセルの配列数を記載します。
low属性は北西端を示す番号です。常に(0,0)です。
high属性は南東端を示す番号です。5mメッシュは(224,149)、10mメッシュは(1124,749)固定です。

  <gml:limits>
      <gml:GridEnvelope>
          <gml:low>0 0</gml:low>
          <gml:high>224 149</gml:high>
      </gml:GridEnvelope>
  </gml:limits>

DataBlock

各グリッドセルの標高データが記載されます。
tupleListタグで囲まれたデータが標高データになります。
水面や標高データがない場合は、標高の値が"-9999."になります。

  <gml:DataBlock>
      <gml:rangeParameters>
          <gml:QuantityList  uom="DEM構成点"/>
      </gml:rangeParameters>
  <gml:tupleList>
  (標高データ)　例) 地表面,1055.33
  </gml:tupleList>
      </gml:DataBlock>

coverageFunction

              <gml:sequenceRule  order="+x-y">Linear</gml:sequenceRule>

sequenceRuleでデータの並びを指定します。
orderの"+x"はxが正の方向(西→東)であることを、"-y"はyは負の方向(北→南)であることを示します。
北西端から南東端に向けて次のような並びで一列にしたデータがtupleListタグ内に配置されます。

  <gml:startPoint>173 0</gml:startPoint>

startPointでグリッドセルの開始点を指定します。
DataBlockに配置されるデータはグリッドセル番号の(173,0)~(224,149)になります。
実際のデータ数は必ずしも(gml:low+1)*(gml:high+1)にならないことに注意してください。

データの取り出し

XMLの内容が分かったのでDEMデータを2次元標高データとして取り出します。
Pythonを使って取り出すことにします。Numpyの2次元ndarray型にしています。
配列サイズは列数(gml:low+1)、行数(gml:high+1)固定にしています。
startPointより前のデータはすべて"-9999."にしています。

hyoukou.py

import sys
import re
import numpy as np
import matplotlib.pyplot as plt

with open(sys.argv[1], "r") as f:
    # mesh
    r = re.compile("<mesh>(.+)</mesh>")
    for ln in f:
        m = r.search(ln)
        if m != None:
            mesh = m.group(1)
            break
    # grid len
    r = re.compile("<gml:high>(.+) (.+)</gml:high>")
    for ln in f:
        m = r.search(ln)
        if m != None:
            xlen = int(m.group(1)) + 1
            ylen = int(m.group(2)) + 1
            break
    # start
    r = re.compile("<gml:tupleList>")
    for ln in f:
        m = r.search(ln)
        if m != None:
            break
    # data
    data = []
    r  = re.compile("</gml:tupleList>")
    r2 = re.compile(",(.+)")
    for ln in f:
        m = r.search(ln)
        if m != None:
            break
        else:
            m = r2.search(ln)
            if m != None:
                data.append(float(m.group(1)))
    # start point
    startx = starty = 0
    r = re.compile("<gml:startPoint>(.+) (.+)</gml:startPoint>")
    for ln in f:
        m = r.search(ln)
        if m != None:
            startx = int(m.group(1))
            starty = int(m.group(2))
            break

data2 = np.empty(xlen*ylen)
start_pos = starty*xlen + startx

for i in range(xlen*ylen):
    if i < start_pos:
        data2[i] = -9999.
    else:
        data2[i] = data[i-start_pos]

data = data2.reshape(ylen, xlen)

#print(data)
plt.imshow(data)
plt.colorbar()
plt.show()

hyoukou.pyの呼び出し方
DEMファイル名を引数にして呼び出します。
2次元標高データをmatplotlibを使ってimageにして表示します。
データ欠損(-9999.)がなければ正常に表示できます。

>python hyoukou.py FG-GML-5338-06-89-DEM5B-20130702.xml

基盤地図情報～数値標高モデルを利用する

はじめに

ファイル名

メッシュ番号

ファイルの内容

データの取り出し