(6)からの続きということにはなるけれど、今回のは本題と寄り道両方。寄り道率のほうが高いかも。
この手のってさ、ある一定レベルまで来たら、多少イジめたほうがいいわけね。ちょっと無理して広域、とか。そのお題として、「大田区民にはかなり不満な」国土地理院の「1:25,000デジタル標高地形図「東京都区部」」：
大田区民は肝心なところが切れてて使えないわけね。なわけで、これの、「東京都区南東部」をやってみようかと。むろん地名とかそんなのはムリだけど、「地形」だけは今の手持ちで出来るから。
やってみてすぐに問題発覚。というかわかってたけど。セル(xml一個)が欠落してると、当然なにも処置してないので例外で処理を終えてしまう。千葉も埼玉も持ってなかったし、どうも東京もダウンロードに失敗してたのもあったみたいで。
ゆえ、converted_loader.py をこうした：
1 # -*- coding: utf-8 -*-
2 import pickle
3 import numpy as np
4
5
6 def _calc_pos0(lon_or_lat, lseg):
7 a = (lon_or_lat % 100) / lseg
8 b = (a - int(a)) / (1 / 8.)
9 c = (b / (1 / 10.)) % 10
10 return (
11 int(a),
12 int(b),
13 int(c),
14 round((int(a) + int(b) / 8. + int(c) / 8. / 10.) * lseg, 9),
15 round((int(a) + int(b) / 8. + (int(c) + 1) / 8. / 10.) * lseg, 9))
16
17
18 def _calc_from_latlon(lat, lon):
19 _fromlat = _calc_pos0(lat, 2/3.)
20 _fromlon = _calc_pos0(lon, 1.)
21 return (
22 "__converted/{}{}-{}{}/{}{}".format(
23 _fromlat[0], _fromlon[0],
24 _fromlat[1], _fromlon[1],
25 _fromlat[2], _fromlon[2]),
26 (_fromlat[3], _fromlon[3] + 100),
27 (_fromlat[4], _fromlon[4] + 100))
28
29
30 def load(lat, lon):
31 base, lc, uc = _calc_from_latlon(lat, lon)
32 shape = (150, 225)
33 elevs = np.full((shape[1] * shape[0], ), np.nan)
34
35 try:
36 _m = pickle.load(open(base + ".meta", "rb"))
37 _raw = pickle.load(open(base + ".data", "rb"))
38 elevs[_m["st"]:len(_raw) + _m["st"]] = _raw
39 except IOError as e:
40 print(str(e))
41
42 elevs = elevs.reshape((shape[0], shape[1]))
43 elevs = np.flipud(elevs)
44 return elevs, lc, uc
45
46
47 def _load_range0(lc_lat, lc_lon, uc_lat, uc_lon):
48 #
49 def _stfun(fun, left, right):
50 if left is None:
51 return right
52 return fun((left, right))
53 #
54 ydis = (2. / 3 / 8 / 10)
55 xdis = (1. / 8 / 10)
56
57 # adjust startpoint to center of mesh
58 # (for special boundary case.)
59 lc_lat = int(lc_lat / ydis) * ydis + ydis / 2.
60 lc_lon = int(lc_lon / xdis) * xdis + xdis / 2.
61
62 ycells = int(np.ceil(uc_lat / ydis) - np.floor(lc_lat / ydis))
63 xcells = int(np.ceil(uc_lon / xdis) - np.floor(lc_lon / xdis))
64 result, lc, uc = None, None, None
65 #
66 for i in range(ycells):
67 tmp = None
68 for j in range(xcells):
69 e, lc_tmp, uc_tmp = load(
70 lc_lat + ydis * i,
71 lc_lon + xdis * j)
72 #
73 if not lc:
74 lc = lc_tmp
75 uc = uc_tmp
76 #
77 tmp = _stfun(np.hstack, tmp, e)
78 #
79 result = _stfun(np.vstack, result, tmp)
80 #
81 return result, lc, uc
82
83
84 def load_range(lc_lat, lc_lon, uc_lat, uc_lon):
85 elevs, lc, uc = _load_range0(lc_lat, lc_lon, uc_lat, uc_lon)
86 X = np.linspace(lc[1], uc[1], elevs.shape[1])
87 Y = np.linspace(lc[0], uc[0], elevs.shape[0])
88 xmin, xmax = np.where(X <= lc_lon)[0][-1], np.where(X >= uc_lon)[0][0]
89 ymin, ymax = np.where(Y <= lc_lat)[0][-1], np.where(Y >= uc_lat)[0][0]
90 return X[xmin:xmax], Y[ymin:ymax], elevs[ymin:ymax, xmin:xmax]
DEM5、つまり 5m メッシュであることがわかっていれば、shape は固定なので、メタデータから読むのをやめた。DEM10 を扱いたい場合の話はあとでいい。DEM5 前提の固定値をほかにもバラまいちゃってるしね。np.full 部分はこの投稿参照。
「問題発覚」のもう一つ、領域が大きくなるとなぜかカラーマップがおかしくなる(真っ赤になる、のような)ことが起こっていた。なんだろうと思っていたら、どうやら「無効値」の意味の -9999.0 が、「-9999.」だったりそうでなかったり、てことだった。なんだよ…。(1)、(2)にはこの件追記しといた。で、仕方ないので、持ってるものを変換し直し。あーぁ…。
で、問題解決したところで。西は東急多摩川駅付近、北は品川駅、南は六郷土手付近、東は羽田空港付け根(海老取川)という範囲で PDF にする：
1 # -*- coding: utf-8 -*-
2 import numpy as np
3 import converted_loader
4
5
6 X, Y, elevs = converted_loader.load_range(
7 35.536068, 139.664211,
8 35.629256, 139.755192)
9
10 elevs[np.isnan(elevs)] = np.nanmin(elevs)
11
12
13 #
14 import matplotlib
15 import matplotlib.cm as cm
16 import matplotlib.pyplot as plt
17 from mpl_toolkits.axes_grid1 import make_axes_locatable
18
19 elevs = np.flipud(elevs)
20
21 fig, ax = plt.subplots()
22 fig.set_size_inches(16.53 * 2, 11.69 * 2)
23
24 from matplotlib.colors import LightSource
25 ls = LightSource(azdeg=180, altdeg=15)
26 rgb = ls.shade(elevs, cm.rainbow)
27 cs = ax.imshow(elevs)
28 ax.imshow(rgb)
29
30 # create an axes on the right side of ax. The width of cax will be 2%
31 # of ax and the padding between cax and ax will be fixed at 0.05 inch.
32 divider = make_axes_locatable(ax)
33 cax = divider.append_axes("right", size="2%", pad=0.05)
34 fig.colorbar(cs, cax=cax)
35
36 ax.set_xticks([])
37 ax.set_yticks([])
38
39 plt.savefig("tokyotokubu-se.pdf", bbox_inches="tight")
40 #
ハイライトした10行目には注目してみて。これにもちょっと関係するんだけど、「nanを除いた最小」を求める方法、です、これは。
make_axes_locatable はカラーバーのサイズ調整関係なので、本題とは関係ないけど、情報としては役に立つことでしょう、アーメン。
set_size_inches はイジめポイント。これ、解像度として「A3 横」の2倍(A2かな?)、としてる。最後に bbox_inches=”tight” するので用紙サイズが「A3 横」の2倍になるわけではないけれど、解像度としてはこのサイズに合ったものになる。つまり「せっかくの 5m メッシュ」を出来る限り精細に再現したい、というわけです。
結果はこれ：
Click here to download
Pop-out すると Google Docs で開くので、それで楽しんでみて。多摩川駅付近で切っているので、多摩川台公園もみえてる。(等々力渓谷は入ってない。)これよりもさらに高解像度のが欲しければどうぞ。
じ～っくり観察すると面白いんだよねぇ。地名とかが入ってないからわかりにくいとは思うけど、多摩川のかつての氾濫原の痕跡らしき地形が、おそらく下丸子あたりに見られる。あと「田園調布南」てのもおそらくかつての氾濫原で、清くて正しい田園調布ではない、てことがよくわかる 
ほかにも、丘陵の縁の堆積地形も面白いね。目黒川の河口近くなんかは、三角州状の地形になっているね。
この三角州状地形と同じ標高を辿った線が、およそかつての海岸線なはず。多摩川の北にある細い川は呑川。呑川が90度に屈曲しているところで細い線が交差してるでしょ? これは JR。その屈曲との交差からやや南にくだったところが JR 蒲田駅。で、この JR の線を北に向かって丘陵にぶつかるところが JR 大森駅なのね。その東に平和島競艇場があるのがわかるでしょ。つまり JR は、品川駅から大森駅くらいまでは、大昔の海岸線に沿って走ってるのね。
なかなか丘陵とかの名前を調べるのに苦労するんだけれど、「丘陵の縁」といっている部分、どうやらこれ、武蔵野扇状地の先端らしい。そうだったのか…。
呑川がその武蔵野扇状地の先端に辿り着くところが、(6.5)で扱った池上本門寺だね。この北西にある池らしきものは池です。洗足池ね。なるほど、池が出来る地形だというのが良くわかるわね。
旧東海道はむしろ京急(というか国道15号(京浜第一))に沿ってるのだけど、これが少なくとも大森くらいまでは江戸時代の海岸線にほぼ沿ってたはずなのだけど、地形にはあんまり現れてないんだねぇ。
ほかにも、かつての川の跡と思われる地形がいくつもありますね。目黒川に合流する2本の旧河川地形が面白い。上の細くて深いほうは大崎駅のちょい南くらいか。確かにあの辺り歩いたことあるけど、ちょっと谷地形っぽかったかも、そういえば。うーん、色々歩いて確認してみたくなるぞ、これは。
なお、今ホットな「福山雅治」(桜坂)も範囲に含まれてるんだけど、わかるかなぁ? 東海道新幹線が特定出来て、それが多摩川を渡る場所がわかれば…。
ついでなので、もう一つお遊び。main.py をちょっと書き換える：
1 # ... (snip) ...
2 cut_min = np.nanmin(elevs) + 15.
3 elevs[elevs < cut_min] = np.nan
4 elevs[np.isnan(elevs)] = np.nanmin(elevs)
5 # ... (snip) ...
Click here to download
こういう東京を見たいか見たくないかは、貴方次第です。「現在の地形のまま海水面が上昇したとしたら」というシミュレーション。(なので、これが過去のたとえば縄文海進を模擬するとかにはならない。)
さて、最後に、「もっとイジめる」。つまり、「1:25,000デジタル標高地形図「東京都区部」」がカバーしている領域を全部網羅しつつ、「大田区民が喜ぶ」ものにする。範囲は「35.536068, 139.572353」～「35.771865, 139.890633」。shape は (4244, 5729) という凄まじいサイズとなり、処理時間は 6m27.543s。こうして出来上がったのが以下：
Click here to download
美しい…。地名等々の情報がない以外は「1:25,000デジタル標高地形図「東京都区部」」に負けてませんな。あ、これはかなりデカくて、Google Docs だとズームレベルの自由度が足りないので、ダウンロードして Acrobat Reader で見たほうがいいです。これについても、さらに高解像度のものあります。ただし凄まじいファイルサイズなので注意(36MB)。(大きさの割にはあんまり印象変わらないと思うよ。)
見てわかる通り、この範囲の最高地点は 80m くらいなのであって、これってのは「日本地図」という広域地図では全部が緑の「真っ平」に思ってしまうような高低差なわけね。多くの港町も、訪れてみると全然平らじゃない場所なんかいくらでもありますね。大きな平野だと、ほかにも帯広周辺(十勝平野)なんかも、行ってみると全然「見渡す限りの平原」ないのでビックリするわけな。
関東平野ってのも、見る場所から見れば確かに「真っ平」ではあって、東北新幹線に乗って関東平野に入ると、随分遠くまで見渡せて、そして山が遠いことには確かに驚きます。仙台平野なんてのは、平らな部分は「本当に広くて真っ平、真性の真っ平」なんだけれどもね、そこってのは水田地帯なわけだし、都市部はなにせ内陸奥深くにあって、人が多く居住する場所からは「山が近い」し、広瀬川の大規模な河岸段丘が日常、だったりして、案外「広く感じない」。そう考えれば関東平野、確かに「だだっ広くて真っ平」、なんだけれども。
「東京がそう平らでもない」ことに気付くのって、やっぱ住んでからかなぁと思う。人間スケールてのは、1m でも 2m でも「高低差」として強く感じるわけであって、広域図には現れない「高低差」は、地上を歩く人間目線にはかなりのものなわけだね。というわけで「見ておわかりの通り」、東京には坂道が多い、のです。タモリ
が好む坂道が。
国道15号(平坦)と国道1号(アップダウンだらけ)を使い分けてみると面白いのよね。前者は海沿いの平らなとこを行き、後者は武蔵野扇状地の先端のあたりを通る。その際何度も河川が彫った谷を上り下りするので、何度もアップダウンする。自動車でもわかるくらいの坂なので、自転車は大変さ。
環八、環七も、西側は結構登ります、下ります。等々力渓谷に向かう際に環八に出てしまったことを書いたけど、本当に「違和感」を感じたのは「環八の坂」でした、ワタシには。自動車生活してないんでね、狭い範囲の環八しか知らんかったの。
2015-10-10 追記
最初の投稿時点で2箇所の問題があった。
converted_loader.py の、「adjust startpoint to center of mesh」部分が10日に追加した箇所。この調整がない場合、緯度経度がちょうどメッシュ境界だった場合の振る舞いがおかしくなる。
もう一つが、同じく converted_loader.py の最後の行、np.where部分。最初の投稿時点では等号が入っていなかった。これが入っていないと IndexError が起こりうる。
