5 Spherical map projections

In this section we present the transformation equations for all spherical map projections likely to be of use in astronomy. Many of these such as the gnomonic, orthographic, zenithal equidistant, Sanson-Flamsteed, Hammer-Aitoff and COBE quadrilateralized spherical cube are in common use. Others with special properties such as the stereographic, Mercator, and the various equal area projections could not be excluded. A selection of the conic and polyconic projections, much favored by cartographers for their minimal distortion properties, has also been included. However, we have omitted numerous other projections which we considered of mathematical interest only. Evenden (1991) presents maps of the Earth for 73 different projections, although without mathematical definition, including most of those described here. These are particularly useful in judging the distortion introduced by the various projections. Snyder (1993) provides fascinating background material on the subject; historical footnotes in this paper, mainly highlighting astronomical connections, are generally taken from this source. It should be evident from the wide variety of projections described here that new projections could readily be accommodated, the main difficulty being in obtaining general recognition for them from the FITS community.

Cartographers have often given different names to special cases of a class of projection. This applies particularly to oblique projections which, as we have seen in Sect. 2, the current formalism handles in a general way. While we have tended to avoid such special cases, the gnomonic, stereographic, and orthographic projections, being specializations of the zenithal perspective projection, are included for conformance with the AIPS convention. It is also true that zenithal and cylindrical projections may be thought of as special cases of conic projections (see Sect. 5.4). However, the limiting forms of the conic equations tend to become intractable and infinite-valued projection parameters may be involved. Even when the conic equations don't have singularities in these limits it is still likely to be less efficient to use them than the simpler special-case equations. Moreover, we felt that it would be unwise to disguise the true nature of simple projections by implementing them as special cases of more general ones. In the same vein, the cylindrical equal area projection, being a specialization of the cylindrical perspective projection, stands on its own right, as does the Sanson-Flamsteed projection which is a limiting case of Bonne's projection. A list of aliases is provided in Appendix A, Table A.1.

Figure 3: (Left) native coordinate system with its pole at the reference point, i.e. $(\phi _0,\theta _0) = (0,90\hbox {$^\circ $ })$, and (right) with the intersection of the equator and prime meridian at the reference point i.e. $(\phi _0,\theta _0) = (0,0)$.

The choice of a projection often depends on particular special properties that it may have. Certain equal area projections (also known as authalic, equiareal, equivalent, homalographic, homolographic, or homeotheric) have the property that equal areas on the sphere are projected as equal areas in the plane of projection. This is obviously a useful property when surface density must be preserved. Mathematically, a projection is equiareal if and only if the Jacobian,

$$\frac{\partial(x,y)}{\partial(\phi\cos\theta,\theta)} \equiv \frac{1}{\cos\theta} \left\vert \begin{array}{l@{\hskip 10pt}l} \... ...{{\textstyle \partial y}}{{\textstyle \partial\theta}} \end{array}\right\vert ,$$

is constant.

Conformality is a property which applies to points in the plane of projection which are locally distortion-free. Practically speaking, this means that the projected meridian and parallel through the point intersect at right angles and are equiscaled. A projection is said to be conformal or orthomorphic if it has this property at all points. Such a projection cannot be equiareal. Conformal projections preserve angle; the angle of intersection of two lines on the sphere is equal to that of their projection. It must be stressed that conformality is a local property, finite regions in conformal projections may be very severely distorted.

A projection is said to be equidistant if the meridians are uniformly, truely, or correctly divided so that the parallels are equispaced. That is, the native latitude is proportional to the distance along the meridian measured from the equator, though the constant of proportionality may differ for different meridians. Equidistance is not a fundamental property. It's main benefit is in facilitating measurement from the graticule since linear interpolation may be used over the whole length of the meridian. This is especially so if the meridians are projected as straight lines which is the case for all equidistant projections presented here.

Zenithal, or azimuthal projections, discussed in Sect. 5.1, give the true azimuth to all points on the map from the reference point at the native pole. By contrast, retroazimuthal projections give the true azimuth from all points on the map to the reference point, measured as an angle from the vertical. The first projection specifically designed with this property, Craig's "Mecca'' projection of 1909, allowed Muslim worshippers to find the direction to Mecca for daily prayers. Such maps have also been used to allow radio operators to determine the bearing to radio transmitters. In practice, however, retroazimuthal projections may be considered mathematical curiosities of questionable value; most are so severely distorted as to be difficult to read, and we have not included any in this work. Instead the stereographic projection (Sect. 5.1.4) can serve the same purpose, except that the azimuth to the reference point must be measured with respect to the, typically curved, inclined meridians, rather than from the vertical.

A number of projections have other special properties and these will be noted for each.

Maps of the Earth are conventionally displayed with terrestrial latitude increasing upwards and longitude to the right, i.e. north up and east to the right, as befits a sphere seen from the outside. On the other hand, since the celestial sphere is seen from the inside, north is conventionally up and east to the left. The AIPS convention arranged that celestial coordinates at points near the reference point should be calculable to first order via the original linear prescription of Wells et al. (1981), i.e. $(\alpha,\delta)\approx(\alpha_0,\delta_0)+(x,y)$. Consequently, the CDELT ia keyword value associated with the right ascension was negative while that for the declination was positive. The handedness of the (x,y) coordinates as calculated by the AIPS convention equivalent of Eq. (1) is therefore opposite to that of the (p₁,p₂) pixel coordinates.

Figure 4: (Left) geometry of the zenithal perspective projections, the point of projection at P is $\mu$ spherical radii from the center of the sphere; (right) the three important special cases. This diagram is at 3:2 scale compared to the graticules of Figs. 8, 9 and 10.

In accordance with the image display convention of Paper I we think of the p₁-pixel coordinate increasing to the right with p₂ increasing upwards, i.e. a right-handed system. This means that the (x,y) coordinates must be left-handed as shown in Fig. 3. Note, however, that the approximation $(\alpha,\delta)\approx(\alpha_0,\delta_0)+(x,y)$ cannot hold unless 1) $(\phi _0,\theta _0)$, and hence $(\alpha _0,\delta _0)$, do actually map to the reference point (Sect. 2.2), 2) $\phi_{\rm p}$assumes its default value (Sect. 2.8), and 3) the projection is scaled true at the reference point (some are not as discussed in Sect. 7.2). Figure 3 also illustrates the orientation of the native coordinate system with respect to the (x,y)coordinate system for the two main cases.

Cartographers, for example Kellaway (1946), think of spherical projections as being a projection of the surface of a sphere onto a plane, this being the forward direction; the deprojection from plane back to sphere is thus the inverse or reverse direction. However, this is at variance with common usage in FITS where the transformation from pixel coordinates to world coordinates is considered the forward direction. We take the cartographic view in this section as being the natural one and trust that any potential ambiguity may readily be resolved by context.

The requirement stated in Sect. 1 that (x,y) coordinates in the plane of projection be measured in "degrees'' begs clarification. Spherical projections are usually defined mathematically in terms of a scale factor, r₀, known as the "radius of the generating sphere''. However, in this work r₀ is set explicitly to $180\hbox{$^\circ$ }/ \pi$ in order to maintain backwards compatibility with the AIPS convention. This effectively sets the circumference of the generating sphere to $360\hbox{$^\circ$ }$ so that arc length is measured naturally in degrees (rather than radians as for a unit sphere). However, this true angular measure on the generating sphere becomes distorted when the sphere is projected onto the plane of projection. So while the "degree'' units of r₀ are notionally carried over by conventional dimensional analysis to the (x,y) they no longer represent a true angle except near the reference point (for most projections).

In addition to the (x,y) coordinates, the native spherical coordinates, $(\phi,\theta)$, celestial coordinates, $(\alpha ,\delta )$, and all other angles in this paper are measured in degrees. In the equations given below, the arguments to all trigonometric functions are in degrees and all inverse trigonometric functions return their result in degrees. Whenever a conversion between radians and degrees is required it is shown explicitly. All of the graticules presented in this section have been drawn to the same scale in (x,y) coordinates in order to represent accurately the exaggeration and foreshortening found in some projections. It will also be apparent that since FITS image planes are rectangular and the boundaries of many projections are curved, there may sometimes be cases when the FITS image must contain pixels that are outside the boundary of the projection. These pixels should be blanked correctly and geometry routines should return a sensible error code to indicate that their celestial coordinates are undefined.

   
5.1 Zenithal (azimuthal) projections

Zenithal or azimuthal projections all map the sphere directly onto a plane. The native coordinate system is chosen to have the polar axis orthogonal to the plane of projection at the reference point as shown on the left side of Fig. 3. Meridians of native longitude are projected as uniformly spaced rays emanating from the reference point and the parallels of native latitude are mapped as concentric circles centered on the same point. Since all zenithal projections are constructed with the pole of the native coordinate system at the reference point we set

$$(\phi_0, \theta_0)_{\rm zenithal} = (0,90\hbox{$^\circ$ }) .$$ (11)

Zenithal projections are completely specified by defining the radius as a function of native latitude, $R_{\theta}$. Rectangular Cartesian coordinates, (x,y), in the plane of projection as defined by Eq. (1), are then given by

  

$$\hphantom{-}R_{\theta}\sin\phi ,$$ (12)

$$-R_{\theta}\cos\phi ,$$ (13)

which may be inverted as

  

$$\phi = \arg~(-y, x) ,$$ (14)

