今回のは、本題ではあるけれど、ただの整理なので、「立体地形図を作りたい」という本題からすれば面白くはないと思う。
だんだん調子に乗ってきて、色んな都道府県のをかき集めだしてるんだけど、データ格納のマズさが気になってきたのね。今、こんな構造で書き出してるわけね：
(1)、(2)ではじめたときは、メタデータだけ読むことが多そうだ、と思ってメタデータとグリッドデータを分離したんだけれどもね。2次メッシュ(画像の -21)内に 10×10 の最大100、つまり meta と data 合わせて 200 って、ダル重いじゃない。
そもそも (3)で判明した通りファイル名に含まれるメッシュ番号から緯度経度範囲は特定出来る(lowerCorner, upperCornerを個々にデータ化する意味がない)し、(7)で割り切った通り shape は DEM5 で固定なので、結局「メタデータ」として必要なのは startPoint (ワタシのスクリプトでの st) だけなのよね。もはやメタデータとして分離されてる意味なんかない。「startPointだけ読み出したい」なんてニーズは想像出来ない。もうファイル一つにまとめてしまいたいな、と。
pickle そのものについても。プロトコルバージョンはデフォルトのバージョン0 を使っちゃってるけど、HIGHEST_PROTOCOLの方がやや効率がいいという話もある。サイズが小さくなるわけではないので、悩みどころではあるけれど。それと、時間効率を考えずにやってたので cPickle 使ってないけれど、Python 2.7 の場合は明示的に cPickle 使った方がいいだろうし。
で、一番気になりだしたのが、「4839-56」という、一次メッシュ-二次メッシュをフラットに置いちゃっていること。今、東京・神奈川・静岡・山梨・千葉・埼玉・宮城・富山・長野の一部を持っているんだけど、500を超えてしまった。Windows では特に「一つのフォルダ内に大量にフォルダやファイルがある」という状態が非常に不愉快なことになるので、このままデータを増やし続けると、そのうちすぐに(少なくとも精神衛生的には)もたなくなる。二次メッシュもフォルダにしてしまいたい。つまり「4839-56」の場合、「4839/56」にしてしまいたい。二次メッシュも最大で 8×8 の 64 あるのだ。
そんなわけで思い切って構造変更。まずは、GML からの変換スクリプトと converted_loader.py を書き換える前に、データ移行スクリプトを先に書いてしまって、一括で変換してしまう：
1 # -*- coding: utf-8 -*-
2 import re
3 import os
4 try:
5 # for CPython 2.7
6 import cPickle as pickle
7 except ImportError:
8 # Python 3.x or not CPython
9 import pickle
10
11 from os import path
12 from distutils.dir_util import create_tree
13
14 #
15 for root, dirs, files in os.walk("__converted"):
16 dn = path.basename(root)
17 if "-" not in dn:
18 continue
19 for fn in files:
20 if ".data" not in fn:
21 continue
22 base = path.splitext(path.join(root, fn))[0].replace(
23 path.sep, "/")
24 elevs = pickle.load(open(base + ".data", "rb"))
25 st = pickle.load(open(base + ".meta", "rb"))["st"]
26
27 # 新しい構成の場所に新しいデータを
28 newdata = {"startPoint": st, "elevations": elevs}
29 newbase = base.replace("-", "/")
30 dest_dir = path.dirname(newbase)
31 dest_filebase = path.basename(newbase)
32 create_tree(dest_dir, [dest_filebase])
33 pickle.dump(
34 newdata,
35 open(newbase + ".data", "wb"),
36 protocol=-1)
37
38 # 元のファイルを消す
39 os.remove(base + ".data")
40 os.remove(base + ".meta")
41
42 # 元のフォルダを、消せれば消す
43 try:
44 os.rmdir(path.dirname(base))
45 except:
46 pass
stとelevsだけなのでタプルでもいいちゃぁいいんだけれど、ほかにも何か入れたくなった場合のことも一応考えて、辞書にした。プロトコルバージョンは、迷ったけど HIGHEST_PROTOCOLに。構造はめでたくこうなった：
あとは GML からの変換スクリプトと converted_loader.py を改造する。前者は converter.py という名前にしてる。
まずは、データ移行がうまくいったことをとっとと確認したいので、converted_loader.py から：
1 # -*- coding: utf-8 -*-
2 try:
3 # for CPython 2.7
4 import cPickle as pickle
5 except ImportError:
6 # Python 3.x or not CPython
7 import pickle
8 import numpy as np
9
10
11 def _calc_pos0(lon_or_lat, lseg):
12 a = (lon_or_lat % 100) / lseg
13 b = (a - int(a)) / (1 / 8.)
14 c = (b / (1 / 10.)) % 10
15 return (
16 int(a),
17 int(b),
18 int(c),
19 round((int(a) + int(b) / 8. + int(c) / 8. / 10.) * lseg, 9),
20 round((int(a) + int(b) / 8. + (int(c) + 1) / 8. / 10.) * lseg, 9))
21
22
23 def _calc_from_latlon(lat, lon):
24 _fromlat = _calc_pos0(lat, 2/3.)
25 _fromlon = _calc_pos0(lon, 1.)
26 return (
27 "__converted/{}{}/{}{}/{}{}".format(
28 _fromlat[0], _fromlon[0],
29 _fromlat[1], _fromlon[1],
30 _fromlat[2], _fromlon[2]),
31 (_fromlat[3], _fromlon[3] + 100),
32 (_fromlat[4], _fromlon[4] + 100))
33
34
35 def load(lat, lon):
36 base, lc, uc = _calc_from_latlon(lat, lon)
37 shape = (150, 225)
38 elevs = np.full((shape[1] * shape[0], ), np.nan)
39
40 try:
41 _d = pickle.load(open(base + ".data", "rb"))
42 _st = _d["startPoint"]
43 _raw = _d["elevations"]
44 elevs[_st:len(_raw) + _st] = _raw
45 except IOError as e:
46 print(str(e))
47
48 elevs = elevs.reshape((shape[0], shape[1]))
49 elevs = np.flipud(elevs)
50 return elevs, lc, uc
51
52
53 def _load_range0(lc_lat, lc_lon, uc_lat, uc_lon):
54 #
55 def _stfun(fun, left, right):
56 if left is None:
57 return right
58 return fun((left, right))
59 #
60 ydis = (2. / 3 / 8 / 10)
61 xdis = (1. / 8 / 10)
62
63 # adjust startpoint to center of mesh
64 # (for special boundary case.)
65 lc_lat = int(lc_lat / ydis) * ydis + ydis / 2.
66 lc_lon = int(lc_lon / xdis) * xdis + xdis / 2.
67
68 ycells = int(np.ceil(uc_lat / ydis) - np.floor(lc_lat / ydis))
69 xcells = int(np.ceil(uc_lon / xdis) - np.floor(lc_lon / xdis))
70 result, lc, uc = None, None, None
71 #
72 for i in range(ycells):
73 tmp = None
74 for j in range(xcells):
75 e, lc_tmp, uc_tmp = load(
76 lc_lat + ydis * i,
77 lc_lon + xdis * j)
78 #
79 if not lc:
80 lc = lc_tmp
81 uc = uc_tmp
82 #
83 tmp = _stfun(np.hstack, tmp, e)
84 #
85 result = _stfun(np.vstack, result, tmp)
86 #
87 return result, lc, uc
88
89
90 def load_range(lc_lat, lc_lon, uc_lat, uc_lon):
91 elevs, lc, uc = _load_range0(lc_lat, lc_lon, uc_lat, uc_lon)
92 X = np.linspace(lc[1], uc[1], elevs.shape[1])
93 Y = np.linspace(lc[0], uc[0], elevs.shape[0])
94 xmin, xmax = np.where(X <= lc_lon)[0][-1], np.where(X >= uc_lon)[0][0]
95 ymin, ymax = np.where(Y <= lc_lat)[0][-1], np.where(Y >= uc_lat)[0][0]
96 return X[xmin:xmax], Y[ymin:ymax], elevs[ymin:ymax, xmin:xmax]
せっかくなので、これまでやってなかった場所で、かつ、上記変換対象になった場所で確認。「38.235130, 140.851829」～「38.336993, 140.954586」。これは、南西に仙台城址付近と若林城址(現宮城刑務所)が両方収まるようにし、北東端が宮城ひとめぼれスタジアムとした範囲：
Click here to download
おけ。
一番右上の三日月上のものは地形ではなく「宮城ひとめぼれスタジアム」。日韓共同開催のサッカーワールドカップで日本がトルコに負けた「あの宮城スタジアム」(「負け試合のスタジアム」)です。ここから少し南西にくだった高いところが「高森山」で、南北朝時代からの城。なんとも長いこと形を留めていたものだと思うが、3.11 でかなり崩れた。身長を超えるような大岩が崩れてきててビビった。その高森山から西へいったところが「鶴が丘団地」で、団地のへりにある山っぽい地形(3つ東西に並んだ山の南の山)が「鶴ヶ城」だったりするんですが、政宗の政敵ともなった政宗の叔父の城で、なんというかちょっと軽視されてる城です 
一応上空から見ると鶴のように見えるから、てことなんだけど、どうかな? (なお、鶴ヶ城の北にある平地は鶴が丘中学校・松森小学校の敷地。)
この高森城やら鶴ヶ城の南に流れる川が七北田川。南西端でうねっているのが無論広瀬川。何段にも重なっている河岸段丘を堪能してみて。真ん中の川は梅田川だけど、仙台市民でも梅田川を意識してる人はあんまりいない気がする。
なんでこんな範囲かというと。ワタシ、高校時代、まさにこの右上らへんから左下付近(仙台城址の下)まで、自転車で毎日通学してたから、なの。低いとこまわり込んだんでしょ、って? 違います。真ん中の丘陵部を超えて通ってたのです。
サンドウィッチマンの伊達の実家のある「八乙女」というあたりから「登って」、台原を越えて「一気に下る」。北西のほうに、斜めに長い池があるよね? この池の南東端の経度から真っ直ぐ七北田川に向かって南下して七北田川にぶつかるところが七北田公園、Perfume のツアーで使われたイズミティ２１はこの間にあります。で、七北田公園の池から東の方に、川の二又の少し上流に線状の高まりがあると思いますが、これが22号線(奥州街道)。これがワタシが通っていた道路。七北田川本流を渡り、さらに支流を渡ったところが八乙女。なお、今の泉中央駅が出来る前は、仙台市営地下鉄は八乙女駅が終点でした。もう一本支流があるけど、これは渡らずに、けれどもこの支流が刻んだ谷を登っていくのが奥州街道。この奥州街道、谷を入ってしばらくは識別出来るね。谷に沿って登っていくと、宅地化の形跡がない自然地形があると思うけど、これが台原森林公園。公園手前あたりで川渡ってるのか。気付いたことは一度もない。自転車でこの谷沿いの街道を走ってる限りは、川を渡ったことに気付かないのと同様に、「あぁ台原森林公園だ」と思うような識別は出来ないんだけれど、とにかくこれを左手に見ながら(見てなくても)登り続ける。ここらへんて、自動車で、特にドライバーじゃなくただ同乗してると「坂だ」と思わないくらいのゆる～い坂でさ、自転車でのツラさ、あんまりわかってもらえないのよねぇ。一度として下らない「無限上昇」が10分以上続くわけですよ。しかも最後の数分だけは自動車でもはっきり「あぁ坂だ!」とわかる急坂。最後の最後に一気にツラい。登りきったらその下りも急で、これはまさに一気に下ってしまう。
「支流が刻んだ谷を登っていく」なんて、あたかも前から知ってたみたいに書いたけど、今こうやって地形を可視化してみてはじめて知ったこと。なるほど、だから最後の登りと上りきってからの下りだけが急だったのか。やっぱり知ってる場所の地形を可視化するのがおもしれーや。
真ん中の一番下にある四角い人工物が若林城址(現宮城刑務所)、その北にある人工物が宮城球場(楽天Koboスタジアム宮城)を含む「宮城野原公園」。このスタジアムの東を走っている線が東北本線の貨物線、ずっと西に行って同じような線が、これが東北本線で、そのちょっと平らな地形にみえてる部分が仙台駅。南西の、範囲内で一番赤領域が広い山が「八木山」で、さま～ずが「モヤモヤさま～ず2」で訪れた八木山ベニーランドのある場所。その谷挟んで北が仙台城址のある青葉山。広瀬川を一番南から北上して最初の緩いカーブのところの左側(右岸)に特徴的な細長い盛り上がりがあるかと思いますが、これが愛宕神社。この愛宕神社の西の似たような特徴的な細長い地形のところにあるのが瑞宝寺だね。政宗墓地のあるところです。実は行ったことがない…。
「宮城野原公園」の北西にあるこんもりした場所は、これは榴岡公園のとこか。こんな地形なのか…。そういやあそこ、坂下るな。思い出したわ。この榴岡公園脇も仙台駅に車で向かうのに良く使う道。45号(仙塩街道)から折れて榴岡公園の通りに出ると、まっすぐ仙台駅東口に着く、のね。
ま、ブラタモリを観た人が楽しめる範囲かなと思う。拡大して、じっくり標高が微妙に高くなってくる箇所を眺めてみるといい。
今気付いたんだけど、仙台城と若林城のちょうど真ん中に愛宕神社、というか一直線に並んでるんだな。狙ったんだろうか? 狙ったならすげーぞ、政宗。愛宕神社なんかは多分仙台市民しか知らないと思いますが…、景色いいとこなんで、仙台を訪れて、時間がある方は行ってみれば? 愛宕神社からの景色は、仙台の印象的な風景写真の結構な割合を占めてたりします。仙台市民はよく「仙台には見るべき観光地はない」と良く言うんですが、んなこたぁないと思うよ。住んでた頃はわからなかったけど。
文章では説明しずれぇなぁ、と思ってたんだけど、良く考えたらアタシ、Acorobat を持ってるんだった。文字入れれるのよね。そうすりゃ良かったと思ったがあとの祭り。ま、持ってる人は自分でそうすればいいさ。(追記： (8.5)でやってみた。)
さて、最後に、converter.py の改造：
1 # -*- coding: utf-8 -*-
2 import re
3 import os
4 try:
5 # for CPython 2.7
6 import cPickle as pickle
7 except ImportError:
8 # Python 3.x or not CPython
9 import pickle
10 from os import path
11 from distutils.dir_util import create_tree
12
13 _RGXES = (
14 (re.compile(r"<gml:lowerCorner>(.*) (.*)</gml:lowerCorner>"), "lc", float),
15 (re.compile(r"<gml:upperCorner>(.*) (.*)</gml:upperCorner>"), "uc", float),
16 (re.compile(r"<gml:low>(.*) (.*)</gml:low>"), "L", int),
17 (re.compile(r"<gml:high>(.*) (.*)</gml:high>"), "H", int),
18 (re.compile(r"<gml:startPoint>(.*) (.*)</gml:startPoint>"), "st", int),
19 )
20
21 def _convert_xml(srcgml):
22 _meta = {}
23 _elevs_raw = []
24
25 with open(srcgml, "r") as fi:
26 # grid infos
27 ri = 0
28 while True:
29 line = fi.readline()
30 m = _RGXES[ri][0].search(line)
31 if m:
32 _meta[_RGXES[ri][1]] = (
33 _RGXES[ri][2](m.group(1)), _RGXES[ri][2](m.group(2)))
34 ri += 1
35 if len(_RGXES) == ri + 1:
36 break
37
38 # find start of grid data
39 while True:
40 line = fi.readline()
41 if "gml:tupleList" in line:
42 break
43
44 # read grid data
45 while True:
46 line = fi.readline()
47 if "gml:tupleList" in line:
48 break
49 ty, d = line.strip().split(",")
50 # we need to treat sea level always NaN rather than zero.
51 if ty.decode('cp932') in (u'データなし', u'海水面'):
52 d = '-9999.'
53 # invalid data are stored with '-9999.0'
54 if float(d) == -9999.0:
55 d = u'NaN'
56 _elevs_raw.append(float(d))
57
58 # read start-point
59 while True:
60 line = fi.readline()
61 m = _RGXES[ri][0].search(line)
62 if m:
63 _meta[_RGXES[ri][1]] = (
64 _RGXES[ri][2](m.group(1)), _RGXES[ri][2](m.group(2)))
65 break
66
67 shape = (_meta["H"][1] + 1, _meta["H"][0] + 1)
68 _meta["st"] = _meta["st"][1] * shape[1] + _meta["st"][0]
69
70 # pickle meta and data
71 rgx = re.compile(r"FG-GML-(.*)-(.*)-(.*)-.*-.*.xml")
72 m = rgx.match(path.basename(srcgml))
73 dest_dir = "__converted/{}/{}".format(m.group(1), m.group(2))
74 dest_filebase = "{}".format(m.group(3))
75 create_tree(dest_dir, [dest_filebase])
76 data = {"startPoint": _meta["st"], "elevations": _elevs_raw}
77 pickle.dump(
78 data,
79 open(path.join(dest_dir, "{}.data".format(dest_filebase)), "wb"),
80 protocol=-1)
81
82
83 #
84 for root, dirs, files in os.walk("__extracted"):
85 for fn in files:
86 srcgml = path.join(root, fn)
87 _convert_xml(srcgml)
88 print(srcgml)
これもせっかくなので、別の場所。茨城は持ってなかったので、今ダウンロードした。筑波山てみる?
Click here to download
おけ。「36.199884, 140.068490」～「36.320861, 140.169969」の、南端が筑波山。
茨城なんてのは、水戸に一度行ったきりの、新幹線で通過するだけの場所、なので、なんにも話せることはないです。いや、「いばらき!」と県民に突っ込まれるくらいの印象しかない、というか。
あー、どうせならワタシが行ったことがある場所がいいな。「小諸城」てみよう：
Click here to download
うーん、谷が深過ぎて、光源調整しないとよくわかんないのよねぇ…。とても印象的な地形(人工的でもあり、自然でもある不思議な)だったので、これで今ひとつ伝わりにくいのが残念。
2015-10-10 追記
最初の投稿時点で2箇所の問題があった。
converted_loader.py の、「adjust startpoint to center of mesh」部分が10日に追加した箇所。この調整がない場合、緯度経度がちょうどメッシュ境界だった場合の振る舞いがおかしくなる。
もう一つが、同じく converted_loader.py の最後の行、np.where部分。最初の投稿時点では等号が入っていなかった。これが入っていないと IndexError が起こりうる。
