Reproducing the AuthaGraph World Map

The AuthaGraph world map by Hajime Narukawa took the media by storm when it was honoured with the Japan Good Design Grand Award in October 2016. After several rounds of copying the news, the journalists dubbed it “the most accurate world map ever”, comparing it with the usual straw man, Mercator projection. Hype aside, the facts are: the map almost avoids the chopping of land, almost has the equal-area property, and does a decent job with land shapes. Its large distortions near the vertices of the basal tetrahedron and drastic discontinuities along the boundary are masked by the ocean.

Since the map projection is not public, I set out to develop its lookalike with PROJ.4. My way is simpler than the AuthaGraph method at the cost of losing the proportions of areas that AuthaGraph claims to preserve. You can compare the results below. Among the most visible changes, my Iceland lies on the European side of the map, Australia properly looks smaller than China, and Brazil’s shape is distorted differently. Given the overall distortions of AuthaGraph, treating Earth as an ellipsoid would not buy us much, so my method assumes it is a sphere. The code is available on Bitbucket.

Original AuthaGraph

pseudo-AuthaGraph

The top border of the original AuthaGraph map corresponds to a great-circle arc that separates Africa from Iceland. If the border passed through the North Pole, it would have to follow a meridian. Any meridian near Africa with westernmost longitude 17° 33′ W and Iceland with easternmost longitude 13° 16′ W must pass through either land. As the border avoids land, it must not pass through the North Pole, which is projected inside the map. Therefore, in spite of the claims about tiling the plane with the map, juxtaposing two maps yields a map with two North Poles close to each other. Incidentally, even the Bering Strait is blocked for a land-avoiding meridian by St. Lawrence Island and Unalaska.

The map of the Arctic below shows the possible locations of the northern vertex with a land-avoiding top border. My northern vertex lies around 73° N 121° E, near the delta of Lena.

arctic

Another constraint is to keep Australia and New Zealand in the middle triangle of the map. The original map does not fulfill it entirely: the Stewart Island lies in its Antarctic triangle. Therefore, I based the projection on a disphenoid, an irregular tetrahedron that can be rotated, scaled, and translated so that its vertices have Cartesian coordinates (+a, +b, +c), (+a, -b, -c), (-a, +b, -c), and (-a, -b, +c). A disphenoid can still be developed to a parallelogram so we can cut a right triangle or trapezoid near the acute-angle vertex of its net and attach it to the opposite side, making the map rectangular.

My approach to projecting a point P from a sphere onto a plane begins with finding the spherical triangle ABC, ABD, ACD, or BCD that contains P.

Given a tetrahedron KLMN inscribed into the unit sphere, the point P lies inside the spherical triangle KLM if \left(P \cdot (K \times L)\right)\left(N \cdot (K \times L)\right) \le 0, \left(P \cdot (K \times M)\right)\left(N \cdot (K \times M)\right) \le 0, and \left(P \cdot (L \times M)\right)\left(N \cdot (L \times M)\right) \le 0.

I considered several ways to map spherical triangles onto planar triangles, keeping in mind that the projection should be continuous on the edges shared by two triangles and smooth in their interiors. Here are my attempts:

  • Barycentric coordinates.

    The barycentric coordinates of point P with respect to a spherical or planar triangle KLM are the ratios of areas of spherical or planar triangles: (\Delta PLM / \Delta KLM, \Delta KPM / \Delta KLM, \Delta KLP / \Delta KLM).

    I was finding a point P' with planar barycentric coordinates corresponding to the spherical barycentric coordinates of P. Unfortunately, this elegant mapping is continuous on the shared edge of two triangles only if they are symmetric with respect to this edge. In an irregular tetrahedron, this is not the case.

  • Gnomonic projection, that is casting from the centre of the sphere onto the plane KLM. Although this projection is continuous on the edges of triangles, I find its distortion unacceptable.
  • Find points K', L', and M' on the sphere, being the intersections of great circles LM with KP, KM with LP, and KL with MP. Express the location of point P as a triple \left(\lvert KM'\rvert / \lvert KL \rvert, \lvert LK'\rvert / \lvert LM \rvert, \lvert ML'\rvert / \lvert MK \rvert\right). On a plane, the three lines that correspond to this triple in general do not cross in one point. Treating the centre of the intersection triangle as P' leads to discontinuities along the edges. I took the incentre of the intersection triangle instead.

The above method produces the map shown below, with sharp turns on the edges of the basal tetrahedron. I smoothed the joints by interpolation.

authagraph1

The cheetah spots on the map below are Tissot indicatrices for the projection, magnified images of infinitesimal circles on Earth’s surface.

The Tissot indicatrix is an ellipse with scale factor along the meridian h = \sqrt{(\partial x/\partial\phi)^2 + (\partial y/\partial\phi)^2}/R, scale factor along the parallel k = \sqrt{(\partial x/\partial\lambda)^2 + (\partial y/\partial\lambda)^2}/(R\cos\phi), and the angle \theta' at which the meridian and parallel intersect given by \sin\theta' = \left((\partial y/\partial\phi)(\partial x/\partial\lambda) - (\partial x/\partial\phi)(\partial y/\partial\lambda)\right)/(R^2 hk\cos\phi). The angular distortion is equal to \lvert\pi/2 - \theta'\rvert.

tissot

Using the Nelder–Mead method, I numerically minimized the total angular distortion in 461 approximately equally spaced points located on land while keeping the borders of the map at least 20 kilometres from any land. With respect to these conditions, the optimal parameters of the projection turned out to be:

  • The coordinates of the northern vertex: 73.10° N 120.96° E.
  • The rotation of the tetrahedron: −23.56°.
  • The coordinates of the vertices: b = 0.763 a, c = 0.963 a.
  • The vertical scale factor: 1.03.

These parameters yield a rectangle whose sides are in ratio 4:1.96. To separate Antarctica from Africa, I cut the map vertically at x = 1.6.

With a forward projection, you can make vector maps. For making raster maps, you need the inverse projection. As inverting the equations of the projection algebraically is a hopeless task, I resorted to the numerical Broyden’s method. The initial guess of (\phi, \lambda) = f^{-1}(x, y) comes from inverting the gnomonic projection. When the accuracy is set to 10-7 R, i.e. about one metre, the method usually finds \phi and \lambda after 5 iterations.

The map below is based on a public domain raster image from Natural Earth. The other maps in this post use public domain shapefiles from the same source.

earth

Advertisements
Reproducing the AuthaGraph World Map