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

# はじめに
国土地理院が公開している数値標高モデル(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次メッシュと同じ考え方になります。

![image](https://qiita-image-store.s3.amazonaws.com/0/84205/3745cea0-91ec-a7b0-0118-d4edf7734f45.png)

# ファイルの内容

主な情報のみ抜粋します

- 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タグ内に配置されます。
![1.JPG](https://qiita-image-store.s3.amazonaws.com/0/84205/1fba80a3-239b-7a8d-4b8b-53110dd0b9c1.jpeg)

>		<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."にしています。

```py3: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
```