$$R_{\theta}= \sqrt{x^2 + y^2} .$$ (15)

   
5.1.1 AZP: Zenithal perspective

Figure 5: Alternate geometries of slant zenithal perspective projections with $\mu = 2$ and $\gamma = 30\hbox {$^\circ $ }$: (left) tilted plane of projection ( AZP), (right) displaced point of projection P ( SZP). Grey lines in each diagram indicate the other point of view. They differ geometrically only by a scale factor, the effect of translating the plane of projection. Each projection has its native pole at the reference point, r and r', but these are geometrically different points. Thus the native latitudes $\theta$ and $\theta '$ of the geometrically equivalent points s and s' differ. Consequently the two sets of projection equations have a different form. This diagram is at 3:2 scale compared to the graticules of Figs. 6 and 7.

Zenithal (azimuthal) perspective projections are generated from a point and carried through the sphere to the plane of projection as illustrated in Fig. 4. We need consider only the case where the plane of projection is tangent to the sphere at its pole; the projection is simply rescaled if the plane intersects some other parallel of native latitude. If the source of the projection is at a distance $\mu$ spherical radii from the center of the sphere with $\mu$ increasing in the direction away from the plane of projection, then it is straightforward to show that

$$R_{\theta}= \frac{180\hbox{$^\circ$ }}{\pi}\frac{(\mu+1) \cos\theta}{\mu+\sin\theta} ,$$ (16)

with $\mu \neq -1$ being the only restriction. When $R_{\theta}$ is given Eq. (16) has two solutions for $\theta$, one for each side of the sphere. The following form of the inverse equation always gives the planeward solution for any $\mu$

$$\theta = \arg~(\rho, 1) - \sin^{-1} \left( \frac{\rho\mu}{\sqrt{\rho^2 + 1}} \right) ,$$ (17)

where

$$\rho = \frac{\pi}{180\hbox{$^\circ$ }}\frac{R_{\theta}}{\mu + 1} \cdot$$ (18)

For $\vert\mu\vert \neq 1$ the sphere is divided by a native parallel at latitude $\theta_{\rm x}$ into two unequal segments that are projected in superposition:

$$\theta_{\rm x} = \left\{ \begin{array}{ll} \sin^{-1}(-1/\mu)... ...u) & \mbox{\ldots $\vert\mu\vert < 1$ } \end{array} \right. .$$ (19)

For $\vert\mu\vert > 1$, the projection is bounded and both segments are projected in the correct orientation, whereas for $\vert\mu\vert \leq 1$ the projection is unbounded and the anti-planewards segment is inverted.

A near-sided perspective projection may be obtained with $\mu < -1$. This correctly represents the image of a sphere, such as a planet, when viewed from a distance $\vert\mu\vert$ times the planetary radius. The coordinates of the reference point may be expressed in planetary longitude and latitude, $(\lambda,\beta)$. Also, the signs of the relevant CDELT ia may be chosen so that longitude increases as appropriate for a sphere seen from the outside rather than from within.

It is particularly with regard to planetary mapping that we must generalize AZP to the case where the plane of projection is tilted with respect to the axis of the generating sphere, as shown on the left side of Fig. 5. It can be shown (Sect. 7.4.1) that this geometry is appropriate for spacecraft imaging with non-zero look-angle, $\gamma$, the angle between the camera's optical axis and the line to the center of the planet.

Such slant zenithal perspective projections are not radially symmetric and their projection equations must be expressed directly in terms of x and y:

  

$$\hphantom{-}R \sin\phi ,$$ (20)

$$- R \sec\gamma \cos\phi ,$$ (21)

where

$$R = \frac{180\hbox{$^\circ$ }}{\pi}\frac{(\mu+1) \cos\theta} {(\mu+\sin\theta) + \cos\theta \cos\phi \tan\gamma} ~$$ (22)

is a function of $\phi$ as well as $\theta$. Equations (20) and (21) reduce to the non-slant case for $\gamma = 0$. The inverse equations are

$$\phi \arg~(-y \cos\gamma, x) ,$$ (23)

$$\theta \left\{ \begin{array}{ll} \psi - \omega \\ \psi + \omega + 180\hbox{$^\circ$ } \end{array} \right. ,$$ (24)

where

$$\psi \arg~(\rho, 1) ,$$ (25)

$$\omega \sin^{-1} \left( \frac{\rho\mu}{\sqrt{\rho^2 + 1}} \right) ,$$ (26)

$$\rho \frac{R}{\frac{180\hbox{$^\circ$ }}{\pi}(\mu + 1) + y \sin\gamma} ,$$ (27)

$$\sqrt{x^2 + y^2\cos^2\gamma} .$$ (28)

For $\vert\mu\vert < 1$ only one of the solutions for $\theta$ will be valid, i.e. lie in the range $[-90\hbox{$^\circ$ }, 90\hbox{$^\circ$ }]$ after normalization. Otherwise there will be two valid solutions; the one closest to $90\hbox{$^\circ$ }$ should be chosen.

Figure 6: Slant zenithal perspective ( AZP) projection with $\mu = 2$ and $\gamma = 30\hbox {$^\circ $ }$ which corresponds to the left-hand side of Fig. 5. The reference point of the corresponding SZP projection is marked at $(\phi ,\theta ) = (0,60\hbox {$^\circ $ })$.

With $\gamma \neq 0$ the projection is not scaled true at the reference point. In fact the x scale is correct but the y scale is magnified by $\sec\gamma$, thus stretching parallels of latitude near the pole into ellipses (see Fig. 6). This also shows the native meridians projected as rays emanating from the pole. For constant $\theta$, each parallel of native latitude defines a cone with apex at the point of projection. This cone intersects the tilted plane of projection in a conic section. Equations (20) and (21) reduce to the parametric equations of an ellipse, parabola, or hyperbola; the quantity

$$C = (\mu + \sin\theta)^2 - \tan^2\gamma \cos^2\theta$$ (29)

determines which:

$$\begin{array}{llll} C & > & 0 & {\rm ...ellipse}, \\ C & =... ......parabola}, \\ C & < & 0 & {\rm ...hyperbola}. \end{array}$$ (30)

If $\vert\mu \cos\gamma\vert \le 1$ then the open conic sections are possible and C = 0 when

$$\theta = \pm \gamma - \sin^{-1}(\mu \cos\gamma) .$$ (31)

For C > 0 the eccentricity of the ellipse is a function of $\theta$, as is the offset of its center in y.

Definition of the perimeter of the projection is more complicated for the slant projection than the orthogonal case. As before, for $\vert\mu\vert > 1$ the sphere is divided into two unequal segments that are projected in superposition. The boundary between these two segments is what would be seen as the limb of the planet in spacecraft photography. It corresponds to native latitude

$$\theta_{\rm x} = \sin^{-1}(-1/\mu) ,$$ (32)

which is projected as an ellipse, parabola, or hyperbola for $\vert\mu \cos\gamma\vert$ greater than, equal to, or less than 1 respectively.

In general, for $\vert\mu \cos\gamma\vert > 1$, the projection is bounded, otherwise it is unbounded. However, the latitude of divergence is now a function of $\phi$:

$$\theta_\infty = \left\{ \begin{array}{ll} \psi - \omega \\ \psi + \omega + 180\hbox{$^\circ$ } \end{array} \right. ,$$ (33)

where

$$\psi - \tan^{-1}( ~ \tan\gamma \cos\phi ) ,$$ (34)

$$\omega \hphantom{-}\sin^{-1} \left( \frac{\mu}{\sqrt{1 + \tan^2\gamma \cos^2\phi}} \right) \cdot$$ (35)

Zero, one, or both of the values of $\theta_\infty$ given by Eq. (33) may be valid, i.e. lie in the range $[-90\hbox{$^\circ$ }, 90\hbox{$^\circ$ }]$ after normalization.

The FITS keywords PV i_1a and PV i_2a, attached to latitude coordinate i, will be used to specify, respectively, $\mu$ in spherical radii and $\gamma$ in degrees, both with default value 0.

   
5.1.2 SZP: Slant zenithal perspective

While the generalization of the AZP projection to tilted planes of projection is useful for certain applications it does have a number of drawbacks, in particular, unequal scaling at the reference point.

Figure 7: Slant zenithal perspective ( SZP) projection with $\mu = 2$ and $(\phi _{\rm c},\theta _{\rm c}) = (180\hbox {$^\circ $ },60\hbox {$^\circ $ })$, which corresponds to the right-hand side of Fig. 5. The reference point of the corresponding AZP projection is marked at $(\phi ,\theta ) = (180\hbox {$^\circ $ },60\hbox {$^\circ $ })$.

Figure 5 shows that moving the point of projection, P, off the axis of the generating sphere is equivalent, to within a scale factor, to tilting the plane of projection. However this approach has the advantage that the plane of projection remains tangent to the sphere. Thus the projection is conformal at the native pole as can be seen by the circle around the native pole in Fig. 7. It is also quite straightforward to formulate the projection equations with P offset in x as well as y.

It is interesting to note that this slant zenithal perspective (SZP) projection also handles the case that corresponds to $\gamma = 90\hbox{$^\circ$ }$ in AZP. AZP fails in this extreme since P falls in the plane of projection - effectively a scale factor of zero is applied to AZP over the corresponding SZP case. One of the more important aspects of SZP is the application of its limiting case with $\mu = \infty$ in aperture synthesis radio astronomy as discussed in Sect. 5.1.5. One minor disadvantage is that the native meridians are projected as curved conic sections rather than straight lines.

If the Cartesian coordinates of P measured in units of the spherical radius are $(x_{\rm p},y_{\rm p},z_{\rm p})$, then

  

$$\hphantom{-}\frac{180\hbox{$^\circ$ }}{\pi}\frac{z_{\rm p} \cos\theta \sin\phi - x_{\rm p} (1 - \sin\theta)} {z_{\rm p} - (1 - \sin\theta)} ,$$ (36)

$$- \frac{180\hbox{$^\circ$ }}{\pi}\frac{z_{\rm p} \cos\theta \cos\phi + y_{\rm p} (1 - \sin\theta)} {z_{\rm p} - (1 - \sin\theta)} \cdot$$ (37)

To invert these equations, compute $\theta$ first via

 

$$\theta \sin^{-1} \left( \frac{-b \pm \sqrt{b^2 - ac}}{a} \right) ,$$ (38)

 

a = X'² + Y'² + 1 , (39)

b = X'(X-X') + Y'(Y-Y') , (40)

