AA Square Equal-area Map Projection
MATTHEW A. PETROFF,
Johns Hopkins UniversityA novel square equal-area map projection is proposed. The projection combines closed-formforward and inverse solutions with relatively low angular distortion and minimal cusps, a com-bination of properties not manifested by any previously published square equal-area projection.Thus, the new projection has lower angular distortion than any previously published squareequal-area projection with a closed-form solution. Utilizing a quincuncial arrangement, the newprojection places the north pole at the center of the square and divides the south pole betweenits four corners; the projection can be seamlessly tiled. The existence of closed-form solutionsmakes the projection suitable for real-time visualization applications, both in cartography andin other areas, such as for the display of panoramic images.CCS Concepts: •
Human-centered computing → Geographic visualization ; •
Computingmethodologies → Image processing.Additional Key Words and Phrases: map projection, world map, equal-area, quincuncial
Although there are a plenitude of map projections [21], there has been relatively littlework done on square equal-area projections. As noted in Gringorten [9], a squareaspect ratio is useful for printed atlases, since it allows for a standard page to beefficiently filled with a map along with a title and descriptive caption. Equal-areaprojections are useful for displaying climatology and population data on the Earth,as well as for displaying cosmology data on the celestial sphere. Mappings from thesphere to the square are also useful in computer graphics applications, such as thedisplay of panoramic images, since textures used by graphic processing units (GPUs)are often square. Quincuncial arrangements, which have previously been proposed fordisplaying panoramic images in an aesthetically-pleasing manner [7], can additionallybe seamlessly tiled, which is a desirable property for texture interpolation. Equal-areaprojections allow for efficient use of resources, since all areas of the sphere have aconstant pixel density.The most notable existing square equal-area projection is that of Gringorten [9],which is a quincuncial projection that was derived using differential equations, al-though other arrangements have also been proposed [4]. It has relatively low angulardistortion and has a smooth derivative, and thus no cusps, everywhere but the equator.However, it is complicated and does not have a closed-form solution, which makes itill-suited for use in real-time graphics display. The other significant existing squareequal-area projection is what will be referred to here as the
Collignon quincuncial projection, which consists of an interrupted Collignon projection [5] for each oc-tant of a spherical octahedron displayed in a quincuncial arrangement. It is calledthe equal-area zenithial orthotriangular projection in Huang et al. [12], the octahe-dral equal area partition in Yan et al. [27], and the triangular octahedral equal area projection in McGlynn et al. [16]. It is briefly described as a possible variant of theCollignon projection, with neither a specific name nor a figure, in Snyder [22] and
Author’s address: Matthew A. Petroff, Johns Hopkins University, Department of Physics & Astronomy,3400 N Charles St, Baltimore, Maryland, 21218, [email protected]. a r X i v : . [ c s . G R ] A ug Matthew A. Petroff is described with neither a specific name, a figure, nor a reference to the Collignonprojection in Maurer [15] (translated to English in Warntz [26]); it is also describedfor a single hemisphere in Roşca [20] and Holhoş and Roşca [11]. This projectionhas the advantage of being mathematically simple and has a closed-form solution,for both forward and inverse mappings. However, it has larger angular distortionsthan the Gringorten projection and has significant visible cusps along the edges ofeach projected octahedron face. Both of these projections, as with all quincuncialprojections, can be thought of as polyhedral map projections for the octahedron withan additional linear rescaling along the axis perpendicular to the equator. Squareequal-area projections with non-quincuncial arrangments, e.g., square aspect ratiocylindrical equal-area projections [1, 24], will not be considered, since they have fewerinterruptions and thus much higher angular distortion.Thus, to improve upon existing square equal-area projections, a new projectionmust combine the low angular distortion of the Gringorten projection and its reducedoccurrence of significant visible cusps with the Collignon quincuncial projection’s ad-vantage of having closed-form solutions. These are the properties of the new projectionpresented in this manuscript. The vertex-oriented great circle projection “slice-and-dice” technique of van Leeuwen and Strebe [25] is generalized for the octahedronsuch that the subdivision of each face can be optimized to minimize angular distortionwhen each octant of a spherical octahedron is mapped to an isosceles right triangleinstead of to an equilateral triangle, while still maintaining the equal-area propertyof the mapping; these isosceles right triangles are then displayed in a quincuncialarrangement. This projection is described in detail in Section 2, the properties of theprojection are discussed in Section 3, and the paper concludes in Section 4.
Here, we wish to create a mapping from latitude, ϕ ∈ [− π / , π / ] , and longitude, λ ∈ [ , π ) , to planar coordinates, x , y ∈ [− , ] . The creation of an equal-area mappingbased on an octahedron can be decomposed into the mapping of one octant of aspherical octahedron onto an equilateral triangle. Furthermore, the octant can besplit into two symmetric right triangles, so only a hexadecant, one-sixteenth of thesphere, has to be considered. This right triangle can be further decomposed into threesmaller triangles that share a vertex along the edge by which the right triangle ismirrored. This decomposition is shown in Figure 1, which serves as a reference forangle, distance, and vertex labels used later in this derivation. In van Leeuwen andStrebe [25], this vertex was chosen to be at the center of the spherical octahedronoctant, which makes the smaller triangles right triangles that are identical up to amirroring about the hypotenuse. While this is optimal when the spherical octant ismapped to an equilateral triangle, it leads to increased average angular distortion whenthis equilateral triangle is scaled to become a right isosceles triangle, as is necessary toform a square map. Thus, the more general case is considered here, where the vertexposition is allowed to move along the octant’s mirroring axis.In the more general case, the smaller triangles differ in area and are not all righttriangles when mapped to the plane. For the map to remain equal-area, the area ratiosof the smaller triangles to the hexadecant must be constant between the sphericaltriangles and the Euclidean triangles. This is accomplished by setting the height of Square Equal-area Map Projection 3 ( r , θ ) A B Cϕ D E π / ( r , θ ) ξ ψ ψ ψ A B C h D E √ Fig. 1. Octant mapping. The octant of the spherical octahedron is shown on the left, and thecorresponding equilateral Euclidean triangle is shown on the right. Angles, distances, andvertices used in the text are labeled. The right halves are highlighted, since the projection ismirrored about the center line. triangle B ′ C ′ D ′ such that the area ratio between triangle B ′ C ′ D ′ and triangle A ′ B ′ C ′ matches the area ratio between spherical triangle BCD and spherical triangle
ABC .This produces the relation h ′ √ / √ / = arcsin (cid:16) / (cid:112) − cos ϕ (cid:17) + arcsin (cid:16) ϕ / (cid:112) − cos ( ϕ ) (cid:17) − π / π / , (1)which results in h ′ = π (cid:34) arcsin (cid:32) (cid:112) − cos ϕ (cid:33) + arcsin (cid:32) ϕ (cid:112) − cos ( ϕ ) (cid:33) − π (cid:35) (2)when solved for h ′ , where ϕ is the latitude of the dividing point on the sphericaltriangle and h ′ is the corresponding distance on the Euclidean triangle. Next, the angle ξ ′ is adjusted such that the area ratio between triangle C ′ D ′ E ′ and triangle A ′ B ′ C ′ matches the area ratio between spherical triangle CDE and spherical triangle
ABC .This produces the relation (cid:0) h ′ + (cid:1) sin (cid:104) π / − arctan (cid:16) h ′ /√ (cid:17)(cid:105) sin (cid:104) π / − arctan (cid:16) h ′ /√ (cid:17) + ξ (cid:105) /( ξ ) √ / = π − (cid:16) / (cid:112) − cos ϕ (cid:17) − arcsin (cid:16) ϕ / (cid:112) − cos ( ϕ ) (cid:17) π / , (3)which results in ξ ′ = arctan (cid:104) π ( h ′ − ) (cid:46) (cid:16) √ (cid:104) π (cid:0) h ′ − h ′ + (cid:1) −
96 arcsin (cid:16) / (cid:112) − cos ϕ (cid:17) −
48 arcsin (cid:16) ϕ / (cid:112) − cos ( ϕ ) (cid:17)(cid:105) (cid:17)(cid:105) (4)when solved for ξ ′ . Since the area ratios of two of the three smaller triangles arenow preserved, the area ratio of the third smaller triangle is also preserved. With the Matthew A. Petroff γα βϵH FaG bcr xy H F α β γ ϵ a u v c x y r G b Fig. 2. Sub-triangle mapping. The spherical triangle is shown on the left, while the Euclideantriangle is shown on the right. Angles, distances, and vertices used in the text are labeled. relative areas of the three smaller triangles preserved, an equal-area mapping of eachof the smaller triangles will result in a mapping that is equal-area over the wholesphere.Next, the mapping within the sub-triangles needs to be defined. Figure 2, serves asa reference for geometry labels for the sub-triangles. Here, we will restrict ourselvesto positive latitudes by setting ϕ c = | ϕ | and center longitudes around the center ofeach octant with λ c = λ − π . The longitude of the dividing point of each octant is q = (cid:22) π λ (cid:23) (5) λ = π q . (6)We then define each point by an angle around and a distance from the dividing pointwith θ = | arctan2 [ cos ϕ c sin ( λ c − λ ) , sin ϕ cos ϕ c cos ( λ c − λ ) − cos ϕ sin ϕ c ]| (7) r = arccos [ sin ϕ sin ϕ c + cos ϕ cos ϕ c cos ( λ c − λ )] , (8)where θ is the counter-clockwise angle around the dividing point starting at zero atthe equator and r is the distance from the dividing point. We then define the angle, β ,between the ray to the target point and the hypotenuse of the sub-triangle with ψ = arcsin (cid:32) (cid:112) − cos ϕ (cid:33) (9) ψ = π − ψ (10) β = ψ − θ θ ≤ ψ θ − ψ ψ < θ ≤ ψ + ψ π − θ ψ + ψ < θ . (11)The hypotenuse of the sub-triangle, c , is Square Equal-area Map Projection 5 c = (cid:40) arccos (cid:16) cos ϕ /√ (cid:17) θ ≤ ψ + ψ π / − ϕ θ > ψ + ψ . (12)The angle of the spherical sub-triangle at the vertex at the dividing point is G = ψ θ ≤ ψ ψ ψ < θ ≤ ψ + ψ ψ ψ + ψ < θ , (13)and G ′ = ψ ′ θ ≤ ψ ψ ′ ψ < θ ≤ ψ + ψ ψ ′ θ > ψ + ψ (14)is the corresponding angle on the Euclidean sub-triangle, where ψ ′ = arctan (cid:32) √ h ′ (cid:33) (15) ψ ′ = π − ψ ′ − ξ ′ (16) ψ ′ = ξ ′ − π . (17)The angle of the spherical sub-triangle at the corner of the octant is F = arcsin (cid:18) ϕ √ − cos ( ϕ ) (cid:19) θ ≤ ψ π / − arcsin (cid:18) ϕ √ − cos ( ϕ ) (cid:19) ψ < θ ≤ ψ + ψ π / θ > ψ + ψ . (18)The length of the side of the Euclidean sub-triangle extending from the mid-point ofthe octant side to the dividing point is a ′ = (cid:40) h ′ θ ≤ ψ √ h ′ + sin [ π / − arctan ( h ′ /√ )] sin ξ ′ θ > ψ , (19)and c ′ = (cid:40) √ h ′ + θ ≤ ψ + ψ − h ′ θ > ψ + ψ (20)is the hypotenuse length of the same Euclidean sub-triangle. The distance from theoctant corner to the target point on the spherical sub-triangle is x = arccos ( cos r cos c + sin r sin c cos β ) , (21)and the angle between that ray and the hypotenuse of the spherical sub-triangle is γ = arcsin (cid:18) sin β sin r sin x (cid:19) . (22) Matthew A. Petroff
If that ray is extended to meet the side of the spherical triangle that extends from themid-point of the octant side to the dividing point, the angle between that ray and theray extending from the intersection point to the dividing point is ϵ = arccos ( sin G sin γ cos c − cos G cos γ ) . (23)With these angles and lengths defined, we can proceed with the “slice-and-dice”approach. First, we divide the side of the Euclidean triangle that extends from themid-point of the octant side to the dividing point into lengths u ′ and v ′ , and we dividethe line that extends from the octant corner to this side into lengths x ′ and y ′ . Wethen define the ratio of the areas of the two slices of the Euclidean sub-triangle tomatch that of the spherical sub-triangle using u ′ ( u ′ + v ′ ) = γ + G + ϵ − πF + G − π . (24)Next, we “dice” these slices by defining the separation point on the slicing line to alsomaintain the area ratios usingcos ( x + y ) = cos (cid:20) arcsin (cid:18) sin G sin c sin ϵ (cid:19)(cid:21) = (cid:115) − (cid:18) sin G sin c sin ϵ (cid:19) (25) x ′ ( x ′ + y ′ ) = (cid:115) − cos x − cos ( x + y ) . (26)Solving for a few remaining lengths and angles, u ′ = ( a ′ ) (cid:18) u ′ ( u ′ + v ′ ) (cid:19) (27) ( x ′ + y ′ ) = √ u ′ + c ′ − u ′ c ′ cos G ′ (28)cos γ ′ = cos (cid:20) arcsin (cid:18) u ′ sin G ′ ( x ′ + y ′ ) (cid:19)(cid:21) = (cid:115) − (cid:18) u ′ sin G ′ ( x ′ + y ′ ) (cid:19) (29) x ′ = ( x ′ + y ′ ) (cid:18) x ′ ( x ′ + y ′ ) (cid:19) (30) y ′ = ( x ′ + y ′ ) − x ′ , (31)we can define the target point on the Euclidean triangle in polar coordinates withrespect to the dividing point with r ′ = (cid:112) x ′ + c ′ − x ′ c ′ cos γ ′ (32) α ′ = arccos (cid:18) y ′ − u ′ − r ′ − u ′ r ′ (cid:19) . (33) Square Equal-area Map Projection 7
We then convert this angle from coordinates relative to the sub-triangle to coordinatesrelative to the octant, θ ′ = α ′ θ ≤ ψ π − ξ ′ − α ′ ψ < θ ≤ ψ + ψ π − ξ ′ + α ′ ψ + ψ < θ . (34)Then, we convert to Cartesian coordinates with x c = sgn ( λ c − λ ) r ′ sin θ ′ (35) y c = h ′ − r ′ cos θ ′ (36)and separate northern and southern hemispheres with y h = y c sgn ϕ − . (37)Finally, we squish the equilateral triangle into a right isosceles triangle and arrangethe octants into a square with ζ = π + π q (38) (cid:18) x m y m (cid:19) = (cid:18) x c cos ζ − y h sin ζ /√ x c sin ζ + y h cos ζ /√ (cid:19) √ √ . (39)As the ratio of areas is an affine invariant, this operation maintains the equal-areaproperty of the projection.The maximum angular distortion at a given point, ω , is calculated using Tissot’sindicatrix [21, 23]. This is defined using h = (cid:115)(cid:18) ∂ x ∂ ϕ (cid:19) + (cid:18) ∂ y ∂ ϕ (cid:19) (40) k = (cid:115)(cid:18) ∂ x ∂ λ (cid:19) + (cid:18) ∂ y ∂ λ (cid:19) (41)sin η ′ = hk cos ϕ (cid:20)(cid:18) ∂ y ∂ ϕ (cid:19) (cid:18) ∂ x ∂ λ (cid:19) − (cid:18) ∂ x ∂ ϕ (cid:19) (cid:18) ∂ y ∂ λ (cid:19)(cid:21) (42) ω = (cid:32)(cid:115) h + k − hk sin η ′ h + k + hk sin η ′ (cid:33) , (43)where h and k are the scale factors along the meridian and parallel, respectively, and η ′ is the angle at which a given meridian and parallel intersect. In order to reduce thisdistortion, numerical optimization was used to find the value of ϕ that minimized theaverage of ω across the entire map. This result was very nearly 3 π /
8, so ϕ was set tothis exact value. A world map constructed using the projection is shown in Figure 3.Although the above equations are not directly invertible, the vertex-oriented greatcircle projection technique does lend itself to a closed-form inverse mapping, and aclosed-form inverse mapping of the projection just presented is given in Appendix A. Matthew A. Petroff
Fig. 3. A world map constructed using the new projection presented in this work. A graticuleis overlaid on the map. Note the lack of significant visible cusps.
As presented in van Leeuwen and Strebe [25], the vertex-oriented great circle projec-tion technique is symmetric around the center of each triangle it projects. Althoughthis is optimal for mapping equilateral triangles, it inflates the angular distortion forisosceles triangles. The generalization of said technique in this work, while main-taining its equal-area property, allows it to be optimized for isosceles triangles, suchas the isosceles right triangles that subdivide a quincuncial projection. In additionalto its use in the quincuncial projection presented in this work, the generalized tech-nique can also be used to produce optimized equal-area star projections, which can beconstructed by mapping to isosceles triangles.A comparison of the new projection presented in this work with the Collignon quin-cuncial and Gringorten projections is shown in Figure 4, which includes a world mapand an angular distortion map, as calculated using Tissot’s indicatrix [equation (43)],
Square Equal-area Map Projection 9
Collignon quincuncial Gringorten This work A n g u l a r D i s t o r t i o n ( r a d ) Fig. 4. Map projection comparison. The new projection presented in this work is compared tothe existing Collignon quincuncial and Gringorten projections. The top row shows a world map,while the bottom row shows the angular distortion ( ω ), as calculated using Tissot’s indicatrix.Table 1. Angular distortion summary. Projection Average (rad) Standard deviation (rad) Max (rad)Collignon quincuncial 0 .
68 0 .
18 1 . .
52 0 .
31 1 . .
54 0 .
27 0 . As can be seen in this comparison, the new projection containsno significant visible cusps, which improves upon the Gringorten projection’s cuspat the equator and the Collignon quincuncial projection’s significant equatorial andmeridional cusps. All three projections can be seamlessly tiled, with adjacent tilesrotated by 180°. To derive distortion statistics, the distortion was calculated at 10 000points placed using a Fibonacci lattice [2]. The average, standard deviation, and max-imum angular distortion values are given in Table 1 for the three projections. Boththe new projection and the Gringorten projection have significantly lower averageangular distortion than the Collignon quincuncial projection. The angular distortionof the new projection and the Gringorten projection are comparable, with the newprojection having lower maximum distortion but the Gringorten projection havingslightly lower average distortion. As all three projections are equal-area, comparisonof area distortions is moot. The distortion was calculated using analytic differentiation for the Collignon quincuncial projection andthe new projection presented in this work but was calculated using a finite-difference method for theGringorten projection, since it needs to be implemented using iterative methods.
Fig. 5. Alternative uses. On the left, a panoramic image is projected with the new projectionpresented in this work, while on the right, the reprocessed Haslam 408 MHz sky map [10, 19],in Galactic coordinates, is shown using the same projection.
Excluding constants, which can be pre-computed, and eliminating trivial identities,e.g., z = cos ( arccos z ) , the given forward mapping equations of the new projectionrequire 12 trigonometric function evaluations and 6 inverse trigonometric functionevaluations; four of the trigonometric function evaluations can be further eliminatedusing identities for the composition of trigonometric and inverse trigonometric func-tions, bringing the total number of required trigonometric function evaluations downto eight. While more computationally complex than the Collignon quincuncial projec-tion, which requires only a single trigonometric function evaluation in its forwardmapping, the new projection still has closed-form forward and inverse solutions,which allow it to be implemented as a GPU shader for real-time visualizations, unlikethe iterative methods required for the Gringorten projection.Alternative uses of the projection, outside of cartography, are demonstrated inFigure 5. Unlike for world maps, where the outside of an approximate sphere isprojected, full spherical panoramas and sky maps seek to project the view from thecenter of the sphere, an analogous operation. As can be seen in the figure, the projectedpanorama appears without extreme distortions, so the projected image is suitablefor direct viewing. At the same time, the closed-form forward solution allows theprojected image to be used as an input for the real-time rendering of rectilinear views,as is done by panorama viewers [17]. The combination of low angular distortion andarea equivalence makes the projection better suited for panorama applications than thecommonly used equirectangular format, which has both extreme area distortions andextreme angle distortions. Sky maps are similar to panoramas except that they showthe celestial sphere instead of a view of nearby objects. Figure 5 also shows a projectionof the reprocessed Haslam 408 MHz sky map [10, 19], in Galactic coordinates, as anexample of this. This view of the radio sky primarily traces emission by Galacticsynchrotron emission, which is brightest along the Galactic plane, in particular in thedirection of the Galactic center. Square Equal-area Map Projection 11
The new equal-area quincuncial projection presented in this work combines the bene-fits of the Collignon quincuncial projection with that of the Gringorten projection. Ithas closed-form forward and inverse solutions, unlike the Gringorten projection, whilemaintaining relatively low angular distortion when compared to the Collignon quin-cuncial projection; it also lacks any significant visible cusps, which is an improvementupon both existing equal-area quincuncial projections. The fully-vectorizable closed-form solutions make the new projection suitable for use in real-time visualizationapplications, including hardware-accelerated GPU implementations. The generalizedvertex-oriented great circle projection technique used in the projection’s constructioncan also be used to produce other map projections, such as equal-area star projections.The example implementation code used for the analysis in this manuscript has beenmade available [18].
ACKNOWLEDGMENTS
The author acknowledges support from the National Science Foundation Division ofAstronomical Sciences under grant numbers 1636634 and 1654494. The preparationof this work made use of JAX [6] for automatic analytic differentiation in orderto calculate Tissot’s indicatrix and Numba [14] to speed up the calculation of theGringorten projection. The world maps in Figures 3 and 4 are based on the publicdomain Natural Earth small-scale land polygons. Asymptote [3], Matplotlib [13],Healpy [8, 28], and D3.js were used in the preparation of the included figures. Thereprocessed Haslam 408 MHz map data used in Figure 5 were downloaded from theLegacy Archive for Microwave Background Data Analysis (LAMBDA), part of theHigh Energy Astrophysics Science Archive Center (HEASARC); HEASARC/LAMBDAis a service of the Astrophysics Science Division at the NASA Goddard Space FlightCenter. REFERENCES [1] Charles Arden-Close. 1947.
Geographical By-ways and some other Geographical Essays . Edward Arnold& Co., 87–88.[2] Sergio Baselga. 2018. Fibonacci lattices for the evaluation and optimization of map projections.
Computers & Geosciences
117 (Aug. 2018), 1–8. https://doi.org/10.1016/j.cageo.2018.04.012[3] John C. Bowman and Andy Hammerlindl. 2008. Asymptote: A vector graphics language.
TUGboat
Car-tography and Geographic Information Science
29, 4 (Jan. 2002), 381–390. https://doi.org/10.1559/152304002782008341[5] Édouard Collignon. 1865. Recherches sur la Représentation Plane de la Surface du Globe Ter-restre.
Journal de l’École Impériale Polytechnique
24 (1865), 73–163. https://hdl.handle.net/2027/coo.31924069342180[6] Roy Frostig, Matthew Johnson, and Chris Leary. 2018. Compiling machine learning programs viahigh-level tracing. In
SysML 2018 . https://mlsys.org/Conferences/2019/doc/2018/146.pdf[7] Daniel M. German, Pablo d’Angelo, Michael Gross, and Bruno Postle. 2007. New Methods to ProjectPanoramas for Practical and Aesthetic Purposes. In
Computational Aesthetics in Graphics, Visualization, https://d3js.org/ and Imaging , Douglas W. Cunningham, Gary Meyer, and Laszlo Neumann (Eds.). The EurographicsAssociation. https://doi.org/10.2312/COMPAESTH/COMPAESTH07/015-022[8] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann. 2005.HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed onthe Sphere. The Astrophysical Journal
Journal of Applied Meteorology
Astronomy & Astrophysics Supplement Series
47 (Jan. 1982), 1–143.[11] Adrian Holhoş and Daniela Roşca. 2014. An octahedral equal area partition of the sphere and nearoptimal configurations of points.
Computers & Mathematics with Applications
67, 5 (March 2014),1092–1107. https://doi.org/10.1016/j.camwa.2014.01.003[12] Shaobo Huang, Ryosuke Shibasaki, Masahiro Kasuya, and Masataka Takagi. 1998. Comparative studyon spherical tessellation schemes for global GIS.
Geocarto International
13, 1 (March 1998), 3–14.https://doi.org/10.1080/10106049809354623[13] John D. Hunter. 2007. Matplotlib: A 2D Graphics Environment.
Computing in Science & Engineering
LLVM ’15 . ACM Press. https://doi.org/10.1145/2833157.2833162[15] Hans Maurer. 1935.
Ebene Kugelbilder: Ein Linnésches System der Kartenentwürfe . Number 221 inPetermanns Mitteilungen. Justus Perthes.[16] Thomas McGlynn, Jonathan Fay, Curtis Wong, and Philip Rosenfield. 2019. Octahedron-based Pro-jections as Intermediate Representations for Computer Imaging: TOAST, TEA, and More.
The Astro-physical Journal Supplement Series
Journal of OpenSource Software
4, 40 (Aug. 2019), 1628. https://doi.org/10.21105/joss.01628[18] Matthew A. Petroff. 2020. Supplement to
A Square Equal-area Map Projection . https://doi.org/10.5281/zenodo.3986153[19] M. Remazeilles, C. Dickinson, A. J. Banday, M. A. Bigot-Sazy, and T. Ghosh. 2015. An improved source-subtracted and destriped 408-MHz all-sky map.
Monthly Notices of the Royal Astronomical Society
Appl. Math. Comput.
Flattening the Earth: Two Thousand Years of Map Projections . University ofChicago Press, 114.[23] A. Tissot. 1881.
Mémoire sur la Représentation des Surfaces et les Projections des Cartes Géographiques .Gauthier-Villars.[24] Waldo Tobler and Zi-tan Chen. 1986. A Quadtree for Global Information Storage.
GeographicalAnalysis
18, 4 (Oct. 1986), 360–371. https://doi.org/10.1111/j.1538-4632.1986.tb00108.x[25] Diederik van Leeuwen and Daniel Strebe. 2006. A “Slice-and-Dice” Approach to Area Equivalencein Polyhedral Map Projections.
Cartography and Geographic Information Science
33, 4 (Jan. 2006),269–286. https://doi.org/10.1559/152304006779500687[26] William Warntz (Ed.). 1968.
Plane Globe Projection: a Linnean System of Map Projection . Number 23 inHarvard Papers in Theoretical Geography: Geography and the Properties of Surfaces. Laboratory forComputer Graphics and Spatial Analysis, Harvard, 166. https://apps.dtic.mil/sti/citations/AD0675812English translation of Maurer (1935).[27] Jin Yan, Xiao Song, and Guanghong Gong. 2016. Averaged ratio between complementary profiles forevaluating shape distortions of map projections and spherical hierarchical tessellations.
Computers &Geosciences
87 (Feb. 2016), 41–55. https://doi.org/10.1016/j.cageo.2015.11.009[28] Andrea Zonca, Leo Singer, Daniel Lenz, Martin Reinecke, Cyrille Rosset, Eric Hivon, and KrzysztofGorski. 2019. Healpy: equal area pixelization and spherical harmonics transforms for data on the spherein Python.
Journal of Open Source Software
4, 35 (March 2019), 1298. https://doi.org/10.21105/joss.01298
Square Equal-area Map Projection 13
A INVERSE MAPPING
Here, the inverse mapping is derived, using the previously defined ϕ . Although theforward mapping equations are not directly invertible, an inverse mapping can bederived by performing the “slice-and-dice” technique in reverse, by matching arearatios on the spherical sub-triangles to those on the Euclidean sub-triangles. Startingwith projected map coordinates x m , y m ∈ [− , ] , q = x m ≥ , y m ≥ x m > , y m > x m ≤ , y m > x m < , y m ≤ ζ = π + π q (45) (cid:18) x c y c (cid:19) = (cid:18) √ ( x m cos ζ + y m sin ζ ) √ ( y m cos ζ − x m sin ζ ) (cid:19) (46)gives Cartesian coordinates relative to the Euclidean equilateral triangle face center.We then separate the northern and southern hemispheres with ϕ sgn = (cid:40) y c ≥ − − yc < − y h = (cid:40) y c y c ≥ − − − y c y c < − . (48)The corresponding polar coordinates are then r ′ = (cid:113) x c + ( h ′ − − y h ) (49) θ ′ = | arctan2 ( x c , h ′ − − y h )| , (50)where h ′ is defined by equation (2). Using ψ ′ , ψ ′ , and ψ ′ defined by equations (15),(16), (17), respectively, we find α ′ = θ ′ θ ′ ≤ ψ ′ π − ψ ′ − θ ′ ψ ′ < θ ′ ≤ ψ ′ + ψ ′ θ ′ + ψ ′ − π ψ ′ + ψ ′ < θ ′ . (51)The fixed sub-triangle lengths c , G , G ′ , F , a ′ , and c ′ are defined by equations (12), (13),(14), (18), (19), and (20), respectively, except with θ replaced with θ ′ , ψ replaced with ψ ′ , and ψ replaced with ψ ′ . The remaining necessary fixed sub-triangle length isdefined with b = π / θ ′ ≤ ψ ′ arctan (cid:16) √ ϕ (cid:17) ψ ′ < θ ′ ≤ ψ ′ + ψ ′ π / − arctan (cid:16) √ ϕ (cid:17) ψ ′ + ψ ′ < θ ′ . (52) Other angles and lengths of the Euclidean sub-triangles can then be found with β ′ = G ′ − α ′ (53) x ′ = (cid:112) r ′ + c ′ − r ′ c ′ cos β ′ (54) γ ′ = arccos (cid:18) r ′ − x ′ − c ′ − x ′ c ′ (cid:19) (55) ϵ ′ = π − G ′ − γ ′ (56) y ′ = r ′ sin α ′ sin ϵ ′ (57) u ′ = (cid:112) c ′ + ( x ′ + y ′ ) − c ′ ( x ′ + y ′ ) cos γ ′ (58) v ′ = a ′ − u ′ . (59)Next, the “slice” operation is performed in reverse to match the area ratio of thespherical sub-triangles to that of the Euclidean sub-triangles, with v ′ a ′ = δ + arccos ( sin δ cos b ) − π / F + G − π / δ = arctan (cid:18) − sin [ v ′ ( F + G − π / )/ a ′ ] cos b − cos [ v ′ ( F + G − π / )/ a ′ ] (cid:19) (61)when solved for δ . Other parameters of the spherical sub-triangles are then foundwith γ = F − δ (62)cos ( x + y ) = cos (cid:20) arctan (cid:18) tan b cos δ (cid:19)(cid:21) = (cid:44)(cid:115) + (cid:18) tan b cos δ (cid:19) . (63)Next, the “dice” step is reversed, with equation (26) yielding x = arccos (cid:34) − (cid:18) x ′ ( x ′ + y ′ ) (cid:19) [ − cos ( x + y )] (cid:35) (64)when solved for x . A distance from and an azimuth around the spherical trianglecenter point are then found with r = arccos ( cos x cos c + sin x sin c cos γ ) (65) β = arcsin (cid:18) sin x sin γ sin r (cid:19) (66) α = ψ − β θ ′ ≤ ψ ′ β − ψ ψ ′ < θ ′ ≤ ψ ′ + ψ ′ π − β ψ ′ + ψ ′ < θ ′ . (67)Finally, these are converted to latitude, ϕ , and longitude, λ , with ϕ h = arcsin ( sin ϕ cos r − cos ϕ sin r cos α ) (68) λ = π + π q (69) Square Equal-area Map Projection 15 ϕ = ϕ sgn ϕ h (70) λ = λ + sgn y c arctan (cid:18) sin α sin r cos ϕ cos r − sin ϕ sin ϕ h (cid:19) ..