はじめに
国土地理院が公開している数値標高モデル(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."にしています。
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