c = (X - X')² + (Y - Y')² - 1 , (41)

$$(X, Y) = \frac{\pi}{180\hbox{$^\circ$ }}(x, y) ,$$ (42)

$$(X',Y') = (X - x_{\rm p}, Y - y_{\rm p}) / z_{\rm p} .$$ (43)

Choose $\theta$ closer to $90\hbox{$^\circ$ }$ if Eq. (38) has two valid solutions; then $\phi$ is given by

$$\phi = \arg~(-(Y - Y'(1 - \sin\theta)), ~ X - X'(1 - \sin\theta)) .$$ (44)

The Cartesian coordinates of P are simply related to the parameters used in AZP. If $\mu$ is the distance of P from the center of the sphere O, and the line through P and O intersects the sphere at $(\phi_{\rm c},\theta_{\rm c})$ on the planewards side (point r in Fig. 5, right), then $\gamma = 90\hbox{$^\circ$ }- \theta_{\rm c}$ and

$$x_{\rm p} - \mu \cos\theta_{\rm c} \sin\phi_{\rm c} ,$$ (45)

$$y_{\rm p} \hphantom{-}\mu \cos\theta_{\rm c} \cos\phi_{\rm c} ,$$ (46)

$$z_{\rm p} \hphantom{-}\mu \sin\theta_{\rm c} + 1 ,$$ (47)

where $\mu$ is positive if P lies on the opposite side of O from r and negative otherwise. For a non-degenerate projection we require $z_{\rm p} \neq 0$ and this is the only constraint on the projection parameters.

For $\vert\mu\vert > 1$ the sphere is divided into two unequal segments that are projected in superposition. The limb is defined by computing the native latitude $\theta_{\rm x}$ as a function of $\phi$

$$\theta_{\rm x} = \left\{ \begin{array}{ll} \psi - \omega \\ \psi + \omega + 180 \end{array} \right. ,$$ (48)

where

$$\psi = \arg~(\rho, \sigma) ,$$ (49)

$$\omega = \sin^{-1} \left( \frac{1}{\sqrt{\rho^2 + \sigma^2}} \right) ,$$ (50)

$$(\rho, \sigma) = (z_{\rm p} - 1, ~ x_{\rm p} \sin\phi - y_{\rm p} \cos\phi) ,$$ (51)

$$\hspace*{9mm} = (\mu \sin\theta_{\rm c}, ~ -\mu \cos\theta_{\rm c} \cos(\phi - \phi_{\rm c}) .$$ (52)

Zero, one, or both of the values of $\theta_{\rm x}$ given by Eq. (48) may be valid, i.e. lie in the range $[-90\hbox{$^\circ$ }, 90\hbox{$^\circ$ }]$ after normalization. A second boundary constraint applies if $\vert 1 - z_{\rm p}\vert \le 1$ in which case the projection diverges at native latitude:

$$\theta_\infty = \sin^{-1}(1 - z_{\rm p}) .$$ (53)

The FITS keywords PV i_1a, PV i_2a, and PV i_3a, attached to latitude coordinate i, will be used to specify, respectively, $\mu$ in spherical radii with default value 0, $\phi_{\rm c}$ in degrees with default value 0, and $\theta_{\rm c}$ in degrees with default value $90\hbox{$^\circ$ }$.

   
5.1.3 TAN: Gnomonic

Figure 8: Gnomonic ( TAN) projection; diverges at $\theta = 0$.

The zenithal perspective projection with $\mu = 0$, the gnomonic projection, is widely used in optical astronomy and was given its own code within the AIPS convention, namely TAN. For $\mu = 0$, Eq. (16) reduces to

$$R_{\theta}= \frac{180\hbox{$^\circ$ }}{\pi}\cot \theta ,$$ (54)

with inverse

 

$$\theta \tan^{-1} \left( \frac{180\hbox{$^\circ$ }}{\pi R_{\theta}} \right) \cdot$$ (55)

The gnomonic projection is illustrated in Fig. 8. Since the projection is from the center of the sphere, all great circles are projected as straight lines. Thus, the shortest distance between two points on the sphere is represented as a straight line interval, which, however, is not uniformly divided. The gnomonic projection diverges at $\theta = 0$, but one may use a gnomonic projection onto the six faces of a cube to display the whole sky. See Sect. 5.6.1 for details.

   
5.1.4 STG: Stereographic

Figure 9: Stereographic ( STG) projection; diverges at $\theta = -90\hbox {$^\circ $ }$.

The stereographic projection, the second important special case of the zenithal perspective projection defined by the AIPS convention, has $\mu = 1$, for which Eq. (16) becomes

$$R_{\theta} \frac{180\hbox{$^\circ$ }}{\pi}\frac{2\cos\theta}{1 + \sin\theta} ,$$ (56)

$${\frac{360\hbox{$^\circ$ }}{\pi}} \tan \left( \frac{90\hbox{$^\circ$ }- \theta}{2} \right) ,$$

with inverse

$$\theta = 90\hbox{$^\circ$ }- 2 \tan^{-1} \left( \frac{\pi R_{\theta}}{360\hbox{$^\circ$ }} \right) \cdot$$ (57)

The stereographic projection illustrated in Fig. 9 is the conformal (orthomorphic) zenithal projection, everywhere satisfying the isoscaling requirement

$$\frac{\partial R_\theta}{\partial \theta} = \frac{-\pi R_\theta}{180\hbox{$^\circ$ }\cos\theta} \cdot$$ (58)

This allows its use as a replacement for retroazimuthal projections, as discussed in Sect. 5.

The stereographic projection also has the amazing property that it maps all circles on the sphere to circles in the plane of projection, although concentric circles on the sphere are not necessarily concentric in the plane of projection. This property made it the projection of choice for Arab astronomers in constructing astrolabes. In more recent times it has been used by the Astrogeology Center for maps of the Moon, Mars, and Mercury containing craters, basins, and other circular features.

   
5.1.5 SIN: Slant orthographic

The zenithal perspective projection with $\mu = \infty$, the orthographic projection, is illustrated in the upper portion of Fig. 10 (at consistent scale). It represents the visual appearance of a sphere, e.g. a planet, when seen from a great distance.

Figure 10: Slant orthographic ( SIN) projection: (top) the orthographic projection, $(\xi ,\eta ) = (0,0)$, north and south sides begin to overlap at $\theta = 0$; (bottom left) $(\phi _{\rm c},\theta _{\rm c}) = (225\hbox {$^\circ $ },60\hbox {$^\circ $ })$, i.e. $(\xi ,\eta ) = (-1/\sqrt {6}, 1/\sqrt {6})$; (bottom right) projection appropriate for an east-west array observing at $\delta _0 = 60\hbox {$^\circ $ }$, $(\phi _{\rm c},\theta _{\rm c}) = (180\hbox {$^\circ $ },60\hbox {$^\circ $ })$, $(\xi ,\eta ) = (0, 1/\sqrt {3})$.

The orthographic projection is widely used in aperture synthesis radio astronomy and was given its own code within the AIPS convention, namely SIN. Use of this projection code obviates the need to specify an infinite value as a parameter of AZP. In this case, Eq. (16) becomes

$$R_{\theta}= \frac{180\hbox{$^\circ$ }}{\pi}\cos\theta , \\$$ (59)

with inverse

$$\theta = \cos^{-1} \left( \frac{\pi}{180\hbox{$^\circ$ }}R_{\theta}\right) \cdot$$ (60)

In fact, use of the orthographic projection in radio interferometry is an approximation, applicable only for small field sizes. However, an exact solution exists where the interferometer baselines are co-planar. It reduces to what Greisen (1983) called the NCP projection for the particular case of an East-West interferometer (Brouw 1971). The projection equations (derived in Appendix C) are

  

$$\hphantom{-}\frac{180\hbox{$^\circ$ }}{\pi}\left[~ \cos\theta \sin\phi + \xi ~ (1 - \sin\theta) \right] ,$$ (61)

$$- \frac{180\hbox{$^\circ$ }}{\pi}\left[~ \cos\theta \cos\phi - \eta ~ (1 - \sin\theta) \right] .$$ (62)

These are the equations of the "slant orthographic'' projection, equivalent to Eqs. (36) and (37) of the SZP projection in the limit $\mu = \infty$, with

$$\xi \hphantom{-}\cot\theta_{\rm c} \sin\phi_{\rm c} ,$$ (63)

$$\eta - \cot\theta_{\rm c} \cos\phi_{\rm c} .$$ (64)

It can be shown that the slant orthographic projection is equivalent to an orthographic projection centered at $(x,y) = \frac{180\hbox{$^\circ$ }}{\pi}(\xi, \eta)$ which has been stretched in the $\phi_{\rm c}$ direction by a factor of ${\rm cosec}~\theta_{\rm c}$. The projection equations may be inverted using Eqs. (38) and (44) except that Eq. (43) is replaced with

$$(X', Y') = (\xi, \eta) .$$ (65)

The outer boundary of the SIN projection is given by Eq. (48) in the limit $\mu = \infty$:

$$\theta_{\rm x} = - \tan^{-1} \left( \xi\sin\phi - \eta\cos\phi \right) .$$ (66)

Two example graticules are illustrated in the lower portion of Fig. 10. We here extend the original SIN projection of the AIPS convention to encompass the slant orthographic projection, with the dimensionless $\xi$ and $\eta$ given by keywords PV i_1a and PV i_2a, respectively, attached to latitude coordinate i, both with default value 0.

   
5.1.6 ARC: Zenithal equidistant

Figure 11: Zenithal equidistant ( ARC) projection; no limits.

Some non-perspective zenithal projections are also of interest in astronomy. The zenithal equidistant projection first appeared in Greisen (1983) as ARC. It is widely used as the approximate projection of Schmidt telescopes. As illustrated in Fig. 11, the native meridians are uniformly divided to give equispaced parallels. Thus

 

$$R_{\theta} 90\hbox{$^\circ$ }- \theta ,$$ = (67)

which is trivially invertible. This projection was also known in antiquity.

   
5.1.7 ZPN: Zenithal polynomial

Figure 12: Zenithal polynomial projection ( ZPN) with parameters, 0.050, 0.975, -0.807, 0.337, -0.065, 0.010, 0.003, -0.001; limits depend upon the parameters.

The zenithal polynomial projection, ZPN, generalizes the ARC projection by adding polynomial terms up to a large degree in the zenith distance. We define it as

$$R_{\theta}= \frac{180\hbox{$^\circ$ }}{\pi}\sum_{m=0}^{20} P_... ...80\hbox{$^\circ$ }}(90\hbox{$^\circ$ }- \theta) \right)^m\cdot$$ (68)

Note the dimensionless units of P_m imparted by $\pi/180\hbox{$^\circ$ }$. Allowance is made for a polynomial of degree up to 20 as a conservative upper limit that should encompass all practical applications. For speed and numerical precision the polynomial should be evaluated in Horner form, i.e.

$$(\ldots(P_{20} \gamma + P_{19})\gamma + \ldots P_2)\gamma + P_1)\gamma + P_0$$

