まずは、元となる世界地図の画像を読み込んで表示してみよう。
%matplotlib inline import matplotlib.pyplot as plt import numpy as np from PIL import Image img = np.array(Image.open("world_map.jpg")) plt.imshow(img)
PILは画像を開くため、matplotlibは画像を表示するために使っている。ほかの行列の計算などのほとんどの処理はnumpyで行っている。
これからは、画像をnumpyの行列として扱っていく。
この画像（行列）には、どういったデータが含まれているのか見てみよう。
> img.shape (1500, 3000, 3)
これは、画像が1500x3000のサイズであることを表している。3つ目の次元は、そのピクセルのRGB値を格納している。
対象となる画像のデータのフォーマットが分かったところで、さっそくプロジェクトを始めよう。
まずは、この世界地図の画像の各ピクセルの、地球上の緯度と経度を求める。この世界地図は、緯度と経度が等間隔で並ぶ正距円筒図法を用いているため、非常に簡単な式で画像上の座標から緯度と経度が求められる。
緯度と経度をラジアンで求める。画像の幅と高さをそれぞれ\(w\)と\(h\)とすると、地図上の点\( (x, y) \)の経度\(\lambda\)、緯度\(\phi\)は、
$$
\lambda = \frac{2\pi}{w} \left( x - \frac{w}{2} \right) \\
\phi = -\frac{\pi}{h} \left( y - \frac{h}{2} \right)
$$
少し複雑に見えるかもしれないが、やっていることは、xとyを、経度は\(-\pi\)から\(\pi\)の、緯度は\(-\frac{\pi}{2}\)から\(\frac{\pi}{2}\)の範囲に縮めているだけである。\(\phi\)にマイナスがついているのは、画像の原点は左上にあり、y座標は下に行くにつれて大きくなるからである。北極の緯度のラジアンを\(\frac{\pi}{2}\)とするために、マイナスをつけるのである。
緯度と経度が分かったので、これからそれを2次元平面に投影する。
3次元座標空間に、原点を中心とする単位球があるとする。北極の座標は(0, 0, 1)、南極の座標は(0, 0, -1)である。本初子午線（経度0°）はx軸の正の方向を向いている。この単位球が地球であるとして、それを北極からxy平面に投影しよう。
\((\phi, \lambda)\)を平面に投影したときの点の極座標\((r, \theta)\)は、
$$
(r, \theta) = \left( \tan \left( \frac{\pi}{4}+\frac{\phi}{2} \right), \lambda \right)
$$
説明は省略するが、図を書いてくれればわかると思う。
直交座標\((x, y)\)はもちろん、
$$
(x, y) = (r \cos \theta, r \sin \theta)
$$
最後に、単位球のサイズに縮められた座標を元のスケールに戻すために、xに\(\frac{w}{2\pi}\)を、yに\(-\frac{h}{\pi}\)をかけてやる。
しかし、投影された画像は元の画像の鏡像（地球を内側から投影しているわけなので、左右が逆になる）となるため、yに-1をかけてやれば正しい方向になるから、結局yには\(\frac{h}{\pi}\)をかけるだけでよい。
今まで言ったことをプログラムしてみよう。
まずは、必要となる変数を定義する。ここでは、numpy.meshgrid関数を使って、すべての座標の組み合わせをxとyに格納している。
h = img.shape[0] w = img.shape[1] x, y = np.array(np.meshgrid(np.arange(0, w, dtype=int), np.arange(0, h, dtype=int), indexing='xy'))
次に、上で導いた公式に従って、緯度（lat）と経度（long）、そして投影された新たな座標を計算する。
long = 2*np.pi/w * (x - w/2) lat = -np.pi/h * (y - h/2) r = np.tan(np.pi/4 + lat/2) x_new = r * np.cos(long) y_new = r * np.sin(long) x_new = np.round(x_new * (w/(2*np.pi))) y_new = np.round(y_new * (h/np.pi))
次に、新たな画像の大きさ（d）を定義し、原点が画像の中心に来るように平面全体をシフトする。そして、画像からはみ出た部分を切り取る。ステレオ投影では、北半球の高緯度のところは原点から大きく離れるので、今回はd=2000として、南極の周り±1000ピクセルを表示する。
d = 2000 x_new = np.int_(x_new + d / 2) y_new = np.int_(y_new + d / 2) x_new = np.minimum(x_new, d - 1) y_new = np.minimum(y_new, d - 1) x_new = np.maximum(x_new, 0) y_new = np.maximum(y_new, 0)
最後に、投影された座標に、元の座標のRGB値を代入し、得られた画像を表示する。
img_new = np.full((d, d, 3), [255, 255, 255]) img_new[y_new, x_new] = img[y, x] plt.imshow(img_new)
見ての通り、南極の周辺は大陸の形がはっきりと美しく見えるが、中心から離れるにつれてピクセルが散らばってゆく。そもそも、有限の表面積を持つ地球上の点を、無限に広がる平面に投影しているのだから、散らばっていくのは当然だ。しかし、正しい形の南極大陸が正しい向きに現れたので、上に示した数学が正しいということが証明された。
他にも、緯度と経度の代わりに球面座標を使って計算する方法も考えた。緯度と経度を使う方法とほとんど変わらないが、少々数式がきれいになるので、興味がある人は挑戦してみるといい。この記事では紹介しないが、GitHub上では公開してある。
次回は、点が散乱しないように投影するにはどうしたらいいか考える。
