生田緑地、石垣一夜城、大森貝塚、八王子城、多摩川台公園と等々力渓谷なんてな S.W. の出歩きは、言ってみれば「独特な地形巡り」でもあって。
たぶん一ヶ月ぶりくらいの「技術ネタ」なんだろうけど、自分的には丸っきりシルバーウィーク期間の出歩きと地続きなのよね。そ、遊びです、遊び。
実際に行って目で見た「楽しい地形を再現してみたい」というのも自然発生的にあったりもするけれど、だいたいこういう史跡とかがあるとこって、立体地形図なんかが置いてあるわけね。八王子城のガイダンス施設もそう。多摩川台公園の古墳資料室もそう。みるたんびに、欲しいなこれ、と思うのよね。丸の内丸善なんかで立体地形図売ってたりすると、欲しくて仕方なくなったりもするしさ。ただでさえ地図が好きなのに、立体となると嬉しくて仕方なくなっちゃってなぁ。
色々やり方は考えられて、たぶん地理院地図(電子国土WEB)の「3D」機能で 3D プリンタ用のデータを取得して、kinko’sで出力、とかがイマドキなんだろうと思う。これはこれでやってみたいとは思うんだけれど、まずは「図工」したいし、それ以前の問題として、ちょっと自由度に欠ける気がすんのよね、今の目的には。スケールに相当の自由度が欲しいわけですよ。富士山とかの山であれば、20万図レベルの広域が欲しいだろうし等高線は 100m 解像度くらいが良いだろうし、八王子城なんかでは1万か2万5000図レベルの、等高線10m刻みくらいは欲しそうだし、田園調布付近の拡大では等高線は 1m 解像度が欲しいだろうし。
なんて色々考えてたんだけど、なんかの拍子に「弁当パックで作る立体地形図」を見かけて、なるほど、と思ったのね。
ダンボールを重ねるとかは考えたけど、切る手間と色塗りの手間まで考えるとちょっと果てしないなと思う。でも「透明の弁当パック」で作るこのやり方は、一日以内の作業で作れそうだものね。なお「透明の弁当パック」については、要は「各層底面が浮いた状態で重ねられる透明なもの」であれば良いので、100円ローソンとかでまとめ売りしてる透明のクリアケースをちょっと加工(端を斜めに折ればいいだろう)しても作れるんじゃないかなと思ってる。こっちの方が安い可能性ある。いくらしたっけな?
そうかそうか、それで作ってみよう、と思い、「図工」部分の想像を思い描くことが出来たので、じゃぁ、自分が作りたい場所の等高線図を作ろう、と。
無論2万5千図とかを書店で買ってきてそれをトレースする、みたいな原始的なやり方も面白いだろうけど、「超イマドキではなくてもややイマドキ」は基盤地図情報の標高データが無償で使える(要ユーザ登録)ので、これを使って等高線図を作っちゃろうと。
狙い撃ちで欲しい場所のデータを取得したかったのだけど、どうやら「重い」らしく、地図から選択の UI が死んでる。というか待てども暮らせども地図が出ない。仕方ないので東京・神奈川の２県分を「まとめてダウンロード」した。が、ダウンロードしたファイルが壊れてると、開けない。なんだろう…。仕方ないので全部チマチマと一つ一つダウンロードした。地図から選択しないもんだから、どこがどこだかわからんので、「全部」取るしかなくて参った。400メガくらいあったかしら。
ひとまずいくつかエクスプローラから開いてみて、ダウンロードした zip の中にどのくらいの GML が含まれてるかをさっと眺めて、「東京・神奈川全ての GML を一個のフォルダに放り込んだらエラいことになる」ことがわかったので、zip の単位でフォルダを作って、そこに展開しておく：
1 me@host: ~$ for i in __arch/*.zip ; do d=__extracted/`basename $i .zip` ; mkdir $d ; unzip.exe $i -d $d ; done
ワタシは仕事関係で海岸線と河川の基盤地図情報を使ったことがあって、基盤地図の GML が「馬鹿げていて扱いにくい」ことは重々承知している。「XMLSchema の UML が美しいので仕事した気分になれる」効果しかない GML 3 であるばかりでなく、「いまさらの Shift-JIS」なことが、輪をかけて扱いにくい、のね。だから最初からマジメに XML として扱おうなんて思ってない。と、「不真面目に XML る」ことだけは心に決めている上で、標高データは初めてなので、まずは構造を眺めておく。げっ…：
1 <?xml version="1.0" encoding="Shift_JIS"?>
2
3 <Dataset xsi:schemaLocation="http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema FGD_GMLSchema.xsd"
4 xmlns:gml="http://www.opengis.net/gml/3.2"
5 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
6 xmlns:xlink="http://www.w3.org/1999/xlink"
7 xmlns="http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema"
8 gml:id="Dataset1">
9 <gml:description>基盤地図情報メタデータ ID=fmdid:13-3101</gml:description>
10 <gml:name>基盤地図情報ダウンロードデータ（GML版）</gml:name>
11
12 <DEM gml:id="DEM001">
13 <fid>fgoid:10-00100-13-60101-48395630</fid>
14 <lfSpanFr gml:id="DEM001-1">
15 <gml:timePosition>2013-07-02</gml:timePosition>
16 </lfSpanFr>
17 <devDate gml:id="DEM001-2">
18 <gml:timePosition>2013-07-02</gml:timePosition>
19 </devDate>
20 <orgGILvl>0</orgGILvl>
21 <orgMDId>H24GC009</orgMDId>
22 <type>5mメッシュ（標高）</type>
23 <mesh>48395630</mesh>
24 <coverage gml:id="DEM001-3">
25 <gml:boundedBy>
26 <gml:Envelope srsName="fguuid:jgd2011.bl">
27 <gml:lowerCorner>32.441666667 139.75</gml:lowerCorner>
28 <gml:upperCorner>32.45 139.7625</gml:upperCorner>
29 </gml:Envelope>
30 </gml:boundedBy>
31 <gml:gridDomain>
32 <gml:Grid dimension="2" gml:id="DEM001-4">
33 <gml:limits>
34 <gml:GridEnvelope>
35 <gml:low>0 0</gml:low>
36 <gml:high>224 149</gml:high>
37 </gml:GridEnvelope>
38 </gml:limits>
39 <gml:axisLabels>x y</gml:axisLabels>
40 </gml:Grid>
41 </gml:gridDomain>
42 <gml:rangeSet>
43 <gml:DataBlock>
44 <gml:rangeParameters>
45 <gml:QuantityList uom="DEM構成点"/>
46 </gml:rangeParameters>
47
48 <gml:tupleList>
49 海水面,-9999.
50 海水面,-9999.
51 ...(省略)...
52 データなし,-9999.
53 データなし,-9999.
54 データなし,-9999.
55 海水面,-9999.
56 海水面,-9999.
57 海水面,-9999.
58 海水面,-9999.
59 </gml:tupleList>
60
61 </gml:DataBlock>
62 </gml:rangeSet>
63 <gml:coverageFunction>
64 <gml:GridFunction>
65 <gml:sequenceRule order="+x-y">Linear</gml:sequenceRule>
66 <gml:startPoint>71 0</gml:startPoint>
67 </gml:GridFunction>
68 </gml:coverageFunction>
69 </coverage>
70 </DEM>
71 </Dataset>
うげぇ。なんつぅ構造じゃ…。
tupleList 部分の格納イメージが想像出来ない。startPoint はひょっとすると、GridEnvelope つまり配列の shape は全ファイル共通で、先頭に無効データ続く場合にスキップするための情報、か? と最初に想像し、いくつかのファイルの tupleList リストの行数とそれらメタ情報との関係を探ってみる。GridEnvelope はいくつかみた感じでは、同じっぽい。どれも high が「224 149」つまり要素数が「250×150 = 37500」。そして tupleList リストの行数はこれと一致していたりしていなかったりする。しかも startPoint ぶん引いてみてもまだ一致しない。てことは、末尾に無効データ続く場合省略? 要素数問題もそうだが、東西南北の関係は?
と、およその想像出来る部分と出来ない部分がありつつも、等高線を描くのには Matplotlib を使おうと思っているし、そうであれば NumPy の ndarray にブチ込めばいいわけなので、であれば 1次元配列に詰め込んでから reshape するなりすれば良かろう、と。startPoint の扱いだって、ndarray の範囲コピーで一撃で解決するであろうし、ということで、完璧な理解をせぬまま作業をはじめた。
結局3回くらいの試行錯誤だけで出来た：
1 # -*- coding: utf-8 -*-
2 import re
3 import numpy as np
4
5 _RGXES = (
6 (re.compile(r"<gml:lowerCorner>(.*) (.*)</gml:lowerCorner>"), "lc", float),
7 (re.compile(r"<gml:upperCorner>(.*) (.*)</gml:upperCorner>"), "uc", float),
8 (re.compile(r"<gml:low>(.*) (.*)</gml:low>"), "L", int),
9 (re.compile(r"<gml:high>(.*) (.*)</gml:high>"), "H", int),
10 (re.compile(r"<gml:startPoint>(.*) (.*)</gml:startPoint>"), "st", int),
11 )
12
13 srcgml = "__extracted/FG-GML-5239-03-DEM5A/FG-GML-5239-03-43-DEM5A-20130702.xml"
14 #srcgml = "__extracted/FG-GML-4839-56-DEM5A/FG-GML-4839-56-30-DEM5A-20130702.xml"
15
16 _grid_info = {}
17 _elevs_raw = []
18
19 with open(srcgml, "r") as fi:
20 # grid infos
21 ri = 0
22 while True:
23 line = fi.readline()
24 m = _RGXES[ri][0].search(line)
25 if m:
26 _grid_info[_RGXES[ri][1]] = (
27 _RGXES[ri][2](m.group(1)), _RGXES[ri][2](m.group(2)))
28 ri += 1
29 if len(_RGXES) == ri + 1:
30 break
31
32 # find start of grid data
33 while True:
34 line = fi.readline()
35 if "gml:tupleList" in line:
36 break
37
38 # read grid data
39 while True:
40 line = fi.readline()
41 if "gml:tupleList" in line:
42 break
43 # we need to treat sea level always NaN rather than zero.
44 # (modified at 2015-10-06)
45 ty, d = line.strip().split(",")
46 if ty.decode('cp932') in (u'データなし', u'海水面'):
47 d = '-9999.'
48 # invalid data are stored with '-9999.0'
49 # (modified at 2015-10-05)
50 if float(d) == -9999.0:
51 #if d == u'-9999.': # No good.
52 d = u'NaN'
53 _elevs_raw.append(float(d))
54
55 # read start-point
56 while True:
57 line = fi.readline()
58 m = _RGXES[ri][0].search(line)
59 if m:
60 _grid_info[_RGXES[ri][1]] = (
61 _RGXES[ri][2](m.group(1)), _RGXES[ri][2](m.group(2)))
62 break
63
64 #print(_grid_info, _elevs_raw)
65 #print(_grid_info)
66 shape = (_grid_info["H"][1] + 1, _grid_info["H"][0] + 1)
67
68 st = _grid_info["st"][1] * shape[1] + _grid_info["st"][0]
69 elevs = np.zeros(shape[1] * shape[0])
70 elevs[st:len(_elevs_raw) + st] = _elevs_raw
71 elevs = elevs.reshape((shape[0], shape[1]))
72 elevs = np.flipud(elevs)
73 print(elevs)
74
75 #
76 import matplotlib
77 import matplotlib.cm as cm
78 import matplotlib.mlab as mlab
79 import matplotlib.pyplot as plt
80 from matplotlib.ticker import FormatStrFormatter
81
82 #fig, ax = plt.figure()
83 fig, ax = plt.subplots()
84 ax.xaxis.set_major_formatter(FormatStrFormatter('%.4f'))
85 ax.yaxis.set_major_formatter(FormatStrFormatter('%.4f'))
86
87 X = np.linspace(_grid_info["lc"][1], _grid_info["uc"][1], shape[1])
88 Y = np.linspace(_grid_info["lc"][0], _grid_info["uc"][0], shape[0])
89 CS = plt.contour(X, Y, elevs)
90 plt.clabel(CS, inline=1, fontsize=10)
91 plt.show()
後半の作図部分は「3回」にカウントしてない。「たぶんほぼ正解」になるまでは経度(X)緯度(Y)とのマッピングも set_major_formatter もやってない。試行錯誤したのはハイライトした 59～64 行目。上下転地(flipud)は、経験的に「きっと (0, 0) が北西端だろ」と、勘で。
「たぶんほぼ正解」になったあとで FAQ をみたら、必要なことは全部書いてあった：
てわけで「合ってんじゃん」、すなわち、「たぶんほぼ正解」は「たぶん正解」とわかった。
そうなんだけど、「日本地図の形がわかる」ようなスケールじゃないからね、それっぽい絵がかけても「正解」感がよくわからんのね。試したファイルも、いったい東京神奈川のどこ指してんのかわからずにやってるし。ので、地理院地図(電子国土WEB)と比較を試みる。あ、これ、三原山(大島)か…。てわけで、ちゃんと正解になってた：
よしよし。第一段階としてはこれでいいな。
ちなみに「ユーザ登録・ダウンロード」は前日に済ませてあったので、本日作業したのは1時間かかってない。作図の仕上げを除いた「たぶんほぼ正解」部分は20分くらいで書いてる。(つまりこの投稿を書いてる時間の方が長いのよ…。)
あとはこの大量のファイルから素早く欲しい範囲を取り出すためのシカケをどうにかして、あと広範囲を扱いたい場合の array の結合とか、高解像過ぎるので平滑してデータセットを小さくするとか、等高線間隔を自由に制御するとか…、と、まだ結構先は長いが…、ま、のんびりやりましょ。
2015-10-05追記
NaN にする -9999.0 判定がまずかった。どうせテンプレートで作ってるんでしょと思い込み、どの -9999.0 も「-9999.」なのだとしてしまってたが、ダメ。やぁねぇ…。
2015-10-06追記
海水面データは -9999.0 だ、と思っていたら、ゼロのがいた。横浜あたり。困る…。ここが統一されてないと、海が海にみえなくなってしまう。