where $\gamma = (180\hbox{$^\circ$ }/\pi)(90\hbox{$^\circ$ }- \theta)$.

Since its inverse cannot be expressed analytically, ZPN should only be used when the geometry of the observations require it. In particular, it should never be used as an nth-degree expansion of one of the standard zenithal projections.

If P₀ is non-zero the native pole is mapped to an open circle centered on the reference point as illustrated in Fig. 12. In other words, $(\phi _0,\theta _0) = (0,90\hbox {$^\circ $ })$ is not at (x,y) = (0,0), which in fact lies outside the boundary of the projection. However, we do not dismiss $P_0 \neq 0$ as a possibility since it is not inconsistent with the formalism presented in Sect. 2.2 and could conceivably be useful for images which do not contain the reference point. Needless to say, care should be exercised in constructing and interpreting such systems particularly in that $(\alpha _0,\delta _0)$ (i.e. the CRVAL ia) do not specify the celestial coordinates of the reference point (the CRPIX ja).

P_m (dimensionless) is given by the keywords PV i_0a, PV i_1a, $\ldots$, PV i_20a, attached to latitude coordinate i, all of which have default values of zero.

   
5.1.8 ZEA: Zenithal equal-area

Figure 13: Zenithal equal area projection ( ZEA); no limits.

Lambert's zenithal equal-area projection illustrated in Fig. 13 is constructed by defining $R_{\theta}$ so that the area enclosed by the native parallel at latitude $\theta$ in the plane of projection is equal to the area of the corresponding spherical cap. It may be generated using

 

$$R_{\theta} \frac{180\hbox{$^\circ$ }}{\pi}\sqrt{2 (1 - \sin\theta)}$$ (69)

$${\frac{360\hbox{$^\circ$ }}{\pi}} \sin \left( \frac{90\hbox{$^\circ$ }- \theta}{2} \right) ,$$

with inverse

$$\theta = 90\hbox{$^\circ$ }- 2 \sin^{-1} \left(\frac{\pi R_{\theta}}{360\hbox{$^\circ$ }}\right) \cdot$$ (70)

   
5.1.9 AIR: Airy projection

The Airy projection minimizes the error for the region within latitude $\theta_{\rm b}$ (Evenden 1991). It is defined by

