Authagraph とは
Authagraphとは世界地図の投影法の一種です。各地域を比較的歪みの少ない形と正しい面積比率で表すことのできる投影法として考案されました。全世界を少ない歪みで描ける投影法としてはDymaxion mapがありますが、Dymaxion mapでは海に切れ目が入ってしまうため世界全体のイメージが掴みにくい感じがします。一方、Authagraphでは大洋が一続きになり、全体が長方形に表示される点が優れています。
投影方法
世界地図は、まず正四面体に投影されて平面に展開されます。その投影方法と計算式は作者である鳴川肇氏の論文「正多面体図法を用いた歪みの少ない長方形世界地図図法の提案」と「オーサグラフ図法の数式化と歪み評価」に記載されています。正四面体の表面は24個の三角形に分割されており、上記論文の数式は各々の三角形の局所座標で記述されています。
描画プログラム
GeoJSONのデータからAuthagraphの世界地図を描くPythonプログラムを作成してみました。地図はSVGフォーマットで出力されます。
プログラムの中で用いているナンバーリングについて下の展開図で説明します。大きい数字で示したのが正四面体の頂点です。例えば0番の頂点の周りの三角形に0から5の番号を振ります。以下同様に24個の三角形に図のように番号を振っています。
使い方：
python authagraph.py foo.geojson bar.svg
ArcGIS HubのGeoJSONデータで検証すると下記のような絵が描けます。
まとめ
Authagraphは素晴らしい投影法なのですが、データの可視化等で使われるケースはほとんどありません。それはQGISやR等のGISソフトウェアでサポートされていないからでしょう。誰かこれを読んだ人が実装してくれないかなあ
authagraph.py
'''
Authagraph world map drawing program
usage: python authagraph.py foo.geojson bar.svg
'''
import sys
import json
import numpy as np
#(degree,minute,second,a) to radian
# positive a:East,North, negative a:West,South
def dms2r(d, m, s, a):
return((d + (m + s/60.)/60.)/180.*np.pi*np.sign(a))
#[longitude, latitude] to 3D coodinate on unit sphere
def ll2c(ll):
x = np.cos(ll[1])*np.cos(ll[0])
y = np.cos(ll[1])*np.sin(ll[0])
z = np.sin(ll[1])
return(np.array([x, y, z]))
#input:x, p0, p1 in 3D coodinate
#output:[longitude, latitude] in the local coodinate
# with p0 on the north pole , p1 on the prime meridian
def locll(x, p0, p1):
p2 = p1 - p0 * np.dot(p0, p1)
p2 = p2 / np.linalg.norm(p2)
x2 = np.cross(p2, x)
c1 = np.dot(p0, x)
c2 = np.dot(p2, x)
s2 = np.dot(p0, x2)
return(np.array([np.arctan2(s2, c2), np.arcsin(c1)]))
#coodinate conversion in 1/24 sphere surface triangle
#input: local coodinate [longitude, latitude]
#output: transformed coodinate, triangle number (0 to 5)
def convtri(a):
lm, ph = a
nn = int((lm + np.pi/3)//(np.pi/3)) % 6
lm = (lm + np.pi/3) % (np.pi*2/3) - np.pi/3
r2 = np.sqrt(2)
r3 = np.sqrt(3)
csph = np.cos(ph)
c = (2+np.cos(lm))*csph/(r2*csph+np.sin(ph))
x = 2/np.pi*c*(lm-np.arcsin(np.sin(lm)/r3))
y = (r2-c)/r3
return(np.array([x, y]), nn)
#global coodinate conversion
#input: global coodinate [longitude, latitude] ,
# authagraph transformation parameters
#output: converted coodinate
def convgl(c, ag):
rr3 = np.sqrt(1/3)
x0 = ll2c(c)
pp, ogp, rtm = ag
i0, c0 = 0, np.dot(x0,pp[0])
for i in range(1,4):
if c0 >= rr3:
break
c1 = np.dot(x0,pp[i])
if c1 > c0:
i0, c0 = i, c1
i1 = (i0 + 1) % 4
a = locll(x0, pp[i0], pp[i1])
x, m = convtri(a)
n = i0 * 6 + m
return(rtm[n] @ x + ogp[n])
def drawpg(pt, svf, ge, ag):
factor, x0, y0, wx, wy, xshift = ge
svf.write(' <polygon stroke="black" fill="white" stroke-width="1"\n')
svf.write(' points="')
for ip in range(len(pt)):
p1 = np.array(pt[ip])/180*np.pi
p2 = convgl(p1, ag)
x = (p2[0] * factor + xshift) % wx + x0
y = (wy - p2[1]*factor) + y0
svf.write(f'{x:g},{y:g} ')
if ip % 8 == 7 :
svf.write('\n ')
svf.write('" />\n')
def main():
vb1, vb2, vb3, vb4 = 0, 0, 1600, 750 # View Box
rw, rh = vb3 / 10, vb4 / 10 #real size
wx = 1500
wy = wx * np.sqrt(3)/4
factor = wx / (np.sqrt(2/3) * 4)
x0, y0 = 50, 50
xshift = -130
#parameters: scale factor, origin, size
ge = [factor, x0, y0, wx, wy, xshift]
#coodinates of the tetrahedron vertices
pp = [ll2c([dms2r(149.,27., 3.56868, 1),dms2r(76.,52.,51.82608, 1)]),
ll2c([dms2r( 97.,21.,25.2126 , 1),dms2r(27.,57., 9.99792,-1)]),
ll2c([dms2r(133.,16.,57.93168,-1),dms2r(22.,55.,41.65104,-1)]),
ll2c([dms2r( 18.,51., 8.037 ,-1),dms2r( 6.,38.,13.37028,-1)])]
#triangle mapping
l0 = np.sqrt(2/3)
v0 = np.array([1.,0.])
v1 = np.array([np.cos(np.pi/3),np.sin(np.pi/3)])
ogp = l0*np.array([v0+v1, v0+v1, 2*v0+v1, 2*v0+v1, 2*v0+2*v1, 2*v1,
0*v0, 2*v0, v0+v1, v0+v1, v1, v1,
3*v0+v1, 3*v0+v1, 2*v0+v1, 2*v0+v1, 2*v0, 4*v0,
2*v1, 2*v0+2*v1, 3*v0+v1, 3*v0+v1, v1, v1])
#rotation angle of triangles
rta = np.pi*np.array([-1/6,-1/6,1/6,1/6,1/2,-1/2,
-1/2,1/2,5/6,5/6,-5/6,-5/6,
5/6,5/6,-5/6,-5/6,-1/2,1/2,
1/2,-1/2,-1/6,-1/6,1/6,1/6])
#triangle rotation matrices
rtl=[]
for i in range(24):
th=rta[i]
rtl.append(np.array([[np.cos(th),-np.sin(th)],
[np.sin(th),np.cos(th)]]))
rtm=np.array(rtl)
#authagraph transformation parameters
ag = [pp, ogp, rtm]
#GeoJSON file
inf = open(sys.argv[1], 'r')
mapjson = json.load(inf)
inf.close()
#SVG header, frame
svf = open(sys.argv[2], 'w')
svf.write('<?xml version="1.0" encoding="UTF-8"?>\n')
svf.write('<svg\n')
svf.write('xmlns="http://www.w3.org/2000/svg"\n')
svf.write('version="1.1"\n')
svf.write(f'viewBox="{vb1:g} {vb2:g} {vb3:g} {vb4:g}"\n')
svf.write(f'width="{rw:g}mm" height="{rh:g}mm">\n')
svf.write(' <rect stroke="black" fill="white" stroke-width="2"')
svf.write(f' x="{x0:g}" y="{y0:g}" width="{wx:g}" height="{wy:g}" />\n')
#draw polygons
nc = len(mapjson['features'])
for ic in range(nc):
tp = mapjson['features'][ic]['geometry']['type']
cd = mapjson['features'][ic]['geometry']['coordinates']
if tp == 'Polygon' :
drawpg(cd[0], svf, ge, ag)
elif tp == 'MultiPolygon' :
na = len(cd)
for ia in range(na):
drawpg(cd[ia][0], svf, ge, ag)
svf.write('</svg>\n')
svf.close()
if __name__ == '__main__':
main()