本題なんだか寄り道なんだかわからない。本題のような気もする。「回転」の話。
河川に沿って、とか、尾根に沿って、という見方をしたいこと、あるでしょう? たとえばこんな領域に興味があるとする：
Click here to download
これは「35.580003, 139.626548」～「35.629057, 139.669465」。けど、「見たい」のは多摩川の北側(都内側)だけで、なおかつ、都内側のもっと北東側を広げてみたいのだ、とする。ワタシの場合は、田園調布・尾山台あたりの変化に富んだ地形が面白いので。つまりは、多摩川が上下に真っ直ぐになるような回転をして、真っ直ぐになった多摩川の右側(左岸)を広げたい、ということ。
緯度経度の tick にこだわったまま回転するのは難しいと思うけれど、imshow は所詮「絵」で、これの回転は簡単…、と言いたいところだけど、numpy だけでやる場合は、結構自力で(回転の公式とかアフィン変換使って)やらないといけない。scipy 使うのが楽。(90度回転は何したって出来るよ。rot90 とか、あるいは transpose (転置 A.T) が目的に適う場合もある。)
1 # -*- coding: utf-8 -*-
2 import numpy as np
3 import converted_loader
4
5 minlat, minlon = 35.580003, 139.626548
6 maxlat, maxlon = 35.629057, 139.669465
7
8 X, Y, elevs = converted_loader.load_range(
9 minlat, minlon,
10 maxlat, maxlon)
11
12 elevs = np.flipud(elevs)
13
14 import scipy.ndimage as ndimage
15 elevs = ndimage.rotate(
16 elevs,
17 -55.,
18 reshape=False,
19 cval=np.nan,
20 prefilter=False)
21 # ...
デフォルトのprefilter=Trueがなんかダメで。reshape は使いやすいほうを使えばいい。reshape=True だと回転で溢れる有効領域全部を収めてくれる。cval は元の領域の外側だった部分の埋め草値。
結果は当然、こういう、期待とは違う絵になる：
Click here to download
とくに多摩川台公園(上の絵の一番下)の右側(実際の方位の北東)方向が欲しいわけだから、領域を広げねば：
1 # -*- coding: utf-8 -*-
2 import numpy as np
3 import converted_loader
4
5 minlat, minlon = 35.580003, 139.626548
6 maxlat, maxlon = 35.629057, 139.669465
7 d1 = (maxlat - minlat)
8 d2 = (maxlon - minlon)
9 d3 = (maxlat - minlat)
10 d4 = (maxlon - minlon)
11 X, Y, elevs = converted_loader.load_range(
12 minlat - d1, minlon - d2,
13 maxlat + d3, maxlon + d4)
14
15 elevs = np.flipud(elevs)
16
17 import scipy.ndimage as ndimage
18 elevs = ndimage.rotate(
19 elevs,
20 -55.,
21 reshape=False,
22 cval=np.nan,
23 prefilter=False)
のでこんな：
Click here to download
もちろんこれでは大き過ぎる。クリップして四角形にしたいのだ。
なんか試行錯誤ってのもダサいけど、こんな：
1 shape = elevs.shape
2 elevs = elevs[
3 int(shape[1] / 2 * 1):int(shape[1] / 5 * 4),
4 int(shape[1] / 5 * 2):int(shape[1] / 3 * 2)]
というわけで：
Click here to download
うむ、満足じゃ。ほかにもね、石垣一夜城とか八王子城なんかさ、尾根筋を上下か左右かに真っ直ぐにしてみたかったりするわけだわ。というか、その真っ直ぐにした尾根を基準にした範囲でみたいわけ。その方が「その尾根に立ったとした場合の景色」が想像つきやすいでしょう?
もうちょっとまともな計算でスマートに出来るとは思うけれど。回転した矩形の内接する矩形を見つける問題、だね。けどまぁ、「この場所を詳しく見たいのだ」という個人的目的でもって「お好きな絶妙な範囲で」やりたいわけだから、多少の試行錯誤は許容出来るであろう。