$$R_{\theta}= -2\frac{180\hbox{$^\circ$ }}{\pi}\left( \frac{\l... ...rac{\ln(\cos\xi_{\rm b})}{\tan^2\xi_{\rm b}}\tan\xi \right) ,$$ (71)

where

$$\xi \frac{90\hbox{$^\circ$ }- \theta}{2} ,$$

$$\xi_{\rm b} \frac{90\hbox{$^\circ$ }- \theta_{\rm b}}{2} \cdot$$

When $\theta_{\rm b}$ approaches $90\hbox{$^\circ$ }$, the second term of Eq. (71) approaches its asymptotic value of $-\frac{1}{2}$. For all $\theta_{\rm b}$, the projection is unbounded at the native south pole. Inversion of Eq. (71), a transcendental equation in $\theta$, must be done via iterative methods.

Figure 14: Airy projection ( AIR) with $\theta _{\rm b} = 45\hbox {$^\circ $ }$; diverges at $\theta = -90\hbox {$^\circ $ }$.

The FITS keyword PV i_1a, attached to latitude coordinate i, will be used to specify $\theta_{\rm b}$ in degrees with a default of  $90\hbox{$^\circ$ }$. This projection is illustrated in Fig. 14.

   
5.2 Cylindrical projections

Figure 15: Geometry of cylindrical projections.

Cylindrical projections are so named because the surface of projection is a cylinder. The native coordinate system is chosen to have its polar axis coincident with the axis of the cylinder. Meridians and parallels are mapped onto a rectangular graticule so that cylindrical projections are described by formulæ which return x and y directly. Since all cylindrical projections are constructed with the native coordinate system origin at the reference point, we set

$$(\phi_0, \theta_0)_{\rm cylindrical} = (0,0) .$$ (72)

Furthermore, all cylindrical projections have

$$x \propto \phi .$$ (73)

Cylindrical projections are often chosen to map the regions adjacent to a great circle, usually the equator, with minimal distortion.

   
5.2.1 CYP: Cylindrical perspective

Figure 15 illustrates the geometry for the construction of cylindrical perspective projections. The sphere is projected onto a cylinder of radius $\lambda$ spherical radii from points in the equatorial plane of the native system at a distance $\mu$ spherical radii measured from the center of the sphere in the direction opposite the projected surface. The cylinder intersects the sphere at latitudes $\theta_{\rm x} = \cos^{-1}\lambda$. It is straightforward to show that

$$\lambda \phi ,$$ (74)

$$\frac{180\hbox{$^\circ$ }}{\pi}\left( \frac{\mu+\lambda}{\mu+\cos\theta}\right) \sin\theta .$$ (75)

This may be inverted as

$$\phi \frac{x}{\lambda} ,$$ (76)

$$\theta \arg~(1,\eta) + \sin^{-1}\left( \frac{\eta\mu}{\sqrt{\eta^2+1}} \right) ,$$ (77)

where

$$\eta = \frac{\pi}{180\hbox{$^\circ$ }}\frac{y}{\mu+\lambda} \cdot$$ (78)

Note that all values of $\mu$ are allowable except $\mu = -\lambda$. For FITS purposes, we define the keywords PV i_1a to convey $\mu$ and PV i_2a for $\lambda$, both measured in spherical radii, both with default value 1, and both attached to latitude coordinate i.

The case with $\mu = \infty$ is covered by the class of cylindrical equal area projections. No other special-cases need be defined since cylindrical perspective projections have not previously been used in FITS. Aliases for a number of special cases are listed in Appendix A, Table A.1. Probably the most important of these is Gall's stereographic projection, which minimizes distortions in the equatorial regions. It has $\mu = 1, \lambda = \sqrt{2}/2$, giving

$$\begin{eqnarray*}x & = & \phi / \sqrt{2} , \\ y & = & \frac{180\hbox{$^\circ$ ... ...\sqrt{2}}{2} \right) \tan \left( \frac{\theta}{2} \right) \cdot \end{eqnarray*}$$

It is illustrated in Fig. 16.

Figure 16: Gall's stereographic projection, CYP with $\mu = 1$, $\theta _{\rm x} = 45\hbox {$^\circ $ }$; no limits.

   
5.2.2 CEA: Cylindrical equal area

The cylindrical equal area projection is the special case of the cylindrical perspective projection with $\mu = \infty$. It is conformal at latitudes $\pm\theta_{\rm c}$ where $\lambda = \cos^2\theta_{\rm c}$. The formulæ are

$$\phi ,$$ (79)

$$\frac{180\hbox{$^\circ$ }}{\pi}\frac{\sin\theta}{\lambda} ,$$ (80)

which reverse to

$$\phi$$ (81)

$$\theta \sin^{-1} \left( \frac{\pi}{180\hbox{$^\circ$ }}\lambda y \right) \cdot$$ (82)

Note that the scaling parameter $\lambda$ is applied to the y-coordinate rather than the x-coordinate as in CYP. The keyword PV i_1a attached to latitude coordinate i is used to specify the dimensionless $\lambda$ with default value 1.

Figure 17: Lambert's equal area projection, CEA with $\lambda = 1$; no limits.

Lambert's equal area projection, the case with $\lambda = 1$, is illustrated in Fig. 17. It shows the extreme compression of the parallels of latitude at the poles typical of all cylindrical equal area projections.

   
5.2.3 CAR: Plate carrée

The equator and all meridians are correctly scaled in the plate carrée projection, whose main virtue is that of simplicity. Its formulæ are

  

$$\phi ,$$ (83)

$$\theta .$$ (84)

The projection is illustrated in Fig. 18.

Figure 18: The plate carrée projection ( CAR); no limits.

   
5.2.4 MER: Mercator

Since the meridians and parallels of all cylindrical projections intersect at right angles the requirement for conformality reduces to that of equiscaling at each point. This is expressed by the differential equation

$$\frac{\partial y}{\partial\theta} = \frac{-1}{\cos\theta} \frac{\partial x}{\partial\phi} ,$$ (85)

the solution of which gives us Mercator's projection:

$$\phi ,$$ (86)

$$\frac{180\hbox{$^\circ$ }}{\pi}\ln \tan \left( \frac{90\hbox{$^\circ$ }+\theta}{2} \right) ,$$ (87)

with inverse

$$\theta = 2 \tan^{-1} \left( {\rm e}^{y\pi/180\hbox{$^\circ$ }} \right) - 90\hbox{$^\circ$ }.$$ (88)

This projection, illustrated in Fig. 19, has been widely used in navigation since it has the property that lines of constant bearing (known as rhumb lines or loxodromes) are projected as straight lines. This is a direct result of its conformality and the fact that its meridians do not converge.

Refer to Sect. 6.1.4 for a discussion of the usage of MER in AIPS.

   
5.3 Pseudocylindrical and related projections

Pseudocylindricals are like cylindrical projections except that the parallels of latitude are projected at diminishing lengths towards the polar regions in order to reduce lateral distortion there. Consequently the meridians are curved. Pseudocylindrical projections lend themselves to the construction of interrupted projections in terrestrial cartography. However, this technique is unlikely to be of use in celestial mapping and is not considered here. Like ordinary cylindrical projections, the pseudocylindricals are constructed with the native coordinate system origin at the reference point. Accordingly we set

$$(\phi_0, \theta_0)_{\rm pseudocylindrical} = (0, 0) .$$ (89)

The Hammer-Aitoff projection is a modified zenithal projection, not a pseudocylindrical, but is presented with this group on account of its superficial resemblance to them.

Figure 19: Mercator's projection ( MER); diverges at $\theta = \pm 90\hbox {$^\circ $ }$.

   
5.3.1 SFL: Sanson-Flamsteed

Bonne's projection (Sect. 5.5.1) reduces to the pseudocylindrical Sanson-Flamsteed projection when $\theta_1 = 0$. Parallels are equispaced and projected at their true length which makes it an equal area projection. The formulæ are

$$\phi \cos\theta ,$$ (90)

$$\theta ,$$ (91)

which reverse into

$$\phi \frac{x}{\cos y} ,$$ (92)

$$\theta$$ (93)

This projection is illustrated in Fig. 20. Refer to Sect. 6.1.4 for a discussion relating SFL to the GLS projection in AIPS.

Figure 20: Sanson-Flamsteed projection ( SFL); no limits.

Figure 21: Parabolic projection ( PAR); no limits.

Figure 22: Mollweide's projection ( MOL); no limits.

Figure 23: Hammer-Aitoff projection ( AIT); no limits.

   
5.3.2 PAR: Parabolic

The parabolic or Craster pseudocylindrical projection is illustrated in Fig. 21. The meridians are projected as parabolic arcs which intersect the poles and correctly divide the equator, and the parallels of latitude are spaced so as to make it an equal area projection. The formulæ are

$$\phi \left( 2\cos \frac{2\theta}{3} - 1 \right) ,$$ (94)

$$180\hbox{$^\circ$ }\sin\frac{\theta}{3} ,$$ (95)

with inverse

$$\theta 3 \sin^{-1} \left( \frac{y}{180\hbox{$^\circ$ }} \right) ,$$ (96)

$$\phi \frac{180\hbox{$^\circ$ }}{\pi}\frac{x}{1 - 4 (y/180\hbox{$^\circ$ })^2} \cdot$$ (97)

   
5.3.3 MOL: Mollweide's

In Mollweide's pseudocylindrical projection, the meridians are projected as ellipses that correctly divide the equator and the parallels are spaced so as to make the projection equal area. The formulæ are

$$\frac{2\sqrt{2}}{\pi} \phi \cos\gamma ,$$ (98)

$$\sqrt{2}\frac{180\hbox{$^\circ$ }}{\pi}\sin\gamma ,$$ (99)

where $\gamma$ is defined as the solution of the transcendental equation

$$\sin\theta = \frac{\gamma}{90\hbox{$^\circ$ }} + \frac{\sin 2\gamma}{\pi} \cdot$$ (100)

The inverse equations may be written directly without reference to $\gamma$,

$$\phi \pi x / \left(2\sqrt{2-(\frac{\pi}{180\hbox{$^\circ$ }}y)^2} ~\right) ,$$ (101)

$$\theta \sin^{-1}\left(\frac{1}{90\hbox{$^\circ$ }}\sin^{-1} \left( \frac{\pi}{180\hbox{$^\circ$ }}\frac{y}{\sqrt{2}}\right) \right.$$

$$\left.+ \frac{y}{180\hbox{$^\circ$ }} \sqrt{2 - (\frac{\pi}{180\hbox{$^\circ$ }}y)^2} ~~\right) \cdot$$ (102)

Mollweide's projection is illustrated in Fig. 22.

   
5.3.4 AIT: Hammer-Aitoff

The Hammer-Aitoff projection illustrated in Fig. 23 is developed from the equatorial case of the zenithal equal area projection by doubling the equatorial scale and longitude coverage. The whole sphere is mapped thereby while preserving the equal area property. Note, however, that the equator is not evenly divided.

This projection reduces distortion in the polar regions compared to pseudocylindricals by making the meridians and parallels more nearly orthogonal. Together with its equal area property this makes it one of most commonly used all-sky projections.

Figure 24: Construction of conic perspective projections ( COP) and the resulting graticules; (left) two-standard projection with $\theta _1 = 20\hbox {$^\circ $ }$, $\theta _2 = 70\hbox {$^\circ $ }$; (right) one-standard projection with $\theta _1 = \theta _2 = 45\hbox {$^\circ $ }$. Both projections have $\theta _{\rm a} = 45\hbox {$^\circ $ }$ and this accounts for their similarity. Both diverge at $\theta = \theta _{\rm a} \pm 90\hbox {$^\circ $ }$.

The formulæ for the projection and its inverse are derived in Greisen (1986) and Calabretta (1992) among others. They are

$$2 \gamma \cos\theta \sin\frac{\phi}{2} ,$$ (103)

$$\gamma \sin \theta ,$$ (104)

where

$$\gamma = \frac{180\hbox{$^\circ$ }}{\pi}\sqrt{\frac{2}{1 + \cos\theta\cos(\phi/2)}} \cdot$$ (105)

The reverse equations are

$$\phi 2 \arg \left(2Z^2 - 1, \frac{\pi}{180\hbox{$^\circ$ }}\frac{Z}{2} x\right) ,$$ (106)

$$\theta \sin^{-1} \left(\frac{\pi}{180\hbox{$^\circ$ }}yZ\right) ,$$ (107)

where

$$\sqrt{1 - \left(\frac{\pi}{180\hbox{$^\circ$ }}\frac{x}{4}\right)^2 - \left(\frac{\pi}{180\hbox{$^\circ$ }}\frac{y}{2}\right)^2} ,$$ (108)

$$\sqrt{\frac{1}{2} \left( 1 + \cos\theta \cos\frac{\phi}{2} \right)} \cdot$$ (109)

Note that $\frac{1}{2} \le Z^2 \le 1$. Refer to Sect. 6.1.4 for a discussion of the usage of AIT in AIPS.

   
5.4 Conic projections

In conic projections the sphere is thought to be projected onto the surface of a cone which is then opened out. The native coordinate system is chosen so that the poles are coincident with the axis of the cone. Native meridians are then projected as uniformly spaced rays that intersect at a point (either directly or by extrapolation), and parallels are projected as equiangular arcs of concentric circles.

Two-standard conic projections are characterized by two latitudes, $\theta_1$ and $\theta_2$, whose parallels are projected at their true length. In the conic perspective projection these are the latitudes at which the cone intersects the sphere. One-standard conic projections have $\theta_1 = \theta_2$ and the cone is tangent to the sphere as shown in Fig. 24. Since conics are designed to minimize distortion in the regions between the two standard parallels they are constructed so that the point on the prime meridian mid-way between the two standard parallels maps to the reference point so we set

$$(\phi_0, \theta_0)_{\rm conic} = (0, \theta_{\rm a}) , \\$$ (110)

where

$$\theta_{\rm a} = (\theta_1 + \theta_2) / 2 .$$ (111)

Being concentric, the parallels may be described by $R_{\theta}$, the radius for latitude $\theta$, and $A_{\phi}$, the angle for longitude $\phi$. An offset in y is also required to force (x,y) = (0,0) at $(\phi, \theta) = (0, \theta_{\rm a})$. All one- and two-standard conics have

$$A_{\phi} C \phi ,$$ = (112)

where C, a constant known as the constant of the cone, is such that the apical angle of the projected cone is $360\hbox{$^\circ$ }~C$. Since the standard parallels are projected as concentric arcs at their true length we have

 

$$C = \frac{180\hbox{$^\circ$ }\cos\theta_1}{\pi R_{\theta_1}} = \frac{180\hbox{$^\circ$ }\cos\theta_2}{\pi R_{\theta_2}} \cdot$$ (113)

Cartesian coordinates in the plane of projection are

$$\hphantom{-}R_{\theta}\sin(C\phi) ,$$ (114)

$$-R_{\theta}\cos(C\phi) + Y_0 ,$$ (115)

and these may be inverted as

  

$$R_{\theta}= \hbox{sign~}\theta_{\rm a} \sqrt{x^2+(Y_0-y)^2} ,$$ (116)

$$\phi = \arg \left( \frac{Y_0-y}{R_{\theta}} , \frac{x}{R_{\theta}} \right) / C .$$ (117)

To complete the inversion the equation for $\theta$ as a function of $R_{\theta}$ is given for each projection. The equations given here correctly invert southern conics, i.e. those with $\theta_{\rm a} < 0$. An example is shown in Fig. 25.

Figure 25: Conic equal area projection ( COE) with $\theta _1 = -20\hbox {$^\circ $ }$, and $\theta _2 = -70\hbox {$^\circ $ }$, also illustrating the inversion of southern hemisphere conics; no limits.

The conics will be parameterized in FITS by $\theta_{\rm a}$ (given by Eq. (111)) and $\eta$ where

$$\eta = \vert \theta_1 - \theta_2 \vert / 2 .$$ (118)

The keywords PV i_1a and PV i_2a attached to latitude coordinate i will be used to specify $\theta_{\rm a}$ and $\eta$respectively, both measured in degrees. PV i_1a has no default while PV i_2a defaults to 0. It is recommended that both keywords always be given. Where $\theta_1$ and $\theta_2$ are required in the equations below they may be computed via

  

$$\theta_1 \theta_{\rm a} - \eta ,$$ (119)

$$\theta_2 \theta_{\rm a} + \eta .$$ (120)

This sets $\theta_2$ to the larger of the two. The order, however, is unimportant.

As noted in Sect. 5, the zenithal projections are special cases of the conics with $\theta_1 = \theta_2 = 90\hbox{$^\circ$ }$. Likewise, the cylindrical projections are conics with $\theta_1 = -\theta_2$. However, we strongly advise against using conics in these cases for the reasons given previously. Nevertheless, the only formal requirement on $\theta_1$ and $\theta_2$ in the equations presented below is that $-90\hbox{$^\circ$ }\leq \theta_1,\theta_2 \leq 90\hbox{$^\circ$ }$.

   
5.4.1 COP: Conic perspective

Development of Colles' conic perspective projection is shown in Fig. 24. The projection formulæ are

$$\sin \theta_{\rm a} ,$$ (121)

$$R_{\theta} \frac{180\hbox{$^\circ$ }}{\pi}\cos\eta~ \left[ \cot\theta_{\rm a} - \tan(\theta-\theta_{\rm a}) \right] ,$$ (122)

$$\frac{180\hbox{$^\circ$ }}{\pi}\cos \eta \cot \theta_{\rm a} .$$ (123)

The inverse may be computed with

$$\theta = \theta_{\rm a} + \tan^{-1}\left( \cot\theta_{\rm a} ... ...{180\hbox{$^\circ$ }}\frac{R_{\theta}}{\cos\eta} \right) \cdot$$ (124)

   
5.4.2 COE: Conic equal area

The standard parallels in Alber's conic equal area projection are projected as concentric arcs at their true length and separated so that the area between them is the same as the corresponding area on the sphere. The other parallels are then drawn as concentric arcs spaced so as to preserve the area. The projection formulæ are

  

$$\gamma/2 ,$$ (125)

$$R_{\theta} \frac{180\hbox{$^\circ$ }}{\pi}\frac{2}{\gamma} \sqrt{1 + \sin\theta_1\sin\theta_2 - \gamma\sin\theta} ,$$ (126)

$$\frac{180\hbox{$^\circ$ }}{\pi}\frac{2}{\gamma} \sqrt{1 + \sin\theta_1\sin\theta_2 - \gamma\sin((\theta_1+\theta_2)/2)} ,$$ (127)

where

$$\gamma = \sin\theta_1 + \sin\theta_2 .$$ (128)

The inverse may be computed with

 

$$\theta \sin^{-1} \left( \frac{1}{\gamma} + \frac{\sin\theta_1\sin\theta_... ...gamma \left( \frac{\pi R_{\theta}}{360\hbox{$^\circ$ }} \right)^2 \right) \cdot$$ (129)

This projection is illustrated in Fig. 25.

   
5.4.3 COD: Conic equidistant

Figure 26: Conic equidistant projection ( COD) with $\theta _1 = 20\hbox {$^\circ $ }$ and $\theta _2 = 70\hbox {$^\circ $ }$; no limits.

In the conic equidistant projection the standard parallels are projected at their true length and at their true separation. The other parallels are then drawn as concentric arcs spaced at their true distance from the standard parallels. The projection formulæ are

$$\frac{180\hbox{$^\circ$ }}{\pi}\frac{\sin\theta_{\rm a}\sin\eta}{\eta} ,$$ (130)

$$R_{\theta} \theta_{\rm a} - \theta + \eta \cot\eta\cot\theta_{\rm a} ,$$ (131)

$$\eta \cot\eta\cot\theta_{\rm a} .$$ (132)

The inverse may be computed with

$$\theta = \theta_{\rm a} + \eta\cot\eta\cot\theta_{\rm a} - R_{\theta}.$$ (133)

For $\theta_1 = \theta_2$ these expressions reduce to

$$C = \sin\theta_{\rm a} ,$$ (134)

$$R_{\theta}= \theta_{\rm a} - \theta + \frac{180\hbox{$^\circ$ }}{\pi}\cot\theta_{\rm a} ,$$ (135)

$$Y_0 = \frac{180\hbox{$^\circ$ }}{\pi}\cot\theta_{\rm a} ,$$ (136)

and inverse

$$\theta = \theta_{\rm a} + \frac{180\hbox{$^\circ$ }}{\pi}\cot\theta_{\rm a} - R_{\theta}~ .$$ (137)

This projection is illustrated in Fig. 26.

   
5.4.4 COO: Conic orthomorphic

The requirement for conformality of conic projections is

$$\frac{\partial R_\theta}{\partial \theta} = \frac{-\pi R_\theta}{180\hbox{$^\circ$ }\cos\theta} C .$$ (138)

Solution of this differential equation gives rise to the formulæ for Lambert's conic orthomorphic projection:

$$\frac{\ln\left(\frac{\cos\theta_2}{\cos\theta_1}\right)} {\ln\lef... ...}{2}\right)}{\tan\left(\frac{90\hbox{$^\circ$ } -\theta_1}{2}\right)}\right]} ,$$ (139)

