これの、のんびりした続き。
前回、基盤地図情報標高データ(5mメッシュ)の東京・神奈川全部をダウンロードしたわけだ。けど「地図的な」インターフェイスで持ってきたわけではないし、ファイル名からはどこを指すのか見当もつかず、開いてみるしかないわけね。まずは欲しい場所から簡単に取り出せるようにせねば。
どんだけあんだ、と、まずファイル数を調べておく：
MSYS bash から。遅い。MSYS は I/O が全般に遅いのよ。
1 me@host: FG-GML$ find __extracted/ -name '*.xml' | wc -l
2 5913
3 me@host: FG-GML$ ls -1 __extracted/FG-GML-5339-40-DEM5A/
4 FG-GML-5339-40-00-DEM5A-20130702.xml
5 FG-GML-5339-40-01-DEM5A-20130702.xml
6 FG-GML-5339-40-02-DEM5A-20130702.xml
7 FG-GML-5339-40-03-DEM5A-20130702.xml
8 FG-GML-5339-40-04-DEM5A-20130702.xml
9 FG-GML-5339-40-05-DEM5A-20130702.xml
10 FG-GML-5339-40-06-DEM5A-20130702.xml
11 FG-GML-5339-40-07-DEM5A-20130702.xml
12 FG-GML-5339-40-08-DEM5A-20130702.xml
13 FG-GML-5339-40-09-DEM5A-20130702.xml
14 FG-GML-5339-40-10-DEM5A-20130702.xml
15 FG-GML-5339-40-11-DEM5A-20130702.xml
16 FG-GML-5339-40-12-DEM5A-20130702.xml
17 FG-GML-5339-40-13-DEM5A-20130702.xml
18 FG-GML-5339-40-14-DEM5A-20130702.xml
19 FG-GML-5339-40-15-DEM5A-20130702.xml
20 FG-GML-5339-40-16-DEM5A-20130702.xml
21 FG-GML-5339-40-17-DEM5A-20130702.xml
22 FG-GML-5339-40-18-DEM5A-20130702.xml
23 FG-GML-5339-40-19-DEM5A-20130702.xml
24 FG-GML-5339-40-20-DEM5A-20130702.xml
25 FG-GML-5339-40-21-DEM5A-20130702.xml
26 FG-GML-5339-40-22-DEM5A-20130702.xml
27 FG-GML-5339-40-23-DEM5A-20130702.xml
28 FG-GML-5339-40-24-DEM5A-20130702.xml
29 FG-GML-5339-40-25-DEM5A-20130702.xml
30 FG-GML-5339-40-26-DEM5A-20130702.xml
31 FG-GML-5339-40-27-DEM5A-20130702.xml
32 FG-GML-5339-40-28-DEM5A-20130702.xml
33 FG-GML-5339-40-29-DEM5A-20130702.xml
34 FG-GML-5339-40-30-DEM5A-20130702.xml
35 FG-GML-5339-40-31-DEM5A-20130702.xml
36 FG-GML-5339-40-32-DEM5A-20130702.xml
37 FG-GML-5339-40-33-DEM5A-20130702.xml
38 FG-GML-5339-40-34-DEM5A-20130702.xml
39 FG-GML-5339-40-35-DEM5A-20130702.xml
40 FG-GML-5339-40-36-DEM5A-20130702.xml
41 FG-GML-5339-40-37-DEM5A-20130702.xml
42 FG-GML-5339-40-38-DEM5A-20130702.xml
43 FG-GML-5339-40-39-DEM5A-20130702.xml
44 FG-GML-5339-40-40-DEM5A-20130702.xml
45 FG-GML-5339-40-41-DEM5A-20130702.xml
46 FG-GML-5339-40-42-DEM5A-20130702.xml
47 FG-GML-5339-40-43-DEM5A-20130702.xml
48 FG-GML-5339-40-44-DEM5A-20130702.xml
49 FG-GML-5339-40-45-DEM5A-20130702.xml
50 FG-GML-5339-40-46-DEM5A-20130702.xml
51 FG-GML-5339-40-47-DEM5A-20130702.xml
52 FG-GML-5339-40-48-DEM5A-20130702.xml
53 FG-GML-5339-40-49-DEM5A-20130702.xml
54 FG-GML-5339-40-50-DEM5A-20130702.xml
55 FG-GML-5339-40-51-DEM5A-20130702.xml
56 FG-GML-5339-40-52-DEM5A-20130702.xml
57 FG-GML-5339-40-53-DEM5A-20130702.xml
58 FG-GML-5339-40-54-DEM5A-20130702.xml
59 FG-GML-5339-40-55-DEM5A-20130702.xml
60 FG-GML-5339-40-56-DEM5A-20130702.xml
61 FG-GML-5339-40-57-DEM5A-20130702.xml
62 FG-GML-5339-40-58-DEM5A-20130702.xml
63 FG-GML-5339-40-59-DEM5A-20130702.xml
64 FG-GML-5339-40-60-DEM5A-20130702.xml
65 FG-GML-5339-40-61-DEM5A-20130702.xml
66 FG-GML-5339-40-62-DEM5A-20130702.xml
67 FG-GML-5339-40-63-DEM5A-20130702.xml
68 FG-GML-5339-40-64-DEM5A-20130702.xml
69 FG-GML-5339-40-65-DEM5A-20130702.xml
70 FG-GML-5339-40-66-DEM5A-20130702.xml
71 FG-GML-5339-40-67-DEM5A-20130702.xml
72 FG-GML-5339-40-68-DEM5A-20130702.xml
73 FG-GML-5339-40-69-DEM5A-20130702.xml
74 FG-GML-5339-40-70-DEM5A-20130702.xml
75 FG-GML-5339-40-71-DEM5A-20130702.xml
76 FG-GML-5339-40-72-DEM5A-20130702.xml
77 FG-GML-5339-40-73-DEM5A-20130702.xml
78 FG-GML-5339-40-74-DEM5A-20130702.xml
79 FG-GML-5339-40-75-DEM5A-20130702.xml
80 FG-GML-5339-40-76-DEM5A-20130702.xml
81 FG-GML-5339-40-77-DEM5A-20130702.xml
82 FG-GML-5339-40-78-DEM5A-20130702.xml
83 FG-GML-5339-40-79-DEM5A-20130702.xml
84 FG-GML-5339-40-80-DEM5A-20130702.xml
85 FG-GML-5339-40-81-DEM5A-20130702.xml
86 FG-GML-5339-40-82-DEM5A-20130702.xml
87 FG-GML-5339-40-83-DEM5A-20130702.xml
88 FG-GML-5339-40-84-DEM5A-20130702.xml
89 FG-GML-5339-40-85-DEM5A-20130702.xml
90 FG-GML-5339-40-86-DEM5A-20130702.xml
91 FG-GML-5339-40-87-DEM5A-20130702.xml
92 FG-GML-5339-40-88-DEM5A-20130702.xml
93 FG-GML-5339-40-89-DEM5A-20130702.xml
94 FG-GML-5339-40-90-DEM5A-20130702.xml
95 FG-GML-5339-40-91-DEM5A-20130702.xml
96 FG-GML-5339-40-92-DEM5A-20130702.xml
97 FG-GML-5339-40-93-DEM5A-20130702.xml
98 FG-GML-5339-40-94-DEM5A-20130702.xml
99 FG-GML-5339-40-95-DEM5A-20130702.xml
100 FG-GML-5339-40-96-DEM5A-20130702.xml
101 FG-GML-5339-40-97-DEM5A-20130702.xml
102 FG-GML-5339-40-98-DEM5A-20130702.xml
103 FG-GML-5339-40-99-DEM5A-20130702.xml
104 me@host: FG-GML$ for i in __extracted/* ; do echo -n `basename $i`" " ; ls $i | wc -l ; done
105 FG-GML-4839-56-DEM5A 11
106 FG-GML-4939-46-DEM5A 23
107 FG-GML-4939-55-DEM5A 21
108 FG-GML-4939-56-DEM5A 62
109 FG-GML-5039-64-DEM5A 27
110 FG-GML-5039-65-DEM5A 5
111 FG-GML-5139-03-DEM5A 13
112 FG-GML-5139-04-DEM5A 25
113 FG-GML-5139-13-DEM5A 8
114 FG-GML-5139-14-DEM5A 32
115 FG-GML-5139-20-DEM5A 7
116 FG-GML-5139-21-DEM5A 34
117 FG-GML-5139-31-DEM5A 11
118 FG-GML-5139-32-DEM5A 5
119 FG-GML-5139-41-DEM5A 14
120 FG-GML-5139-42-DEM5A 32
121 FG-GML-5139-52-DEM5A 9
122 FG-GML-5139-62-DEM5A 11
123 FG-GML-5238-67-DEM5A 100
124 FG-GML-5238-77-DEM5A 100
125 FG-GML-5238-77-DEM5B 11
126 FG-GML-5239-02-DEM5A 14
127 FG-GML-5239-03-DEM5A 58
128 FG-GML-5239-12-DEM5A 18
129 FG-GML-5239-13-DEM5A 29
130 FG-GML-5239-50-DEM5A 84
131 FG-GML-5239-51-DEM5B 11
132 FG-GML-5239-54-DEM5A 10
133 FG-GML-5239-55-DEM5A 20
134 FG-GML-5239-55-DEM5B 1
135 FG-GML-5239-60-DEM5A 87
136 FG-GML-5239-60-DEM5B 8
137 FG-GML-5239-61-DEM5A 7
138 FG-GML-5239-61-DEM5B 16
139 FG-GML-5239-64-DEM5A 24
140 FG-GML-5239-64-DEM5B 3
141 FG-GML-5239-65-DEM5A 66
142 FG-GML-5239-65-DEM5B 6
143 FG-GML-5239-70-DEM5A 100
144 FG-GML-5239-71-DEM5A 87
145 FG-GML-5239-72-DEM5A 39
146 FG-GML-5239-72-DEM5B 5
147 FG-GML-5239-73-DEM5A 31
148 FG-GML-5239-73-DEM5B 7
149 FG-GML-5239-74-DEM5A 71
150 FG-GML-5239-74-DEM5B 4
151 FG-GML-5239-75-DEM5A 56
152 FG-GML-5239-75-DEM5B 9
153 FG-GML-5338-07-DEM5A 97
154 FG-GML-5338-07-DEM5B 17
155 FG-GML-5338-17-DEM5A 85
156 FG-GML-5338-57-DEM5A 91
157 FG-GML-5338-67-DEM5A 18
158 FG-GML-5339-00-DEM5A 92
159 FG-GML-5339-00-DEM5B 30
160 FG-GML-5339-01-DEM5A 100
161 FG-GML-5339-01-DEM5B 14
162 FG-GML-5339-02-DEM5A 100
163 FG-GML-5339-02-DEM5B 10
164 FG-GML-5339-03-DEM5A 100
165 FG-GML-5339-03-DEM5B 10
166 FG-GML-5339-04-DEM5A 100
167 FG-GML-5339-05-DEM5A 38
168 FG-GML-5339-05-DEM5B 8
169 FG-GML-5339-10-DEM5A 100
170 FG-GML-5339-11-DEM5A 100
171 FG-GML-5339-12-DEM5A 100
172 FG-GML-5339-12-DEM5B 9
173 FG-GML-5339-13-DEM5A 100
174 FG-GML-5339-13-DEM5B 10
175 FG-GML-5339-14-DEM5A 100
176 FG-GML-5339-15-DEM5A 78
177 FG-GML-5339-15-DEM5B 4
178 FG-GML-5339-16-DEM5B 75
179 FG-GML-5339-20-DEM5A 100
180 FG-GML-5339-21-DEM5A 100
181 FG-GML-5339-22-DEM5A 100
182 FG-GML-5339-23-DEM5A 100
183 FG-GML-5339-25-DEM5A 100
184 FG-GML-5339-26-DEM5A 44
185 FG-GML-5339-30-DEM5A 100
186 FG-GML-5339-31-DEM5A 100
187 FG-GML-5339-32-DEM5A 100
188 FG-GML-5339-34-DEM5A 100
189 FG-GML-5339-35-DEM5A 100
190 FG-GML-5339-36-DEM5B 13
191 FG-GML-5339-37-DEM5A 37
192 FG-GML-5339-37-DEM5B 65
193 FG-GML-5339-40-DEM5A 100
194 FG-GML-5339-41-DEM5A 100
195 FG-GML-5339-42-DEM5A 100
196 FG-GML-5339-43-DEM5A 100
197 FG-GML-5339-44-DEM5A 100
198 FG-GML-5339-45-DEM5A 100
199 FG-GML-5339-46-DEM5A 0
200 FG-GML-5339-47-DEM5A 100
201 FG-GML-5339-50-DEM5A 61
202 FG-GML-5339-51-DEM5A 72
203 FG-GML-5339-51-DEM5B 68
204 FG-GML-5339-52-DEM5A 100
205 FG-GML-5339-53-DEM5A 100
206 FG-GML-5339-54-DEM5A 100
207 FG-GML-5339-55-DEM5A 100
208 FG-GML-5339-56-DEM5A 100
209 FG-GML-5339-57-DEM5A 100
210 FG-GML-5339-60-DEM5A 9
211 FG-GML-5339-61-DEM5A 36
212 FG-GML-5339-61-DEM5B 90
213 FG-GML-5339-62-DEM5A 100
214 me@host: FG-GML$
予想通り凄い数だ。一つの zip に最大100なのか。なお、全体が「地表面」データを含まない領域は xml そのものが存在しないので、緯度経度範囲だけから計算でファイル数を求めることは出来ない。
このもとの zip の単位は維持するのが楽かな。それと緯度経度範囲から「機械的に」ファイル名を計算出来るようにしたいし。辞書やデータベース的な取り扱いよりは「計算」が速いはずなのでね。
「緯度経度範囲から「機械的に」ファイル名を計算」のためにはファイル名の規則で最初から邪魔なのがあるのよね。「DEM5A」「DEM5B」てやつね。これ、
標高モデルについて、5mメッシュは、航空レーザースキャナ測量又は写真測量の成果から作成しています。どちらの成果で作成されたかはファイル名で判断することができ、前者には”A”が後者には”B”が含まれます。
同様に、10mメッシュは火山基本図の等高線より作成したものと1/25,000地形図の等高線から作成したものがあり、ファイル名には前者には”A”が、 後者には”B”が含まれます。
てことなんだけど、単に立体地形図を作りたいだけのワタシにはまったく不要な情報であって、邪魔以外の何者でもなくて。それとな、前回も書いたけれど、GML のまんま扱うことそのものがバカらしいのね。
うーん、「zip の単位は維持」しつつ、「DEM5A」「DEM5B」とか末尾の日付情報みたいな邪魔な情報は捨てつつ、python データを pickle、という感じか、まずは。ファイル名で必要なのはメッシュ番号部分だけなのね。
「pythonデータ」についてはどうしようか。X, Y メッシュを永続化しとくとバカでかくなるので、データの持ち方はもとの GML に近いのがいいな。前回スクリプトの、numpy データにする前のものをそのまま pickle するのが良さそうだ。グリッドとメタデータは、でもファイル分けた方がいいな。メタデータだけ使うのが最初は多そうだし。緯度経度範囲からファイル名を検索するのはひとまずあとで考えよう。
てわけで、こんな：
雑なので真似しなくていいぞ
1 # -*- coding: utf-8 -*-
2 import re
3 import os
4 import pickle
5 from os import path
6 from distutils.dir_util import create_tree
7
8 _RGXES = (
9 (re.compile(r"<gml:lowerCorner>(.*) (.*)</gml:lowerCorner>"), "lc", float),
10 (re.compile(r"<gml:upperCorner>(.*) (.*)</gml:upperCorner>"), "uc", float),
11 (re.compile(r"<gml:low>(.*) (.*)</gml:low>"), "L", int),
12 (re.compile(r"<gml:high>(.*) (.*)</gml:high>"), "H", int),
13 (re.compile(r"<gml:startPoint>(.*) (.*)</gml:startPoint>"), "st", int),
14 )
15
16 def _convert_xml(srcgml):
17 _meta = {}
18 _elevs_raw = []
19
20 with open(srcgml, "r") as fi:
21 # grid infos
22 ri = 0
23 while True:
24 line = fi.readline()
25 m = _RGXES[ri][0].search(line)
26 if m:
27 _meta[_RGXES[ri][1]] = (
28 _RGXES[ri][2](m.group(1)), _RGXES[ri][2](m.group(2)))
29 ri += 1
30 if len(_RGXES) == ri + 1:
31 break
32
33 # find start of grid data
34 while True:
35 line = fi.readline()
36 if "gml:tupleList" in line:
37 break
38
39 # read grid data
40 while True:
41 line = fi.readline()
42 if "gml:tupleList" in line:
43 break
44 # we need to treat sea level always NaN rather than zero.
45 # (modified at 2015-10-06)
46 ty, d = line.strip().split(",")
47 if ty.decode('cp932') in (u'データなし', u'海水面'):
48 d = '-9999.'
49 # invalid data are stored with '-9999.0'
50 # (modified at 2015-10-05)
51 if float(d) == -9999.0:
52 #if d == u'-9999.': # No good.
53 d = u'NaN'
54 _elevs_raw.append(float(d))
55
56 # read start-point
57 while True:
58 line = fi.readline()
59 m = _RGXES[ri][0].search(line)
60 if m:
61 _meta[_RGXES[ri][1]] = (
62 _RGXES[ri][2](m.group(1)), _RGXES[ri][2](m.group(2)))
63 break
64
65 shape = (_meta["H"][1] + 1, _meta["H"][0] + 1)
66 _meta["shape"] = shape
67 _meta["st"] = _meta["st"][1] * shape[1] + _meta["st"][0]
68 del _meta["H"]
69 del _meta["L"]
70
71 # pickle meta and data
72 rgx = re.compile(r"FG-GML-(.*)-(.*)-(.*)-.*-.*.xml")
73 m = rgx.match(path.basename(srcgml))
74 dest_dir = "__converted/{}-{}".format(m.group(1), m.group(2))
75 dest_filebase = "{}".format(m.group(3))
76 create_tree(dest_dir, [dest_filebase])
77 pickle.dump(
78 _meta,
79 open(path.join(dest_dir, "{}.meta".format(dest_filebase)), "wb"))
80 pickle.dump(
81 _elevs_raw,
82 open(path.join(dest_dir, "{}.data".format(dest_filebase)), "wb"))
83
84
85 #
86 for root, dirs, files in os.walk("__extracted/"):
87 for fn in files:
88 srcgml = path.join(root, fn)
89 _convert_xml(srcgml)
90 print(srcgml)
変換後はこんな感じ：
以後は変換後のこっちしか使わない、っと。で、この新構造に従って前回と同じ場所を絵にしてみる：
1 # -*- coding: utf-8 -*-
2 import pickle
3 import numpy as np
4
5 base = "__converted/5239-03/43"
6 #base = "__converted/4839-56/30"
7
8 _m = pickle.load(open(base + ".meta", "rb"))
9 _raw = pickle.load(open(base + ".data", "rb"))
10
11 elevs = np.zeros(_m["shape"][1] * _m["shape"][0])
12 elevs[_m["st"]:len(_raw) + _m["st"]] = _raw
13 elevs = elevs.reshape((_m["shape"][0], _m["shape"][1]))
14 elevs = np.flipud(elevs)
15 print(elevs)
16
17 #
18 import matplotlib
19 import matplotlib.cm as cm
20 import matplotlib.pyplot as plt
21 from matplotlib.ticker import FormatStrFormatter
22
23 #fig, ax = plt.figure()
24 fig, ax = plt.subplots()
25 ax.xaxis.set_major_formatter(FormatStrFormatter('%.4f'))
26 ax.yaxis.set_major_formatter(FormatStrFormatter('%.4f'))
27
28 X = np.linspace(_m["lc"][1], _m["uc"][1], _m["shape"][1])
29 Y = np.linspace(_m["lc"][0], _m["uc"][0], _m["shape"][0])
30 CS = plt.contour(X, Y, elevs)
31 plt.clabel(CS, inline=1, fontsize=10)
32 plt.show()
おけ。
さてと。あとはメッシュ番号と緯度経度の対応を考えないとね。計算式が導出出来るといいんだけどな…。データなし領域もメッシュ番号を持ってるかどうか次第。雰囲気的には期待通りになってるような気がするが?
…とデータを観察してたが、あと一歩わからんなぁ、と思って、前回諦めた「地図からダウンロード」の UI を訪れてみたところ…。あ、今日は行ける：
ネットワークが不調だったのかな、前回は。
一番大きな「5339」の次の単位はこんななのね：
つまり 8×8 の 64 てことか。なるほど。で、「5339-00」の中が、東西方向10、南北方向10の100分割(00が南西端)になってるのは、地図を見ずに確認出来てたんで、あとは「5339」の単位の謎だけ解ければ良いのだけれど、ただ、「5339」を基点に、北上すると100増えて、東に行くと1増えるという関係がみて取れるので、まぁあとはどうにかなるかな。幸い「5339-00」が手許にあるし。
2015-10-05追記
NaN にする -9999.0 判定がまずかった。どうせテンプレートで作ってるんでしょと思い込み、どの -9999.0 も「-9999.」なのだとしてしまってたが、ダメ。やぁねぇ…。
2015-10-06追記
海水面データは -9999.0 だ、と思っていたら、ゼロのがいた。横浜あたり。困る…。ここが統一されてないと、海が海にみえなくなってしまう。
