(5)の、本題の続き。
(5)で「ひとまず」やってみた隣接セルの結合を整理するのが最初の目標。(5)で作った converted_loader.py の末尾に以下を追加：
1 def _load_range0(lc_lat, lc_lon, uc_lat, uc_lon):
2 #
3 def _stfun(fun, left, right):
4 if left is None:
5 return right
6 return fun((left, right))
7 #
8 ydis = (2. / 3 / 8 / 10)
9 xdis = (1. / 8 / 10)
10
11 # adjust startpoint to center of mesh
12 # (for special boundary case.)
13 lc_lat = int(lc_lat / ydis) * ydis + ydis / 2.
14 lc_lon = int(lc_lon / xdis) * xdis + xdis / 2.
15
16 ycells = int(np.ceil(uc_lat / ydis) - np.floor(lc_lat / ydis))
17 xcells = int(np.ceil(uc_lon / xdis) - np.floor(lc_lon / xdis))
18 result, lc, uc = None, None, None
19 #
20 for i in range(ycells):
21 tmp = None
22 for j in range(xcells):
23 e, lc_tmp, uc_tmp = load(
24 lc_lat + ydis * i,
25 lc_lon + xdis * j)
26 #
27 if not lc:
28 lc = lc_tmp
29 uc = uc_tmp
30 #
31 tmp = _stfun(np.hstack, tmp, e)
32 #
33 result = _stfun(np.vstack, result, tmp)
34 #
35 return result, lc, uc
書き始めるまで結構時間かかるかなぁと思ったけど、書き始めたら10分くらいで書けてしまった。
範囲指定のインターフェイスは悩ましくて、理想形は「中心と距離で」とかかもしれないのだが、緯度経度差の距離換算はここでなくても出来るので、欲しい矩形の西南端・北東端を渡す形にした。_load_range0 という「プライベートちっくな」名前なのは、これが目的の最終形ではないから。このままであれば、(5)と同じことしか出来ない。そりゃそうだ、単に前回のを「すま～とに」書いただけだもん。ともあれこれを使って、(5)の最後のと同じ矩形を：
1 # -*- coding: utf-8 -*-
2 import numpy as np
3 import converted_loader
4
5
6 lat, lon = 35.224167, 139.025373 # 箱根駒ケ岳
7 elevs, lc, uc = converted_loader._load_range0(
8 lat, lon - 1 / 8. / 10.,
9 lat + 2. / 3 / 8 / 10, lon)
10 #elevs, lc, uc = converted_loader._load_range0(
11 # lat, lon - 6 * 1 / 8. / 10.,
12 # lat + 6 * 2. / 3 / 8 / 10, lon)
13
14
15 #
16 import matplotlib
17 import matplotlib.cm as cm
18 import matplotlib.pyplot as plt
19 from matplotlib.ticker import FormatStrFormatter
20
21 fig, ax = plt.subplots()
22 ax.xaxis.set_major_formatter(FormatStrFormatter('%.4f'))
23 ax.yaxis.set_major_formatter(FormatStrFormatter('%.4f'))
24
25 X = np.linspace(lc[1], uc[1], elevs.shape[1])
26 Y = np.linspace(lc[0], uc[0], elevs.shape[0])
27 CS = plt.contour(X, Y, elevs)
28 plt.clabel(CS, inline=1, fontsize=10)
29 plt.show()
同じ矩形から同じ絵は描けたのでそれは貼り付けない。かわりに上の引用でコメントアウトしてある「6 * 」のほう：
よし、っと。ここまでで、「自在なセル結合」が出来るようになった、というわけである。
そうなんだけど、当然「与えられたセルの倍数」でしか扱えないのは困る。ので、ゴールはこれをさらにクリッピングすることである。_load_range0 という名前はそれを反映してたわけだわ。このクリッピングまで行うものをこそ load_range と呼びたいわけです。
てなわけで、新しい「load_range」を含む、新 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
33 _m = pickle.load(open(base + ".meta", "rb"))
34 _raw = pickle.load(open(base + ".data", "rb"))
35
36 elevs = np.zeros(_m["shape"][1] * _m["shape"][0])
37 elevs[_m["st"]:len(_raw) + _m["st"]] = _raw
38 elevs = elevs.reshape((_m["shape"][0], _m["shape"][1]))
39 elevs = np.flipud(elevs)
40 return elevs, lc, uc
41
42
43 def _load_range0(lc_lat, lc_lon, uc_lat, uc_lon):
44 #
45 def _stfun(fun, left, right):
46 if left is None:
47 return right
48 return fun((left, right))
49 #
50 ydis = (2. / 3 / 8 / 10)
51 xdis = (1. / 8 / 10)
52
53 # adjust startpoint to center of mesh
54 # (for special boundary case.)
55 lc_lat = int(lc_lat / ydis) * ydis + ydis / 2.
56 lc_lon = int(lc_lon / xdis) * xdis + xdis / 2.
57
58 ycells = int(np.ceil(uc_lat / ydis) - np.floor(lc_lat / ydis))
59 xcells = int(np.ceil(uc_lon / xdis) - np.floor(lc_lon / xdis))
60 result, lc, uc = None, None, None
61 #
62 for i in range(ycells):
63 tmp = None
64 for j in range(xcells):
65 e, lc_tmp, uc_tmp = load(
66 lc_lat + ydis * i,
67 lc_lon + xdis * j)
68 #
69 if not lc:
70 lc = lc_tmp
71 uc = uc_tmp
72 #
73 tmp = _stfun(np.hstack, tmp, e)
74 #
75 result = _stfun(np.vstack, result, tmp)
76 #
77 return result, lc, uc
78
79
80 def load_range(lc_lat, lc_lon, uc_lat, uc_lon):
81 elevs, lc, uc = _load_range0(lc_lat, lc_lon, uc_lat, uc_lon)
82 X = np.linspace(lc[1], uc[1], elevs.shape[1])
83 Y = np.linspace(lc[0], uc[0], elevs.shape[0])
84 xmin, xmax = np.where(X <= lc_lon)[0][-1], np.where(X >= uc_lon)[0][0]
85 ymin, ymax = np.where(Y <= lc_lat)[0][-1], np.where(Y >= uc_lat)[0][0]
86 return X[xmin:xmax], Y[ymin:ymax], elevs[ymin:ymax, xmin:xmax]
今度の load_range は X, Y を返す。lower corner, upper corner のかわりに。なぜって、elevs クリップの根拠になる X, Y の各々クリップと一緒にやっちまえるからさ。
これを使って箱根駒ケ岳：
1 # -*- coding: utf-8 -*-
2 import numpy as np
3 import converted_loader
4
5
6 lat, lon = 35.224167, 139.025373 # 箱根駒ケ岳
7 # この lat, lon が中心に来るように…
8 lon += 1 / 8. / 10. / 2
9 lat -= 2. / 3 / 8 / 10 / 2
10
11 X, Y, elevs = converted_loader.load_range(
12 lat, lon - 1 / 8. / 10.,
13 lat + 2. / 3 / 8 / 10, lon)
14
15
16 #
17 import matplotlib
18 import matplotlib.cm as cm
19 import matplotlib.pyplot as plt
20 from matplotlib.ticker import FormatStrFormatter
21
22 fig, ax = plt.subplots()
23 ax.xaxis.set_major_formatter(FormatStrFormatter('%.4f'))
24 ax.yaxis.set_major_formatter(FormatStrFormatter('%.4f'))
25
26 CS = plt.contour(X, Y, elevs)
27 plt.clabel(CS, inline=1, fontsize=10)
28 plt.show()
出来た絵はこう：
アスペクト比の制御をしてないんでわかりにくいんだけれど、正方である必要はなく、たとえば 35.650601, 139.249413、35.657680, 139.272501 (八王子城付近) のように横長を指定してこう：
ちなみにこの範囲を 5.5でやってみた LightSource するとカッチョいいぞ：
さてさて。これでかなり自由に欲しい場所が取れるようにはなった。データの取り扱いについては、おおよそ完成間近。残るは、やってもやらなくてもいいかもしれない「データ量削減」。かなりの広域を扱おうとすると、当然メモリがもたないし、遅すぎて処理が終わらないなんてことも起こるが、その場合でも単なる間引きではなく平均を取るとかしたい、よね。けどこれは当面やらない。当座、小さな範囲の地形のほうに興味いってるから。
データそのものの扱い以外で残っているのはまさに「お絵描き」のほう。表示等高線間隔を自在に操りたいわけよね。それによって、立体模型を作る際の精密さ(もしくは製作の難易度)を制御したいわけです。ここいらは純粋に matplotlib の世界。これが次なる課題、ね。
2015-10-10 追記
最初の投稿時点で2箇所の問題があった。
converted_loader.py の、「adjust startpoint to center of mesh」部分が10日に追加した箇所。この調整がない場合、緯度経度がちょうどメッシュ境界だった場合の振る舞いがおかしくなる。
もう一つが、同じく converted_loader.py の最後の行、np.where部分。最初の投稿時点では等号が入っていなかった。これが入っていないと IndexError が起こりうる。