$$R_{\theta} \psi \left[ \tan\left(\frac{90\hbox{$^\circ$ }-\theta}{2} \right) \right] ^C ,$$ (140)

$$\psi \left[ \tan\left(\frac{90\hbox{$^\circ$ }-\theta_{\rm a}}{2} \right) \right] ^C ,$$ (141)

where

$$\psi \frac{180\hbox{$^\circ$ }}{\pi}\frac{\cos\theta_1}{C\left[ \tan\left( \frac{90\hbox{$^\circ$ }-\theta_1}{2} \right) \right] ^C} ,$$ (142)

$$\frac{180\hbox{$^\circ$ }}{\pi}\frac{\cos\theta_2}{C\left[ \tan\left( \frac{90\hbox{$^\circ$ }-\theta_2}{2} \right) \right] ^C} \cdot$$ (143)

The inverse may be computed with

$$\theta = 90\hbox{$^\circ$ }- 2 \tan^{-1} \left( \left[ \frac{R_{\theta}}{\psi} \right]^\frac{1}{C} \right) \cdot$$ (144)

When $\theta_1 = \theta_2$ the expression for C may be replaced with $C = \sin\theta_1$. This projection is illustrated in Fig. 27.

Figure 27: Conic orthomorphic projection ( COO) with $\theta _1 = 20\hbox {$^\circ $ }$ and $\theta _2 = 70\hbox {$^\circ $ }$; diverges at $\theta = -90\hbox {$^\circ $ }$.

   
5.5 Polyconic and pseudoconic projections

