はじめに
GeoJSONデータからAuthagraph世界地図を描くという記事を書きました。そこで少しDymaxion Map について触れましたが、ほとんど同じプログラムでこれも描く事ができるので、ついでに書いてみました。
Dymaxion Mapとは
Richard Buckminster Fullerが考案した地図の投影法です。世界地図を心射図法で正二十面体に投影し、それを平面に展開したものです。シンプルな発想で、そこそこ少ない歪みで全世界を平面に描くことができます。
正二十面体はどこを切って展開しても構わないのですが、大陸に切れ目が入らないような標準的な切り方が提案されています1。この場合、海には切れ目が入ってばらばらになってしまうので、その広さや形は分かりにくくなっています。また、北極点を中心に左右に展開されるので、方向感とか東西方向の距離感が混乱します。そういった欠点もあるものの、陸地の形をかなり正確に再現できる優れた投影法です。
描画プログラム
GeoJSONのデータからDymaxion Mapを描くPythonプログラムを作成してみました。地図はSVGフォーマットで出力されます。
プログラムの中で用いているナンバーリングについて下の展開図で説明します。
丸囲みの数字が頂点の番号で、三角形の中に記載したものが面の番号です。同じ番号が振られた頂点は立体になった際に重なる点です。12面と16面は2つずつの領域に分割されます。それぞれの面は下記のように6個ずつの小さい三角形に分割されます。
使い方：
python dymaxion.py foo.geojson bar.svg
ArcGIS HubのGeoJSONデータで検証すると下記のような絵が描けます。
まとめ
Dymaxion Mapは輪郭の形の主張が強すぎるので、統計データの可視化などの用途には向かないような気がします。ただ、これほど公平な視点から世界を眺められる地図は他になくて、国連の旗は正距方位図法で描かれていますがDymaxion Mapの方がいいんじゃないかなどと思ってしまいます。
dymaxion.py
'''
Dymaxion world map drawing program
usage: python dymaxion.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/120 sphere surface triangle
#input: local coodinate [longitude, latitude]
#output: transformed coodinate, triangle number (0 to 5)
def convtri(a):
lm, ph = a
n = int(lm//(np.pi/3)) % 6
c = (np.sqrt(5)+3)/(2*np.sqrt(3))/np.sin(ph)*np.cos(ph)
x = np.cos(lm)*c
y = np.sin(lm)*c
return(np.array([x, y]), n)
#global coodinate conversion
#input: global coodinate [longitude, latitude] ,
# dymaxion map transformation parameters
#output: converted coodinate
def convgl(c, ag):
cc = np.sqrt(np.sqrt(5)/6 + 1/2)
x0 = ll2c(c)
pp, fc, cf, ogp, rtm = ag
i0, c0 = 0, np.dot(x0,cf[0])
for i in range(1,20):
if c0 >= cc :
break
c1 = np.dot(x0,cf[i])
if c1 > c0:
i0, c0 = i, c1
i1 = fc[i0, 0]
a = locll(x0, cf[i0], pp[i1])
x, j = convtri(a)
n = i0 * 6 + j
return(rtm[n] @ x + ogp[n])
def drawpg(pt, svf, ge, ag):
factor, x0, y0, wx, wy = 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 + 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, 800 # View Box
rw, rh = vb3 / 10, vb4 / 10 #real size
wx = 1500
wy = wx * 3*np.sqrt(3)/11
factor = wx / 11
x0, y0 = 50, 50
#parameters: scale factor, origin, size
ge = [factor, x0, y0, wx, wy]
#coodinates of the icosahedron vertices
pp = np.array([[0.420152, 0.078145, 0.904083],
[0.995005, -0.091348, 0.040147],
[0.518837, 0.83542, 0.181332],
[-0.414682, 0.655962, 0.630676],
[-0.515456, -0.381717, 0.767201],
[0.355781, -0.84358, 0.402234],
[0.414682, -0.655962, -0.630676],
[0.515456, 0.381717, -0.767201],
[-0.355781, 0.84358, -0.402234],
[-0.995009, 0.091348, -0.040147],
[-0.518837, -0.83542, -0.181332],
[-0.420152, -0.078145, -0.904083]])
fc = np.array([[0,1,2],[0,2,3],[0,3,4],[0,4,5],[0,5,1],
[2,1,7],[3,2,8],[4,3,9],[5,4,10],[1,5,6],
[1,6,7],[2,7,8],[3,8,9],[4,9,10],[5,10,6],
[7,6,11],[8,7,11],[9,8,11],[10,9,11],[6,10,11]])
cf = np.empty([20,3])
for i in range(20):
cf[i,:] = (pp[fc[i,0],:]+pp[fc[i,1],:]+pp[fc[i,2],:])/3
cf = cf * np.sqrt(3*np.sqrt(5)+15)/np.sqrt(3*np.sqrt(5)+7)
#triangle mapping
v0 = np.array([1.,0.])
v1 = np.array([np.cos(np.pi/3),np.sin(np.pi/3)])
vc0 = (v0 + v1) *2/3
vc1 = (2*v1 - v0) *2/3
ogp = np.array([v0+4*v1+vc0]*6 + [3*v0+2*v1+vc1]*6 + [3*v0+2*v1+vc0]*6 +
[5*v0+2*v1+vc1]*6 + [3*v0+4*v1+vc1]*6 +
[v0+4*v1+vc1]*6 + [v0+2*v1+vc0]*6 + [5*v0+vc1]*6 +
[5*v0+2*v1+vc0]*6 + [5*v0+4*v1+vc0]*6 +
[7*v0+4*v1+vc0]*6 + [v0+2*v1+vc1]*6 +
[3*v0+vc1]*4+[3*v0+vc0]*2 +
[7*v0+vc1]*6 + [7*v0+2*v1+vc1]*6 +
[9*v0+2*v1+vc1]*6 + [v0+vc1]*2+[9*v0+2*v1+vc0]*3+[v0+vc1] +
[v0+vc0]*6 + [7*v0+vc0]*6 + [7*v0+2*v1+vc0]*6)
#rotation angle of triangles
rta = np.pi*np.array([-1/6]*6 + [1/6]*6 + [1/2]*6 + [5/6]*6 + [-1/2]*6 +
[-1/2]*6 + [-1/6]*6 + [1/6]*6 + [1/2]*6 + [1/2]*6 +
[1/2]*6 + [1/6]*6 + [1/6]*4+[1/2]*2 + [5/6]*6 + [5/6]*6 +
[1/6]*6 + [1/6]*2+[-1/6]*3+[1/6] +
[-1/6]*6 + [1/2]*6 + [1/2]*6)
#triangle rotation matrices
rtm=np.array([np.empty([2,2])]*120)
for i in range(120):
th=rta[i]
rtm[i]=np.array([[np.cos(th),-np.sin(th)],
[np.sin(th), np.cos(th)]])
#dymaxion map transformation parameters
ag = [pp, fc, cf, ogp, rtm]
outline = np.array([-v0+2*v1] + [v0+v1] + [v0] + [3*v0] + [3*v0+vc1] +
[3*v0+2*v1] + [3*v0+vc0] + [5*v0] + [5*v0+2*v1] +
[7*v0] + [9*v0] + [7*v0+2*v1] + [10*v0+2*v1] +
[9*v0+4*v1] + [7*v0+6*v1] + [7*v0+4*v1] + [5*v0+6*v1]+
[5*v0+4*v1] + [3*v0+4*v1] + [3*v0+6*v1] + [-v0+6*v1] +
[v0+4*v1] + [-v0+4*v1] + [v0+2*v1] + [-v0+2*v1])
#GeoJSON file
inf = open(sys.argv[1], 'r')
mapjson = json.load(inf)
inf.close()
#SVG header
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')
#draw outline
svf.write(' <polygon stroke="black" fill="none" stroke-width="1"\n')
svf.write(' points="')
for ip in range(25):
x = outline[ip,0] * factor + x0
y = (wy - outline[ip,1]*factor) + y0
svf.write(f'{x:g},{y:g} ')
svf.write('" />\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()