Polyconics are generalizations of the standard conic projections; the parallels of latitude are projected as circular arcs which may or may not be concentric, and meridians are curved rather than straight as in the standard conics. Pseudoconics are a sub-class with concentric parallels. The two polyconics presented here have parallels projected at their true length and use the fact that for a cone tangent to the sphere at latitude $\theta_1$, as shown in Fig. 24, we have $R_{\theta_1} = \frac{180\hbox{$^\circ$ }}{\pi}\cot\theta_1$. Since both are constructed with the native coordinate system origin at the reference point we set

$$(\phi_0, \theta_0)_{\rm polyconic} = (0, 0) .$$ (145)

   
5.5.1 BON: Bonne's equal area

In Bonne's pseudoconic projection all parallels are projected as concentric equidistant arcs of circles of true length and true spacing. This is sufficient to guarantee that it is an equal area projection. It is parameterized by the latitude $\theta_1$ for which $R_{\theta_1} = \frac{180\hbox{$^\circ$ }}{\pi}\cot\theta_1$. The projection is conformal at this latitude and along the central meridian. The equations for Bonne's projection become divergent for $\theta_1 = 0$ and this special case is handled as the Sanson-Flamsteed projection. The projection formulæ are

$$\hphantom{-}R_{\theta}\sin A_{\phi},$$ (146)

$$-R_{\theta}\cos A_{\phi}+ Y_0 ,$$ (147)

where

$$A_{\phi} \frac{180\hbox{$^\circ$ }}{\pi R_{\theta}} \phi\cos\theta ,$$ (148)

$$R_{\theta} Y_0 - \theta ,$$ (149)

$$\frac{180\hbox{$^\circ$ }}{\pi}\cot\theta_1 + \theta_1 .$$ (150)

The inverse formulæ are then

$$\theta Y_0 - R_{\theta},$$ (151)

$$\phi \frac{\pi}{180\hbox{$^\circ$ }}A_{\phi}R_{\theta}~ / \cos\theta ,$$ (152)

where

$$R_{\theta} \hbox{sign~}\theta_1 \sqrt{x^2 + (Y_0-y)^2} ,$$ (153)

$$A_{\phi} \arg \left(\frac{Y_0-y}{R_{\theta}}, \frac{x}{R_{\theta}} \right) \cdot$$ (154)

This projection is illustrated in Fig. 28. The keyword PV i_1a attached to latitude coordinate i will be used to give $\theta_1$in degrees with no default value.

Figure 28: Bonne's projection ( BON) with $\theta _1 = 45\hbox {$^\circ $ }$; no limits.

   
5.5.2 PCO: Polyconic

Each parallel in Hassler's polyconic projection is projected as an arc of a circle of radius $\frac{180\hbox{$^\circ$ }}{\pi}\cot\theta$ at its true length, $360\hbox{$^\circ$ }\cos\theta$, and correctly divided. The scale along the central meridian is true and consequently the parallels are not concentric. The projection formulæ are

$$\frac{180\hbox{$^\circ$ }}{\pi}\cot\theta \sin(\phi\sin\theta) ,$$ (155)

$$\theta + \frac{180\hbox{$^\circ$ }}{\pi}\cot\theta \left[1 - \cos(\phi\sin\theta)\right] .$$ (156)

Inversion requires iterative solution for $\theta$ of the equation

$$x^2 - \frac{360\hbox{$^\circ$ }}{\pi}(y-\theta)\cot\theta + (y-\theta)^2 = 0 .$$ (157)

Once $\theta$ is known $\phi$ is given by

$$\phi = \frac{1}{\sin\theta} \arg \left( \frac{180\hbox{$^\circ$ }}{\pi}- (y-\theta)\tan\theta, x \tan\theta \right) .$$ (158)

The polyconic projection is illustrated in Fig. 29.

   
5.6 Quad-cube projections

Quadrilateralized spherical cube (quad-cube) projections belong to the class of polyhedral projections in which the sphere is projected onto the surface of an enclosing polyhedron, typically one of the five Platonic polyhedra: the tetrahedron, hexahedron (cube), octahedron, dodecahedron, and icosahedron.

The starting point for polyhedral projections is the gnomonic projection of the sphere onto the faces of an enclosing polyhedron. This may then be modified to make the projection equal area or impart other special properties. However, minimal distortion claims often made for polyhedral projections should be balanced against the many interruptions of the flattened polyhedron; the breaks should be counted as extreme distortions. Their importance in astronomy is not so much in visual representation but in solving the problem of distributing N points as uniformly as possible over the sphere. This may be of particular importance in optimizing computationally intensive applications. The general problem may be tackled by pixelizations such as HEALPix (Hierarchical Equal Area isoLatitude PIXelization, Górski & Hivon 1999), which define a distribution of points on the sphere but do not relate these to points in a plane of projection and are therefore outside the scope of this paper.

Figure 29: Polyconic projection ( PCO); no limits.

The quad-cubes have been used extensively in the COBE project and are described by Chan & O'Neill (1975) and O'Neill & Laubscher (1976). The icosahedral case has also been studied by Tegmark (1996). It is close to optimal, providing a 10% improvement over the cubic case. However, we have not included it here since it relies on image pixels being organized in an hexagonal close-packed arrangement rather than the simple rectangular arrangement supported by FITS.

The six faces of quad-cube projections are numbered and laid out as

$$\begin{array}{cccccccc} \hskip 90pt & & & & 0 \\ & 4 & 3 & 2 & 1 & 4 & 3 & 2 \\ & & & & 5 \end{array}$$

where faces 2, 3 and 4 may appear on one side or the other (or both). The layout used depends only on the FITS writer's choice of $\phi_{\rm c}$ in Table 4. It is also permissible to split faces between sides, for example to put half of face 3 to the left and half to the right to create a symmetric layout. FITS readers should have no difficulty determining the layout since the origin of (x,y) coordinates is at the center of face 1 and each face is $90\hbox{$^\circ$ }$ on a side. The range of x therefore determines the layout. Other arrangements are possible and Snyder (1993) illustrates the "saw-tooth'' layout of Reichard's 1803 map of the Earth. While these could conceivably have benefit in celestial mapping we judged the additional complication of representing them in FITS to be unwarranted. The layout used in the COBE project itself has faces 2, 3, and 4 to the left.

 

 

Table 4: Assignment of parametric variables and central longitude and latitude by face number for quadrilateralized spherical cube projections.

Face	$\hphantom{-}\xi$	$\hphantom{-}\eta$	$\hphantom{-}\zeta$		$\phi_{\rm c}$		$\theta_{\rm c}$
0	$\hphantom{-}m$	-l	$\hphantom{-}n$		$0\hbox{$^\circ$ }$		$90\hbox{$^\circ$ }$
1	$\hphantom{-}m$	$\hphantom{-}n$	$\hphantom{-}l$		$0\hbox{$^\circ$ }$		$0\hbox{$^\circ$ }$
2	-l	$\hphantom{-}n$	$\hphantom{-}m$	$\hphantom{-}-270\hbox{$^\circ$ }$	${\rm or}$	$90\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$
3	-m	$\hphantom{-}n$	-l	$-180\hbox{$^\circ$ }$	${\rm or}$	$180\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$
4	$\hphantom{-}l$	$\hphantom{-}n$	-m	$-90\hbox{$^\circ$ }$	${\rm or}$	$270\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$
5	$\hphantom{-}m$	$\hphantom{-}l$	-n		$0\hbox{$^\circ$ }$		$-90\hbox{$^\circ$ }$

The native coordinate system has its pole at the center of face 0 and origin at the center of face 1 (see Fig. 30) which is the reference point whence

$$(\phi_0, \theta_0)_{\rm quad-cube} = (0, 0) .$$ (159)

The face number may be determined from the native longitude and latitude by computing the direction cosines:

 

$$\cos \theta \cos\phi ,$$

$$\cos \theta \sin\phi ,$$ (160)

$$\sin \theta .$$

The face number is that which maximizes the value of $\zeta$ in Table 4. That is, if $\zeta$ is the largest of n, l, m, -l, -m, and -n, then the face number is 0 through 5, respectively. Each face may then be given a Cartesian coordinate system, $(\xi,\eta)$, with origin in the center of each face as per Table 4. The formulæ for quad-cubes are often couched in terms of the variables

$$\chi \xi / \zeta ,$$ (161)

$$\psi \eta / \zeta .$$ (162)

Being composed of six square faces the quad-cubes admit the possibility of efficient data storage in FITS. By stacking them on the third axis of a three-dimensional data structure the storage required for an all-sky map may be halved. This axis will be denoted by a CTYPE ia value of CUBEFACE. In this case the value of (x,y) computed via Eq. (1) for the center of each face must be (0,0) and the FITS interpreter must increment this by $(\phi_{\rm c},\theta_{\rm c})$ using its choice of layout. Since the CUBEFACE axis type is purely a storage mechanism the linear transformation of Eq. (1) must preserve the CUBEFACE axis pixel coordinates.

   
5.6.1 TSC: Tangential spherical cube

While perspective quad-cube projections could be developed by projecting a sphere onto an enclosing cube from any point of projection, inside or outside the sphere, it is clear that only by projecting from the center of the sphere will every face be treated equally. Thus the tangential spherical cube projection (TSC) consists of six faces each of which is a gnomonic projection of a portion of the sphere. As discussed in Sect. 5.1.1, gnomonic projections map great circles as straight lines but unfortunately diverge very rapidly away from the poles and can only represent a portion of the sphere without extreme distortion. The TSC projection partly alleviates this by projecting great circles as piecewise straight lines. To compute the forward projection first determine $\chi$ and $\psi$ as described above, then

$$\phi_{\rm c} + 45\hbox{$^\circ$ }\chi ,$$ (163)

$$\theta_{\rm c} + 45\hbox{$^\circ$ }\psi .$$ (164)

To invert these first determine to which face the (x,y) coordinates refer, then compute

$$\chi (x - \phi_{\rm c}) / 45\hbox{$^\circ$ },$$ (165)

$$\psi (y - \theta_{\rm c}) / 45\hbox{$^\circ$ },$$ (166)

then

$$\zeta 1/\sqrt{1 + \chi^2 + \psi^2} .$$ (167)

Once $\zeta$ is known $\xi$ and $\eta$ are obtained via

$$\xi \chi \zeta ,$$ (168)

$$\eta \psi \zeta .$$ (169)

The direction cosines (l,m,n) may be identified with $(\xi,\eta,\zeta)$ with with the aid of Table 4, whence $(\phi,\theta)$ may readily be computed. The projection is illustrated in Fig. 30 for the full sphere.

   
5.6.2 CSC: COBE quadrilateralized spherical cube

Figure 30: Tangential spherical cube projection ( TSC); no limits.

The COBE quadrilateralized spherical cube projection illustrated in Fig. 31 modifies the tangential spherical cube projection in such a way as to make it approximately equal area. The forward equations are

$$\phi_{\rm c} + 45\hbox{$^\circ$ }~ F (\chi, \psi) ,$$ (170)

$$\theta_{\rm c} + 45\hbox{$^\circ$ }~ F (\psi, \chi) ,$$ (171)

where the function F is given by

 

$$F(\chi,\psi) \chi \gamma^* + \chi^3 (1 - \gamma^*)$$

$$+\chi \psi^2 (1 - \chi^2) \left[ ~ \Gamma + (M - \Gamma) \chi^2 \vphantom{\sum_{i=0}^{\infty}} \right.$$

$$+ \left. (1 - \psi^2) \sum_{i=0}^{\infty} \sum_{j=0}^{\infty} C_{ij} \chi^{2i}\psi^{2j} \right]$$

$$+\chi^3 (1 - \chi^2) \left[ ~ \Omega_1 - (1 - \chi^2) \sum_{i=0}^{\infty} D_i \chi^{2i} \right] .$$ (172)

C_ij and D_i are derived from c_ij^* and d_i^* as given by Chan & O'Neill (1975). The other parameters are given by exact formulæ developed by O'Neill & Laubscher (1976), who provide the numeric values of their parameters in tables and software listings. Both disagree with their formulæ, but the software listings do contain the actual numeric parameters still in use for the COBE Project (Immanuel Freedman, private communication, 1993). They are

$$\begin{array}{lcl} \gamma^* & = & \hphantom{-}1.37484847732 ... ...m{-}0.0759196200467 \\ D_1 & = & -0.0217762490699 \end{array}$$

Chan & O'Neill (1975) actually defined the projection via the inverse equations.

Figure 31: COBE quadrilateralized spherical cube projection ( CSC); no limits.

Their formulation may be rewritten in a more convenient form which is now the current usage in the COBE Project (Immanuel Freedman, private communication, 1993):

$$\chi f(x-\phi_{\rm c}, y-\theta_{\rm c}) ,$$ (173)

$$\psi f(y-\theta_{\rm c}, x-\phi_{\rm c}) ,$$ (174)

where

 

$$f(x-\phi_{\rm c},y-\theta_{\rm c}) = X + X \left(1-X^2\right) \sum_{j=0}^N \sum_{i=0}^{N-j} P_{ij} X^{2i} Y^{2j} ,$$ (175)

and

$$(x - \phi_{\rm c}) / 45\hbox{$^\circ$ },$$

$$(y-\theta_{\rm c})/45\hbox{$^\circ$ } .$$

For COBE, N = 6 and the best-fit parameters have been taken to be

$$\begin{array}{lcrlcr} & P_{00} =& -0.27292696 &\qquad &P_{04... ...P_{13} =& 1.50880086 &\qquad &P_{06} =& 0.14381585. \end{array}$$

Given the face number, $\chi$, and $\psi$, the native coordinates $(\phi,\theta)$ may be computed as for the tangential spherical cube projection.

Equations (172) and (175), the forward and reverse projection equations used by COBE, are not exact inverses. Each set could of course be inverted to any required degree of precision via iterative methods (in that case Eq. (175) should be taken to define the projection). However, the aim here is to describe the projection in use within the COBE project. One may evaluate the closure error in transforming $(x-\phi_{\rm c},y-\theta_{\rm c})$ to $(\chi,\psi)$ with Eq. (175) and then transforming back to $(x-\phi_{\rm c},y-\theta_{\rm c})$ with Eq. (172), i.e.

$$\begin{eqnarray*}E_{ij}^2 & = & \left(F(f(x_i,y_j),f(y_j,x_i)) - x_i\right)^2\\ & & \mbox{}+\left(F(f(y_j,x_i),f(x_i,y_j)) - y_j\right)^2 . \end{eqnarray*}$$

The COBE parameterization produces an average error of 4.7 arcsec over the full field. The root mean square and peak errors are 6.6 and 24 arcsec, respectively. In the central parts of the image ( $\vert X\vert, \vert Y\vert \leq 0.8$), the average and root mean square errors are 5.9 and 8.0 arcsec, larger than for the full field.

Measures of equal-area conformance obtained for Eq. (175) show that the rms deviation is 1.06% over the full face and 0.6% over the inner 64% of the area of each face. The maximum deviation is +13.7% and -4.1% at the edges of the face and only $\pm1.3$% within the inner 64% of the face.

   
5.6.3 QSC: Quadrilateralized spherical cube

Figure 32: Quadrilateralized spherical cube projection ( QSC); no limits.

O'Neill & Laubscher (1976) derived an exact expression for an equal-area transformation from a sphere to the six faces of a cube. At that time, their formulation was thought to be computationally intractable, but today, with modern computers and telescopes of higher angular resolution than COBE, their formulation has come into use. Fred Patt (1993, private communication) has provided us with the inverse of the O'Neill & Laubscher formula and their expression in Cartesian coordinates.

O'Neill & Laubscher's derivation applies only in the quadrant $-45\hbox{$^\circ$ }\leq \phi \leq 45\hbox{$^\circ$ }$ and must be reflected into the other quadrants. This has the effect of making the projection non-differentiable along the diagonals as is evident in Fig. 32. To compute the forward projection first identify the face and find $(\xi,\eta,\zeta)$ and $(\phi_{\rm c},\theta_{\rm c})$ from Table 4. Then

$$(x,y) = (\phi_{\rm c},\theta_{\rm c}) + \left\{ \begin{array... ...a\vert$ } \\ (v,u) & \mbox{otherwise} \end{array} \right. ,$$ (176)

where

$$45\hbox{$^\circ$ }S~ \sqrt{\frac{1-\zeta}{1-1/\sqrt{2+\omega^2}}} ,$$ (177)

$$\frac{u}{15\hbox{$^\circ$ }} \left[ \tan^{-1}\left( \omega \right) - \sin^{-1} \left( \frac{\omega}{\sqrt{2(1+\omega^2)}} \right) \right] ,$$ (178)

$$\omega \left\{ \begin{array}{ll} \eta / \xi & \mbox{if $\vert\xi\vert > \vert\eta\vert$ } \\ \xi / \eta & \mbox{otherwise} \end{array} \right. ,$$

$$\left\{ \begin{array}{ll} +1 & \mbox{if $\xi > \vert\eta\vert$\space or $\eta > \vert\xi\vert$ } \\ -1 & \mbox{otherwise} \end{array} \right. .$$

To compute the inverse first identify the face from the (x,y) coordinates, then determine (u,v) via

$$(u,v) = \left\{ \begin{array}{ll} (x-\phi_{\rm c},y-\theta_{... ... c},x-\phi_{\rm c}) & \mbox{otherwise} \end{array} \right. .$$ (179)

Then

$$\zeta = 1 - \left(\frac{u}{45\hbox{$^\circ$ }}\right)^2 \left(1 - \frac{1}{\sqrt{2+\omega^2}} \right) ,$$ (180)

where

$$\omega = \frac{\sin(15\hbox{$^\circ$ }v/u)}{\cos(15\hbox{$^\circ$ }v/u) - 1/\sqrt{2}} \cdot$$ (181)

If $\vert x-\phi_{\rm c}\vert>\vert y-\theta_{\rm c}\vert$ then

$$\xi \sqrt{\frac{1-\zeta^2}{1+\omega^2}} ,$$ (182)

$$\eta \xi\omega ,$$ (183)

otherwise

$$\eta \sqrt{\frac{1-\zeta^2}{1+\omega^2}} ,$$ (184)

$$\xi \eta\omega .$$ (185)

Given the face number and $(\xi,\eta,\zeta)$, the native coordinates $(\phi,\theta)$ may be computed with reference to Table 4 as for the tangential spherical cube projection.