A&A 395, 1077-1122 (2002)
DOI: 10.1051/0004-6361:20021327

Representations of celestial coordinates in FITS

M. R. Calabretta¹ - E. W. Greisen²

1 - Australia Telescope National Facility, PO Box 76, Epping, NSW 1710, Australia
2 - National Radio Astronomy Observatory, PO Box O, Socorro, NM 87801-0387, USA

Received 24 July 2002 / Accepted 9 September 2002

Abstract
In Paper I, Greisen & Calabretta (2002) describe a generalized method for assigning physical coordinates to FITS image pixels. This paper implements this method for all spherical map projections likely to be of interest in astronomy. The new methods encompass existing informal FITS spherical coordinate conventions and translations from them are described. Detailed examples of header interpretation and construction are given.

Key words: methods: data analysis - techniques: image processing - astronomical data bases: miscellaneous - astrometry

   
1 Introduction

This paper is the second in a series which establishes conventions by which world coordinates may be associated with FITS (Hanisch et al. 2001) image, random groups, and table data. Paper I (Greisen & Calabretta 2002) lays the groundwork by developing general constructs and related FITS header keywords and the rules for their usage in recording coordinate information. In Paper III, Greisen et al. (2002) apply these methods to spectral coordinates. Paper IV (Calabretta et al. 2002) extends the formalism to deal with general distortions of the coordinate grid. This paper, Paper II, addresses the specific problem of describing celestial coordinates in a two-dimensional projection of the sky. As such it generalizes the informal but widely used conventions established by Greisen (1983, 1986) for the Astronomical Image Processing System, hereinafter referred to as the AIPS convention.

Paper I describes the computation of world coordinates as a multi-step process. Pixel coordinates are linearly transformed to intermediate world coordinates that in the final step are transformed into the required world coordinates.

Figure 1: Conversion of pixel coordinates to celestial coordinates. The INTERMEDIATE WORLD COORDINATES of Paper I, Fig. 1 are here interpreted as PROJECTION PLANE COORDINATES, i.e. Cartesian coordinates in the plane of projection, and the multiple steps required to produce them have been condensed into one. This paper is concerned in particular with the steps enclosed in the dotted box. For later reference, the mathematical symbols associated with each step are shown in the box at right (see also Tables 1 and 13).

Open with DEXTER

In this paper we associate particular elements of the intermediate world coordinates with Cartesian coordinates in the plane of the spherical projection. Figure 1, by analogy with Fig. 1 of Paper I, focuses on the transformation as it applies to these projection plane coordinates. The final step is here divided into two sub-steps, a spherical projection defined in terms of a convenient coordinate system which we refer to as native spherical coordinates, followed by a spherical rotation of these native coordinates to the required celestial coordinate system.

The original FITS paper by Wells et al. (1981) introduced the CRPIX ja keyword to define the pixel coordinates $(r_1,r_2,r_3,\ldots)$ of a coordinate reference point. Paper I retains this but replaces the coordinate rotation keywords CROTA i with a linear transformation matrix. Thus, the transformation of pixel coordinates $(p_1,p_2,p_3,\ldots)$ to intermediate world coordinates $(x_1,x_2,x_3,\ldots)$ becomes

 

$$s_i \sum_{j=1}^{N} m_{ij} (p_j - r_j) ,$$ (1)

$$\hphantom{s_i} \sum_{j=1}^{N} (s_i m_{ij}) (p_j - r_j) ,$$

where N is the number of axes given by the NAXIS keyword. As suggested by the two forms of the equation, the scales s_i, and matrix elements m_ij may be represented either separately or in combination. In the first form s_i is given by CDELT ia and m_ij by PC i_ja; in the second, the product s_i m_ij is given by CD i_ja. The two forms may not coexist in one coordinate representation.

Equation (1) establishes that the reference point is the origin of intermediate world coordinates. We require that the linear transformation be constructed so that the plane of projection is defined by two axes of the x_i coordinate space. We will refer to intermediate world coordinates in this plane as projection plane coordinates, (x,y), thus with reference point at (x,y) = (0,0). Note that this does not necessarily correspond to any plane defined by the p_j axes since the linear transformation matrix may introduce rotation and/or skew.

Wells et al. (1981) established that all angles in FITS were to be measured in degrees and this has been entrenched by the AIPS convention and confirmed in the IAU-endorsed FITS standard (Hanisch et al. 2001). Paper I introduced the CUNIT ia keyword to define the units of CRVAL ia and CDELT ia. Accordingly, we require CUNIT ia = 'deg' for the celestial CRVAL ia and CDELT ia, whether given explicitly or not. Consequently, the (x,y) coordinates in the plane of projection are measured in degrees. For consistency, we use degree measure for native and celestial spherical coordinates and for all other angular measures in this paper.

For linear coordinate systems Wells et al. (1981) prescribed that world coordinates should be computed simply by adding the relative world coordinates, x_i, to the coordinate value at the reference point given by CRVAL ia. Paper I extends this by providing that particular values of CTYPE ia may be established by convention to denote non-linear transformations involving predefined functions of the x_i parameterized by the CRVAL ia keyword values and possibly other parameters.

In Sects. 2, 3 and 5 of this paper we will define the functions for the transformation from (x,y)coordinates in the plane of projection to celestial spherical coordinates for all spherical map projections likely to be of use in astronomy.

The FITS header keywords discussed within the main body of this paper apply to the primary image header and image extension headers. Image fragments within binary tables extensions defined by Cotton et al. (1995) have additional nomenclature requirements, a solution for which was proposed in Paper I. Coordinate descriptions may also be associated with the random groups data format defined by Greisen & Harten (1981). These issues will be expanded upon in Sect. 4.

Section 6 considers the translation of older AIPS convention FITS headers to the new system and provisions that may be made to support older FITS-reading programs. Section 7 discusses the concepts presented here, including worked examples of header interpretation and construction.

   
2 Basic concepts

   
2.1 Spherical projection

As indicated in Fig. 1, the first step in transforming (x,y)coordinates in the plane of projection to celestial coordinates is to convert them to native longitude and latitude, $(\phi,\theta)$. The equations for the transformation depend on the particular projection and this will be specified via the CTYPE ia keyword. Paper I defined "4-3'' form for such purposes; the rightmost three-characters are used as an algorithm code that in this paper will specify the projection. For example, the stereographic projection will be coded as STG. Some projections require additional parameters that will be specified by the FITS keywords PV i_ma for $m = 0,1,2,\ldots$, also introduced in Paper I. These parameters may be associated with the longitude and/or latitude coordinate as specified for each projection. However, definition of the three-letter codes for the projections and the equations, their inverses and the parameters which define them, form a large part of this work and will be discussed in Sect. 5. The leftmost four characters of CTYPE ia are used to identify the celestial coordinate system and will be discussed in Sect. 3.

   
2.2 Reference point of the projection

The last step in the chain of transformations shown in Fig. 1 is the spherical rotation from native coordinates, $(\phi,\theta)$, to celestial coordinates $(\alpha ,\delta )$. Since a spherical rotation is completely specified by three Euler angles it remains only to define them.

In principle, specifying the celestial coordinates of any particular native coordinate pair provides two of the Euler angles (either directly or indirectly). In the AIPS convention, the CRVAL ia keyword values for the celestial axes specify the celestial coordinates of the reference point and this in turn is associated with a particular point on the projection. For zenithal projections that point is the sphere's point of tangency to the plane of projection and this is the pole of the native coordinate system. Thus the AIPS convention links a celestial coordinate pair to a native coordinate pair via the reference point. Note that this association via the reference point is purely conventional; it has benefits which are discussed in Sect. 5 but in principle any other point could have been chosen.

Section 5 presents the projection equations for the transformation of (x,y) to $(\phi,\theta)$. The native coordinates of the reference point would therefore be those obtained for (x,y) = (0,0). However, it may happen that this point lies outside the boundary of the projection, for example as for the ZPN projection of Sect. 5.1.7. Therefore, while this work follows the AIPS approach it must of necessity generalize it.

 

 

Table 1: Summary of important variable names and other symbols used throughout the paper.

Variable(s)	Meaning	Related FITS keywords (if any)
i	Index variable for world coordinates
j	Index variable for pixel coordinates
a	Alternate version code, blank or A to Z
pj	Pixel coordinates
rj	Reference pixel coordinates	CRPIX ja
mij	Linear transformation matrix	CD i_ja or PC i_ja
si	Coordinate scales	CDELT ia
xi	Intermediate world coordinates (in general)
(x,y)	Projection plane coordinates
$(\phi,\theta)$	Native longitude and latitude
$(\alpha ,\delta )$	Celestial longitude and latitude
$(\phi _0,\theta _0)$	Native longitude and latitude of the fiducial point	PV i_1a$^\dag$, PV i_2a$^\dag$
$(\alpha _0,\delta _0)$	Celestial longitude and latitude of the fiducial point	CRVAL ia
$(\phi_{\rm p},\theta_{\rm p})$	Native longitude and latitude of the celestial pole	LONPOLE a ( = PV i_3a$^\dag$), LATPOLE a, ( = PV i_4a$^\dag$)
$(\alpha_{\rm p},\delta_{\rm p})$	Celestial longitude and latitude of the native pole $(\delta_{\rm p} = \theta_{\rm p})$
$\arg()$	Inverse tangent function that returns the correct quadrant

$^\dag$ Associated with longitude axis i.

Accordingly we specify only that a fiducial celestial coordinate pair $(\alpha _0,\delta _0)$ given by the CRVAL ia will be associated with a fiducial native coordinate pair $(\phi _0,\theta _0)$ defined explicitly for each projection. For example, zenithal projections all have $(\phi _0,\theta _0) = (0,90\hbox {$^\circ $ })$, while cylindricals have $(\phi _0,\theta _0) = (0,0)$. The AIPS convention has been honored here as far as practicable by constructing the projection equations so that $(\phi _0,\theta _0)$ transforms to the reference point, (x,y) = (0,0). Thus, apart from the one exception noted, the fiducial celestial and native coordinates are the celestial and native coordinates of the reference point and we will not normally draw a distinction.

It is important to understand why $(\phi _0,\theta _0)$ differs for different projection types. There are two main reasons; the first makes it difficult for it to be the same, while the second makes it desirable that it differs. Of the former, some projections such as Mercator's, diverge at the native pole, therefore they cannot have the reference point there because that would imply infinite values for CRPIX ja. On the other hand, the gnomonic projection diverges at the equator so it can't have the reference point there for the same reason. Possibly $(\phi _0,\theta _0)$ chosen at some mid-latitude could satisfy all projections, but that leads us to the second reason.

Different projection types are best suited to different purposes. For example, zenithal projections are best for mapping the region in the vicinity of a point, often a pole; cylindrical projections are appropriate for the neighborhood of a great circle, usually an equator; and the conics are suitable for small circles such as parallels of latitude. Thus, it would be awkward if a cylindrical used to map, say, a few degrees on either side of the galactic plane, had its reference point, and thus CRPIX ja and CRVAL ia, at the native pole, way outside the map boundary. In formulating the projection equations themselves the native coordinate system is chosen to simplify the geometry as much as possible. For the zenithals the natural formulation has (x,y) = (0,0) at the native pole, whereas for the cylindricals the equations are simplest if (x,y) = (0,0) at a point on the equator.

As discussed above, a third Euler angle must be specified and this will be given by the native longitude of the celestial pole, $\phi_{\rm p}$, specified by the new FITS keyword

LONPOLE a (floating-valued).

The default value of LONPOLE a will be 0 for $\delta_0\ge\theta_0$ or $180\hbox{$^\circ$ }$ for $\delta_0<\theta_0$. This is the condition for the celestial latitude to increase in the same direction as the native latitude at the reference point. Thus, for example, in zenithal projections the default is always $180\hbox{$^\circ$ }$ (unless $\delta_0 = 90\hbox{$^\circ$ }$) since $\theta_0 = 90\hbox{$^\circ$ }$. In cylindrical projections, where $\theta_0 = 0$, the default value for LONPOLE a is 0 for $\delta_0 \ge 0$, but it is $180\hbox{$^\circ$ }$ for $\delta_0 < 0$.

2.3 Spherical coordinate rotation

Since $(\phi _0,\theta _0)$ differs for different projections it is apparent that the relationship between $(\alpha _0,\delta _0)$ and the required Euler angles also differs.

For zenithal projections, $(\phi _0,\theta _0) = (0,90\hbox {$^\circ $ })$ so the CRVAL ia specify the celestial coordinates of the native pole, i.e. $(\alpha_0,\delta_0) = (\alpha_{\rm p},\delta_{\rm p})$. There is a simple relationship between the Euler angles for consecutive rotations about the Z-, X-, and Z-axes and $\alpha_{\rm p}$, $\delta_{\rm p}$ and $\phi_{\rm p}$; the ZXZ Euler angles are $(\alpha_{\rm p}+90\hbox{$^\circ$ }, 90\hbox{$^\circ$ }-\delta_{\rm p}, \phi_{\rm p}-90\hbox{$^\circ$ })$. Given this close correspondence it is convenient to write the Euler angle transformation formulæ directly in terms of $\alpha_{\rm p}$, $\delta_{\rm p}$ and $\phi_{\rm p}$:

$$\begin{array}{rcl} \alpha & = & \alpha_{\rm p} + \arg~(\sin\... ...eta \cos\delta_{\rm p}\cos(\phi-\phi_{\rm p}) ) , \end{array}$$ (2)

where $\arg~()$ is an inverse tangent function that returns the correct quadrant, i.e. if $(x,y) = (r\cos\beta,r\sin\beta)$ with r > 0 then $\arg~(x,y) = \beta$. Note that, if $\delta_{\rm p} = 90\hbox{$^\circ$ }$, Eqs. (2) become

$$\begin{array}{rcl} \alpha & = & \alpha_{\rm p} + \phi - \phi... ... - 180\hbox{$^\circ$ }, \\ \delta & = & \theta , \end{array}$$ (3)

which may be used to define a simple change in the origin of longitude. Likewise for $\delta_{\rm p} = -90\hbox{$^\circ$ }$

$$\begin{array}{rcl} \alpha & = & \alpha_{\rm p} - \phi + \phi_{\rm p} , \\ \delta & = & -\theta . \end{array}$$ (4)

The inverse equations are

$$\begin{array}{rcl} \phi & = & \phi_{\rm p} + \arg~(\sin\delt... ... \cos\delta_{\rm p}\cos(\alpha-\alpha_{\rm p})) . \end{array}$$ (5)

Useful relations derived from Eqs. (2) and (5) are

 

$$\cos\delta \cos(\alpha - \alpha_{\rm p}) \sin\theta\cos\delta_{\rm p} - \cos\theta \sin\delta_{\rm p} \cos(\phi-\phi_{\rm p}) ,$$ =

$$\cos\delta \sin(\alpha - \alpha_{\rm p}) -\cos\theta \sin(\phi - \phi_{\rm p}) ,$$ = (6)

 

$$\cos\theta \cos(\phi - \phi_{\rm p}) \sin\delta \cos\delta_{\rm p} -\cos\delta \sin\delta _{\rm p} \cos(\alpha - \alpha_{\rm p}) ,$$ =

$$\cos\theta \sin(\phi - \phi_{\rm p}) -\cos\delta \sin(\alpha - \alpha_{\rm p}) .$$ = (7)

A matrix method of handling the spherical coordinate rotation is described in Appendix B.

   
2.4 Non-polar $(\phi _0,\theta _0)$

Projections such as the cylindricals and conics for which $(\phi_0,\theta_0) \neq (0,90\hbox{$^\circ$ })$ are handled by providing formulae to compute $(\alpha_{\rm p},\delta_{\rm p})$ from $(\alpha _0,\delta _0)$ whence the above equations may be used.

Given that $(\alpha _0,\delta _0)$ are the celestial coordinates of the point with native coordinates $(\phi _0,\theta _0)$, Eqs. (6) and (7) may be inverted to obtain

 

$$\delta_{\rm p} \arg~(\cos\theta_0 \cos(\phi_{\rm p}-\phi_0) , \sin\theta_0) \;$$

$$\pm \cos^{-1} \left( \frac{\sin\delta_0} {\sqrt{1 - \cos^2\theta_0\sin^2(\phi_{\rm p}-\phi_0)}} \right) ,$$ (8)

  

$$\sin(\alpha_0-\alpha_{\rm p}) \sin(\phi_{\rm p}-\phi_0)\cos\theta_0 / \cos\delta_0 ,$$ (9)

$$\cos(\alpha_0-\alpha_{\rm p}) \frac{\sin\theta_0 - \sin\delta_{\rm p}\sin\delta_0} {\cos\delta_{\rm p}\cos\delta_0} ,$$ (10)

whence Eqs. (2) may be used to determine the celestial coordinates. Note that Eq. (8) contains an ambiguity in the sign of the inverse cosine and that all three indicate that some combinations of $\phi _0$, $\theta _0$, $\alpha_0$, $\delta_0$, and $\phi_{\rm p}$ are not allowed. For these projections, we must therefore adopt additional conventions:

1.

Equations (9) and (10) indicate that $\alpha_{\rm p}$ is undefined when $\delta_0 = \pm90\hbox{$^\circ$ }$. This simply represents the longitude singularity at the pole and forces us to define $\alpha_{\rm p} = \alpha_0$ in this case.

2.

If $\delta_{\rm p} = \pm90\hbox{$^\circ$ }$ then the longitude at the native pole is $\alpha_{\rm p} = \alpha_0 + \phi_{\rm p} - \phi_0 - 180\hbox{$^\circ$ }$ for $\delta_{\rm p} = 90\hbox{$^\circ$ }$ and $\alpha_{\rm p} = \alpha_0 - \phi_{\rm p} + \phi_0$ for $\delta_{\rm p} = -90\hbox{$^\circ$ }$.

3.

Some combinations of $\phi _0$, $\theta _0$, $\delta_0$, and $\phi_{\rm p}$ produce an invalid argument for the $\cos^{-1}()$ in Eq. (8). This is indicative of an inconsistency for which there is no solution for $\delta_{\rm p}$. Otherwise Eq. (8) produces two solutions for $\delta_{\rm p}$. Valid solutions are ones that lie in the range $-90\hbox{$^\circ$ }$ to $+90\hbox{$^\circ$ }$, and it is possible in some cases that neither solution is valid.

Note, however, that if $\phi_0 = 0$, as is usual, then when LONPOLE a ( $\equiv \phi_{\rm p}$) takes its default value of 0 or $180\hbox{$^\circ$ }$ (depending on $\theta _0$) then any combination of $\delta_0$ and $\theta _0$ produces a valid argument to the $\cos^{-1}()$ in Eq. (8), and at least one of the solutions is valid.

4.

Where Eq. (8) has two valid solutions the one closest to the value of the new FITS keyword

LATPOLE a (floating-valued)

is chosen. It is acceptable to set LATPOLE a to a number greater than $+90\hbox{$^\circ$ }$ to choose the northerly solution (the default if LATPOLE a is omitted), or a number less than $-90\hbox{$^\circ$ }$ to select the southern solution.

5.

Equation (8) often only has one valid solution (because the other lies outside the range $-90\hbox{$^\circ$ }$ to $+90\hbox{$^\circ$ }$). In this case LATPOLE a is ignored.

6.

For the special case where $\theta_0 = 0$, $\delta_0 = 0$, and $\phi_{\rm p} - \phi_0 =\pm90\hbox{$^\circ$ }$ then $\delta_{\rm p}$ is not determined and LATPOLE a specifies it completely. LATPOLE a has no default value in this case.

These rules governing the application of Eqs. (8-10) are certainly the most complex of this formalism. FITS writers are well advised to check the values of $\phi _0$, $\theta _0$, $\alpha_0$, $\delta_0$, and $\phi_{\rm p}$ against them to ensure their validity.

   
User-specified $(\phi _0,\theta _0)$

In Sect. 2.2 we formally decoupled $(\alpha _0,\delta _0)$ from the reference point and associated it with $(\phi _0,\theta _0)$. One implication of this is that it should be possible to allow $(\phi _0,\theta _0)$to be user-specifiable. This may be useful in some circumstances, mainly to allow CRVAL ia to match a point of interest rather than some predefined point which may lie well outside the image and be of no particular interest. We therefore reserve keywords PV i_1a and PV i_2a attached to longitude coordinate i to specify $\phi _0$ and $\theta _0$ respectively.

By itself, this prescription discards the AIPS convention and lacks utility because it breaks the connection between CRVAL ia and any point whose pixel coordinates are given in the FITS header. New keywords could be invented to define these pixel coordinates but this would introduce additional complexity and still not satisfy the AIPS convention. The solution adopted here is to provide an option to force (x,y) = (0,0) at $(\phi _0,\theta _0)$ by introducing an implied offset (x₀,y₀) which is computed for $(\phi _0,\theta _0)$ from the relevant projection equations given in Sect. 5. This is to be applied to the (x,y) coordinates when converting to or from pixel coordinates. The operation is controlled by the value of PV i_0a attached to longitude coordinate i; the offset is to be applied only when this differs from its default value of zero.

This construct should be considered advanced usage, of which Figs. 33 and 34 provide an example. Normally we expect that PV i_1a and PV i_2a will either be omitted or set to the projection-specific defaults given in Sect. 5.

   
2.6 Encapsulation

So that all required transformation parameters can be contained completely within the recognized world coordinate system (WCS) header cards, the values of LONPOLE a and LATPOLE a may be recorded as PV i_3a and PV i_4a, attached to longitude coordinate i, and these take precedence where a conflict arises.

We recommend that FITS writers include the PV i_1a, PV i_2a, PV i_3a, and PV i_4a cards in the header, even if only to denote the correct use of the default values.

Note carefully that these are associated with the longitude coordinate, whereas the projection parameters defined later are all associated with the latitude coordinate.

Figure 2: A linear equatorial coordinate system (top) defined via the methods of Wells et al. (1981), and the corresponding oblique system constructed using the methods of this paper. The reference coordinate $(\alpha _0,\delta _0)$ for each is at right ascension $8^{\rm hr}$, declination $+60\hbox {$^\circ $ }$ (marked). The two sets of FITS header cards differ only in their CTYPE ia keyword values. The non-oblique graticule could be obtained in the current system by setting $(\alpha _0,\delta _0,\phi _{\rm p})=(120\hbox {$^\circ $ },0,0)$.

Open with DEXTER

2.7 Change of coordinate system

A change of coordinate system may be effected in a straightforward way if the transformation from the original system, $(\alpha ,\delta )$, to the new system, $(\alpha',\delta')$, and its inverse are known. The new coordinates of the $(\phi _0,\theta _0)$, namely $(\alpha'_0,\delta'_0)$, are obtained simply by transforming $(\alpha _0,\delta _0)$. To obtain $\phi'_{\rm p}$, first transform the coordinates of the pole of the new system to the original celestial system and then transform the result to native coordinates via Eq. (5) to obtain $(\phi'_{\rm p},\theta'_{\rm p})$. As a check, compute $\delta'_{\rm p}$ via Eq. (8) and verify that $\theta'_{\rm p} = \delta'_{\rm p}$.

   
2.8 Comparison with linear coordinate systems

It must be stressed that the coordinate transformation described here differs from the linear transformation defined by Wells et al. (1981) even for some simple projections where at first glance they may appear to be the same. Consider the plate carrée projection defined in Sect. 5.2.3 with $\phi = x$, $\theta = y$ and illustrated in Fig. 18. Figure 1 shows that while the transformation from (x,y) to $(\phi,\theta)$ may be linear (in fact identical), there still remains the non-linear transformation from $(\phi,\theta)$ to $(\alpha ,\delta )$. Hence the linear coordinate description defined by the unqualified CTYPE ia pair of RA, DEC which uses the Wells et al. prescription will generally differ from that of RA--CAR, DEC--CAR with the same CRVAL ia, etc. If LONPOLE a assumes its default value then they will agree to first order at points near the reference point but gradually diverge at points away from it.

Figure 2 illustrates this point for a plate carrée projection with reference coordinates of $8^{\rm hr}$ right ascension and $+60\hbox {$^\circ $ }$ declination and with $\phi_{\rm p} = 0$. It is evident that since the plate carrée has $(\phi _0,\theta _0) = (0,0)$, a non-oblique graticule may only be obtained by setting $\delta_0 = 0$ with $\phi_{\rm p} = 0$. It should also be noted that where a larger map is to be composed of tiled submaps the coordinate description of a submap should only differ in the value of its reference pixel coordinate.

   
3 Celestial coordinate systems

As mentioned in Sect. 2.1, Paper I defined "4-3'' form for the CTYPE ia keyword value; i.e., the first four characters specify the coordinate type, the fifth character is a "-'', and the remaining three characters specify an algorithm code for computing the world coordinate value. Thus while the right half specifies the transformation to be applied in computing the spherical coordinate pair, $(\alpha ,\delta )$, the left half simply identifies the celestial system with which $(\alpha ,\delta )$ are to be associated. In this sense CTYPE ia contains an active part which drives the transformation and a passive part which describes the results.

Consistent with past practice, an equatorial coordinate pair is denoted by RA- and DEC-, and other celestial systems are of the formxLON and xLAT for longitude and latitude pairs, where x = G for galactic, E for ecliptic, H for helioecliptic, and S for supergalactic coordinates. Since representation of planetary, lunar and solar coordinate systems could exceed the 26 possibilities afforded by the single character x, we also here allow yzLN and yzLT pairs. Additional values of x and yz will undoubtedly be defined. A basic requirement, however, is that the coordinate system be right-handed and that the pole be at latitude $+90\hbox{$^\circ$ }$. A coordinate pair such as azimuth and zenith distance would have to be represented as a negative azimuth and elevation with implied conversion. In some situations negation of the longitude coordinate may be implemented via a sign reversal of the appropriate CDELT ia. It remains the responsibility of the authors of new coordinate system types to define them properly and to gain general recognition for them from the FITS community. However, FITS interpreters should be able to associate coordinate pairs even if the particular coordinate system is not recognized.

Paper I clarified that, while the subscript of CRPIX ja and the jsubscript of PC i_ja and CD i_ja refer to pixel coordinate elements, the i subscript of PC i_ja and CD i_ja and the subscripts of CDELT ia, CTYPE ia and CRVAL ia refer to world coordinate elements. However we now have three different sets of world coordinates, (x,y), $(\phi,\theta)$, and $(\alpha ,\delta )$. This leads us to associate x, $\phi$, and $\alpha$ on the one hand, and y, $\theta$, and $\delta$ on the other. This simply means, for example, that if CTYPE3 = 'GLON-AIT', then the third element of the intermediate world coordinate calculated via Eq. (1) corresponds to what we have been calling the x-coordinate in the plane of projection, the association being between $\alpha$ and x. In this way pairs of CTYPE ia with complementary left halves and matching right halves define which elements of the intermediate world coordinate vector form the plane of projection.

   
3.1 Equatorial and ecliptic coordinates

Several systems of equatorial coordinates (right ascension and declination) are in common use. Apart from the International Celestial Reference System (ICRS, IAU 1998), the axes of which are by definition fixed with respect to the celestial sphere, each system is parameterized by time. In particular, mean equatorial coordinates (Hohenkerk et al. 1992) are defined in terms of the epoch (i.e. instant of time) of the mean equator and equinox (i.e. pole and origin of right ascension). The same applies for ecliptic coordinate systems. Several additional keywords are therefore required to specify these systems fully. We introduce the new keyword

RADESYS a (character-valued)

to specify the particular system. Recognized values are given in Table 2. Apart from FK4-NO-E these keywords are applicable to ecliptic as well as equatorial coordinates.

 

 

Table 2: Allowed values of RADESYS a.

RADESYS a	Definition
ICRS	International Celestial Reference System
FK5	mean place,
	new (IAU 1984) system
FK4	mean place,
	old (Bessell-Newcomb) system
FK4-NO-E	mean place,
	old system but without e-terms
GAPPT	geocentric apparent place,
	IAU 1984 system

Wells et al. (1981) introduced the keyword EPOCH to mean the epoch of the mean equator and equinox. However we here replace it with

EQUINOX a (floating-valued),

since the word "epoch'' is often used in astrometry to refer to the time of observation. The new keyword takes preference over EPOCH if both are given. Note that EQUINOX a applies to ecliptic as well as to equatorial coordinates.

For RADESYS a values of FK4 and FK4-NO-E, any stated equinox is Besselian and, if neither EQUINOX a nor EPOCH are given, a default of 1950.0 is to be taken. For FK5, any stated equinox is Julian and, if neither keyword is given, it defaults to 2000.0.

If the EQUINOX a keyword is given it should always be accompanied by RADESYS a. However, if it should happen to appear by itself then RADESYS a defaults to FK4 if EQUINOX a < 1984.0, or to FK5 if EQUINOX a $\ge\ 1984.0$. Note that these defaults, while probably true of older files using the EPOCH keyword, are not required of them.

RADESYSa defaults to ICRS if both it and EQUINOX a are absent.

Geocentric apparent equatorial and ecliptic coordinates (GAPPT) require the epoch of the equator and equinox of date. This will be taken as the time of observation rather than EQUINOX a. The time of observation may also be required for other astrometric purposes in addition to the usual astrophysical uses, for example, to specify when the mean place was correct in accounting for proper motion, including "fictitious'' proper motions in the conversion between the FK4 and FK5 systems. The old DATE-OBS keyword may be used for this purpose. However, to provide a more convenient specification we here introduce the new keyword

MJD-OBS (floating-valued),

that provides DATE-OBS as a Modified Julian Date ( ${\rm JD} - 2~400~000.5$) but is identical to it in all other respects. MJD-OBS does not have a version code since there can only be one time of observation. Following the year-2000 conventions for DATE keywords (Bunclark & Rots 1996), this time refers by default to the start of the observation unless another interpretation is clearly explained in the comment field. In the present case the distinction is not important. We leave it to future agreements to clarify systems of time measurement and other matters related to time.

The combination of CTYPE ia, RADESYS a, and EQUINOX a define the coordinate system of the CRVAL ia and of the celestial coordinates that result from the sequence of transformations summarized by Fig. 1. However, FK4 coordinates are not strictly spherical since they include a contribution from the elliptic terms of aberration, the so-called "e-terms'' which amount to $\leq$343 milliarcsec. Strictly speaking, therefore, a map obtained from, say, a radio synthesis telescope, should be regarded as FK4-NO-E unless it has been appropriately resampled or a distortion correction provided (Paper IV). In common usage, however, CRVAL ia for such maps is usually given in FK4 coordinates. In doing so, the e-terms are effectively corrected to first order only.

   
4 Alternate FITS image representations

4.1 Random groups visibility data

The random-groups extension to FITS (Greisen & Harten 1981) has been used to transmit interferometer fringe visibility sample data. It has been customary among users of this format to convey the suggested projection type as the last four characters of the random parameter types (PTYPE n) of the Fourier plane coordinates (u,v,w). Any rotation of these coordinates was carried, however, by CROTA i associated with the one-pixel array axis used to convey the field declination. Within the new convention users of random groups for visibility data should be prepared to read, carry forward, and use as needed, the new keywords PC i_ja, CD i_ja, LONPOLE a, LATPOLE a, PV i_ja, MJD-OBS, EQUINOX a, and RADESYS a. In particular, any rotation of the (u,v) coordinate system should be represented via PC i_ja or CD i_ja (but not both) attached to the two degenerate axes used to convey the celestial coordinates of the field center.

4.2 Pixel list and image array column elements

Paper I defined the coordinate keywords used to describe a multi-dimensional image array in a single element of a FITS binary table and a tabulated list of pixel coordinates in a FITS ASCII or binary table. In this section we extend this to the keywords specific to celestial coordinates.

4.2.1 Keyword naming convention

Table 3 lists the corresponding set of coordinate system keywords for use with each type of FITS image representation.

The data type of the alternate keyword matches that of the primary keyword and the allowed values are the same. The following notes apply to the naming conventions used in Table 3:

• a is a 1-character coordinate version code and may be blank (for the primary version of the coordinate description) or any single uppercase character from A through Z.
• n is an integer table column number without any leading zeros (1 - 999).

When using the BINTABLE image array format (Cotton et al. 1995), if the table only contains a single image column or if there are multiple image columns but they all have the same value for any of the keywords in Table 3 then the simpler form of the keyword name, as used for primary arrays, may be used. For example, if all the images in the table have the same epoch then one may use a single EQUINOX a keyword rather than multiple EQUI na keywords. The other keywords, however, must always be specified using the more complex keyword name with the column number suffix and the axis number prefix.

 

 

Table 3: Alternate coordinate keywords; the data type of the alternate keyword matches that of the primary keyword.

Keyword	Primary	BINTABLE	Pixel
Description	Array	Array	List
Coord Rotation	LONPOLE a	LONP na	LONP na
Coord Rotation	LATPOLE a	LATP na	LATP na
Coord Epoch	EQUINOX a	EQUI na	EQUI na
Date of Obs	MJD-OBS	MJDOB n	MJDOB n
Reference Frame	RADESYS a	RADE na	RADE na

In principle, more than one pixel list image can be stored in a single FITS table by defining more than one pair of p₁ and p₂ pixel coordinate columns. Under the convention defined here, however, all the images must share the same values for the keywords in Table 3.

Example binary table and pixel list headers are given in Sect. 7.3.2.

   
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)$.
Open with DEXTER

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.
Open with DEXTER

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.

Open with DEXTER

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 $ })$.
Open with DEXTER

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 $ })$.
Open with DEXTER

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$.
Open with DEXTER

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 $ }$.
Open with DEXTER

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})$.

Open with DEXTER

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.
Open with DEXTER

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.
Open with DEXTER

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.
Open with DEXTER

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 $ }$.
Open with DEXTER

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.
Open with DEXTER

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.
Open with DEXTER

   
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.
Open with DEXTER

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.
Open with DEXTER

   
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 $ }$.
Open with DEXTER

   
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.
Open with DEXTER

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

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

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

   
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 $ }$.

Open with DEXTER

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.
Open with DEXTER

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.
Open with DEXTER

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 $ }$.
Open with DEXTER

   
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.
Open with DEXTER

   
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.
Open with DEXTER

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.
Open with DEXTER

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.
Open with DEXTER

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.
Open with DEXTER

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.

   
6 Support for the AIPS convention

A large number of FITS images have been written using the AIPS coordinate convention and a substantial body of software exists to interpret it. Consequently, the AIPS convention has acquired the status of a de facto standard and FITS interpreters will need to support it indefinitely in order to obey the maxim "once FITS always FITS''. Translations between the old and new system are therefore required.

6.1 Interpreting old headers

In the AIPS convention, CROTA i assigned to the latitude axis was used to define a bulk rotation of the image plane. Since this rotation was applied after CDELT i the translation to the current formalism follows from

 

$${ \left( \begin{array}{cc} \hbox{{\tt CDELT1}} & 0 \\ 0 & \hbox{... ...t PC1\_2}} \\ \hbox{{\tt PC2\_1}} & \hbox{{\tt PC2\_2}} \end{array}\right) = }$$

$$\hspace{70pt} \left( \begin{array}{rr} \cos\rho & -\sin\rho \\ \... ...y}{cc} \hbox{{\tt CDELT1}} & 0 \\ 0 & \hbox{{\tt CDELT2}} \end{array}\right) ,$$

where we have used subscript 1 for the longitude axis, 2 for latitude, and written $\rho$ for the value of CROTA2. Equation (186), which includes the added constraint of preserving CDELT i in the translation, is readily solved for the elements of the PC i_ja matrix

$$\left( \begin{array}{cc} \hbox{{\tt PC1\_1}} & \hbox{{\tt P... ... \frac{1}{\lambda} \sin\rho & \cos\rho \end{array} \right) ,$$ (187)

where

$$\lambda = \hbox{{\tt CDELT2}} / \hbox{{\tt CDELT1}} .$$ (188)

Note that Eq. (187) defines a rotation if and only if $\lambda = \pm 1$, which is often the case. In fact, the operations of scaling and rotation are commutative if and only if the scaling is isotropic, i.e. $\lambda = +1$; for $\lambda = -1$ the direction of the rotation is reversed. However, whatever the value of $\lambda$, the interpretation of Eq. (186) as that of a scale followed by a rotation is preserved.

The translation for CD i_ja is simpler, effectively because the CDELT i have an implied value of unity and the constraint on preserving them in the translation is dropped:

 

$$\left( \begin{array}{cc} \hbox{{\tt CD1\_1}} & \hbox{{\tt CD1\_2}... ...box{{\tt CDELT1}}~\sin\rho & \hbox{{\tt CDELT2}}~\cos\rho \end{array} \right) .$$ (189)

The expressions in Hanisch & Wells (1988) and Geldzahler & Schlesinger (1997) yield the same results as Eq. (189) for the usual left-handed sky coordinates and right-handed pixel coordinates, but can lead to an incorrect interpretation (namely, possible sign errors for the off-diagonal elements) for other configurations of the coordinate systems. The Hanisch & Wells draft and the coordinate portions of Geldzahler & Schlesinger are superseded by this paper.

6.1.1 SIN

The SIN projection defined by Greisen (1983) is here generalized by the addition of projection parameters. However, these parameters assume default values which reduce to the simple orthographic projection of the AIPS convention. Therefore no translation is required.

   
6.1.2 NCP

The "north-celestial-pole'' projection defined by Greisen (1983) is a special case of the new generalized SIN projection. The old header cards

$$\begin{eqnarray*}\hbox{{\tt CTYPE1}} & = & \hbox{{\tt 'RA---NCP'}} , \\ \hbox{{\tt CTYPE2}} & = & \hbox{{\tt 'DEC--NCP'}} , \end{eqnarray*}$$

should be translated to the current formalism as

$$\begin{eqnarray*}\hbox{{\tt CTYPE1}} & = & \hbox{{\tt 'RA---SIN'}} , \\ \hbox{... ...PV2\_1}} & = & 0 , \\ \hbox{{\tt PV2\_2}} & = & \cot\delta_0 . \end{eqnarray*}$$

6.1.3 TAN, ARC and STG

The TAN, ARC and STG projections defined by Greisen (1983, 1986) are directly equivalent to those defined here and no translation is required.

   
6.1.4 AIT, GLS and MER

Special care is required in interpreting the AIT (Hammer-Aitoff), GLS (Sanson-Flamsteed), and MER (Mercator) projections in the AIPS convention as defined by Greisen (1986). As explained in Sect. 7.1, the AIPS convention cannot represent oblique celestial coordinate graticules such as the one shown in Fig. 2. CRVAL i for these projections in AIPS does not correspond to the celestial coordinates $(\alpha _0,\delta _0)$ of the reference point, as understood in this formalism, unless they are both zero in which case no translation is required.

A translation into the new formalism exists for non-zero CRVAL i but only if CROTA i is zero. It consists of setting CRVAL i to zero and adjusting CRPIX j and CDELT i accordingly in the AIPS header whereupon the above situation is obtained. The corrections to CRPIX j are obtained by computing the pixel coordinates of $(\alpha,\delta) = (0,0)$ within the AIPS convention. For AIT and MER (but not GLS), CDELT i must also be corrected for the scaling factors $f_\alpha$ and $f_\delta$incorporated into the AIPS projection equations.

Of the three projections only GLS is known to have been used with non-zero CRVAL i. Consequently we have renamed it as SFL as a warning that translation is required.

6.2 Supporting old interpreters

As mentioned in Sect. 6, FITS interpreters will need to recognize the AIPS convention virtually forever. It stands to reason, therefore, that if modern FITS-writers wish to assist older FITS interpreters they may continue to write older style headers, assuming of course that it is possible to express the coordinate system in the AIPS convention.

Modern FITS-writers must not attempt to help older interpreters by including CROTA i together with the new keyword values (assuming the combination of CDELT i and PC i_ja matrix, or CD i_ja matrix, is amenable to such translation). We make this requirement primarily to minimize confusion.

Assuming that a header has been developed using the present formalism the following test may be applied to determine whether the combination of CDELT ia and PC i_ja matrix represents a scale followed by a rotation as in Eq. (186). Firstly write

 

$$\left( \begin{array}{cc} \hbox{{\tt CD1\_1}} & \hbox{{\tt CD1\_2}... ...\tt PC1\_2}} \\ \hbox{{\tt PC2\_1}} & \hbox{{\tt PC2\_2}} \end{array}\right) ,$$ (190)

where 1 is the longitude coordinate and 2 the latitude coordinate, then evaluate $\rho_{\rm a}$ and $\rho_{\rm b}$ as

$$\rho_{\rm a} \left\{ \begin{array}{ll} \arg~(\hphantom{-}\hbox{{\tt CD1\_1}}, ... ...D2\_1}}) & \quad\mbox{if \quad $\hbox{{\tt CD2\_1}} < 0$ } \end{array}\right. ,$$

$$\rho_{\rm b} \left\{ \begin{array}{ll} \arg~( - \hbox{{\tt CD2\_2}}, \hphantom... ...1\_2}}) & \quad \mbox{if \quad $\hbox{{\tt CD1\_2}} < 0$ } \end{array}\right. .$$ (191)

If $\rho_{\rm a} = \rho_{\rm b}$ to within reasonable precision (Geldzahler & Schlesinger 1997), then compute

$$\rho = (\rho_{\rm a} + \rho_{\rm b}) / 2$$ (192)

as the best estimate of the rotation angle, the older keywords are

 

$$\hbox{{\tt CDELT1}} \hbox{{\tt CD1\_1}} / \cos\rho ,$$ =

$$\hbox{{\tt CDELT2}} \hbox{{\tt CD2\_2}} / \cos\rho ,$$ = (193)

$$\hbox{{\tt CROTA2}} \rho .$$ =

Note that the translated values of CDELT i in Eqs. (193) may differ from the starting values in Eq. (190).

Solutions for CROTA2 come in pairs separated by $180\hbox{$^\circ$ }$. The above formulæ give the solution which falls in the half-open interval $[0,180\hbox{$^\circ$ })$. The other solution is obtained by subtracting $180\hbox{$^\circ$ }$ from CROTA2 and negating CDELT1 and CDELT2. While each solution is equally valid, if one makes $\hbox{{\tt CDELT1}} < 0$ and $\hbox{{\tt CDELT2}} > 0$ then it would normally be the one chosen.

Of course, the projection must be one of those supported by the AIPS convention, which only recognizes SIN, NCP, TAN, ARC, STG, AIT, GLS and MER. Of these, we strongly recommend that the AIPS version of AIT, GLS, and MER not be written because of the problems described in Sect. 6.1.4. It is interesting to note that a translation does exist for the slant orthographic (SIN) projection defined in Sect. 5.1.5 to the simple orthographic projection of AIPS. However, we advise against such translation because of the likelihood of creating confusion and so we do not define it here. The exception is where the SIN projection may be translated as NCP as defined in Sect. 6.1.2.

   
7 Discussion

   
7.1 Oblique projections

The term oblique projection is often used when a projection's native coordinate system does not coincide with the coordinate system of interest. Texts on terrestrial cartography often separately name and derive projection equations for particular oblique projections. Thus the Cassini, Transverse Mercator, and Bartholomew nordic projections are nothing more than the plate carrée, Mercator, and Mollweide projections in disguise.

The view of obliqueness as being a property of a projection arose mainly because of the computational difficulty of producing oblique graticules in the days before electronic computers. The particular aspects chosen were those for which geometrical construction was possible or for which the mathematical formulation had a simple form, and this tended to entrench particular oblique projections as separate entities. In addition, terrestrial longitude and latitude are so closely tied to the visual representation of the Earth's surface that it is shown predominantly in the usual north-south orientation. The only other natural terrestrial coordinate system is that defined by the magnetic pole but the difference between it and the geographic graticule is sufficiently small that it is usually treated as a local correction to the magnetic bearing. The special view of obliqueness was probably also reinforced by traditional methods of map making using sextant and chronometer, which were based on geographic longitude and latitude.

The situation in astronomical cartography is quite different. The celestial sphere has a variety of natural coordinate systems - equatorial, ecliptic, galactic, etc. - and oblique and non-oblique graticules are often plotted together on the same map. It wouldn't make sense to describe such a projection as being simultaneously oblique and non-oblique; clearly it is the particular coordinate graticule which may be oblique, not the projection. Visually, an area of the sky may be seen in different orientations depending on whether it's rising, transiting, or setting, or whether seen from the northern or southern hemisphere. Moreover, oblique graticules arise as a normal feature of observations in optical, infrared, radio, and other branches of astronomy, just as they do now in terrestrial mapping based on aerial and satellite photography. The center of the field of view, wherever it may happen to be, typically corresponds to the native pole of a zenithal projection. Thus the aim of this paper has been to provide a formalism whereby obliquity may be handled in a general way.

Figures 33 and 34 illustrate the same four oblique celestial graticules for the zenithal equal area, conic equidistant, Hammer-Aitoff, and COBE quadrilateralized spherical cube projections (at variable (x,y) scale). For the sake of intercomparison, these graticules are defined in terms of the celestial coordinates of the native pole, $(\alpha_{\rm p},\delta_{\rm p})$, together with $\phi_{\rm p}$ by setting $(\phi _0,\theta _0) = (0,90\hbox {$^\circ $ })$ as described in Sect. 2.5.

The first graticule, A, when compared to the non-oblique native coordinate graticules presented earlier for each projection, illustrates the effect of changing $\delta_{\rm p}$ (and hence $\delta_0$). Comparison of graticules A and B shows that changing $\alpha_{\rm p}$ (and hence $\alpha_0$) results in a simple change in origin of longitude. Graticules A, C, and D show the more interesting effect of varying $\phi_{\rm p}$ (conveyed by LONPOLE a). For the zenithal projections the result is indistinguishable from a bulk rotation of the image plane, but this is not the case for any other class of projection. This explains why the role of LONPOLE a was covered by CROTA i in the AIPS convention for the zenithal projections introduced by Greisen (1983), but why this does not work for any other class of projection.

Figure 33: Oblique $(\alpha ,\delta )$ graticules plotted for the zenithal equal area ( ZEA) and conic equidistant ( COD) projections with parameters (A) $\alpha _{\rm p} = 0$, $\delta _{\rm p} = 30\hbox {$^\circ $ }$, $\phi _{\rm p} = 180\hbox {$^\circ $ }$; (B) $\alpha _{\rm p} = 45\hbox {$^\circ $ }$, $\delta _{\rm p} = 30\hbox {$^\circ $ }$, $\phi _{\rm p} = 180\hbox {$^\circ $ }$; (C) $\alpha _{\rm p} = 0$, $\delta _{\rm p} = 30\hbox {$^\circ $ }$, $\phi _{\rm p} = 150\hbox {$^\circ $ }$; (D) $\alpha _{\rm p} = 0$, $\delta _{\rm p} = 30\hbox {$^\circ $ }$, $\phi _{\rm p} = 75\hbox {$^\circ $ }$.

Open with DEXTER

	Figure 34: Oblique ( $\alpha ,\delta$) graticules plotted for the Hammer-Aitoff ( AIT) and COBE quadrilateralized spherical cube ( CSC) projections using the same parameters as Fig. 33.
Open with DEXTER

   
7.2 Choice of projection

The projected meridians and parallels in Figs. 33 and 34 also serve to illustrate the distortions introduced by the various projections. In particular, the quad-cube projection, while doing a good job within each face, is very distorted over the sphere as a whole, especially where the faces join, so much so that it is difficult even to trace the path of some of the meridians and parallels. However, as indicated previously, this projection is designed for efficient numerical computation rather than visual representation of the sphere.

This leads us to the question of the choice of projection for a particular application. In some cases the projection is the natural result of the geometry of the observation and there is no choice. For example, maps produced by rotational synthesis radio telescopes are orthographic (SIN) projections with native pole at the phase center of the observation. Photographic plates produced by Schmidt telescopes are best described by the zenithal equidistant (ARC) projection, while the field of view of other optical telescopes is closer to a gnomonic (TAN) projection. On the other hand, the great circle scanning technique of the Sloan Digital Sky Survey produces a cylindrical projection. Similarly a map of the surface of the Moon as it appears from Earth requires a zenithal perspective (AZP) projection with $\mu \approx -220$. The same is true for spacecraft generated images of distant moons and planets. All-sky cameras should be well served by a ZPN projection with empirically determined parameters.

Sometimes observational data will be regridded onto a projection chosen for a particular purpose. Equal area projections are often used since they preserve surface density and allow integration, whether numerical or visual, to be performed with reasonable accuracy by summing pixel values. The Hammer-Aitoff (AIT) projection is probably the most widely used for all-sky maps. However, the Sanson-Flamsteed (SFL) may be preferred as being easier to take measurements from. The humble plate carrée (CAR) excels in this regard and may be considered adequate, say, for mapping a few degrees on either side of the galactic plane. In general terms, zenithal projections are good for mapping the region in the vicinity of a point, often a pole; cylindrical projections are good for the neighborhood of a great circle, usually an equator; and the conics are suitable for small circles such as parallels of latitude. A conic projection might be a good choice for mapping a hemisphere; distortion at most points is reduced compared to a zenithal projection although at the expense of the break between native longitudes $\pm180\hbox{$^\circ$ }$. However, in some cases this break might be considered an extreme distortion which outweighs other reductions in distortion and, depending on the application, may mandate the use of a zenithal projection. Cartographers typically favor conics and polyconics for "Australia-sized'' regions of the sphere and at this they excel. Oblique forms of the zenithal equal-area projection are also commonly used. Celestial applications might include large scale maps of the Magellanic clouds. Tables 3.1 and 4.1 of Snyder (1993) summarize actual usage in a variety of 19th and 20th century atlases.

As mentioned in Sect. 5, some projections are not scaled true at the reference point. In most cases this is deliberate, usually because the projection is designed to minimize distortion over a wide area; the scale will be true at some other place on the projection, typically along a parallel of latitude. If the projection is required to have $(\alpha,\delta)\approx(\alpha_0,\delta_0)+(x,y)$ then a number of projections will serve provided also that $\phi_{\rm p}$ takes its default value. The following are always scaled true at the reference point: SZP, TAN, SIN, STG, ARC, ZEA, CAR, MER, BON, PCO, SFL, and AIT. The following are scaled true at the reference point for the particular conditions indicated: AZP ( $\gamma = 0$), ZPN ( P₀=0, P₁=1), AIR ( $\theta_b=90\hbox{$^\circ$ }$), CYP ($\lambda = 1$), CEA ($\lambda = 1$), COP ($\eta=0$), COD ($\eta=0$), COE ($\eta=0$), and COO ($\eta=0$). The following are never scaled true at the reference point: PAR (but is scaled true in x), MOL, TSC, CSC, and QSC; the latter three are equiscaled at the reference point.

Of course CDELT ia may be used to control scaling of (x,y) with respect to (p₁,p₂). Thus, for example, the plate carrée projection (CAR) also serves as the more general equirectangular projection.

   
7.3 Header interpretation examples

We now consider three examples chosen to illustrate how a FITS reader would interpret a celestial coordinate header.

Example 1 is a simple header whose interpretation is quite straightforward.

Example 2 is more complicated; it also serves to illustrate the WCS header cards for 1) image arrays in binary tables, and 2) pixel lists.

Example 3 highlights a subtle problem introduced into a header by the FITS writer and considers how this should be corrected. As such, it introduces the generally more difficult task of composing WCS headers, considered in greater detail in Sect. 7.4.

   
7.3.1 Header interpretation example 1

 

 

Table 5: FITS header for example 1.

123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789

NAXIS = 4 / 4-dimensional cube

NAXIS1 = 512 / x axis (fastest)

NAXIS2 = 512 / y axis (2nd fastest)

NAXIS3 = 196 / z axis (planes)

NAXIS4 = 1 / Dummy to give a coordinate

CRPIX1 = 256 / Pixel coordinate of reference point

CDELT1 = -0.003 / 10.8 arcsec per pixel

CTYPE1 = 'RA---TAN' / Gnomonic projection

CRVAL1 = 45.83 / RA at reference point

CUNIT1 = 'deg ' / Angles are degrees always

CRPIX2 = 257 / Pixel coordinate of reference point

CDELT2 = 0.003 / 10.8 arcsec per pixel

CTYPE2 = 'DEC--TAN' / Gnomonic projection

CRVAL2 = 63.57 / Dec at reference point

CUNIT2 = 'deg ' / Angles are degrees always

CRPIX3 = 1 / Pixel coordinate of reference point

CDELT3 = 7128.3 / Velocity increment

CTYPE3 = 'VELOCITY' / Each plane at a velocity

CRVAL3 = 500000.0 / Velocity in m/s

CUNIT3 = 'm/s ' / metres per second

CRPIX4 = 1 / Pixel coordinate of reference point

CDELT4 = 1 / Required here.

CTYPE4 = 'STOKES ' / Polarization

CRVAL4 = 1 / Unpolarized

CUNIT4 = ' ' / Conventional unitless = I pol

LONPOLE = 180 / Native longitude of celestial pole

RADESYS = 'FK5 ' / Mean IAU 1984 equatorial coordinates

EQUINOX = 2000.0 / Equator and equinox of J2000.0

Consider as the first example the relatively simple optical image whose header is given in Table 5. The NAXIS and NAXIS j keywords indicate that we have a four-dimensional image consisting of 512 columns, 512 rows, 196 planes, and one polarization. The degenerate STOKES axis is used simply to convey a coordinate value which applies to every pixel in the image. The CRPIX j keywords tell us that the reference point is at pixel coordinate (256,257,1,1). The PC i_ja keywords default to the unit matrix which indicates that no bulk rotation or shear is applied to the pixel coordinates. Intermediate world coordinates may thus be computed from Eq. (1) as

$$\begin{eqnarray*}\left( \begin{array}{c} x \\ y \\ z \\ s \end{array} \... ...- & 257 \\ p_3 & - & 1 \\ p_4 & - & 1 \end{array} \right) . \end{eqnarray*}$$

Since VELOCITY and STOKES are linear axis types, the velocity and Stokes value of each point are found simply by adding the coordinate value at the reference point to the relative coordinate. Thus,

$$\begin{eqnarray*}{\rm Velocity} & = & 500000.0 + 7128.3~(p_3 - 1)\ {\rm m~s}^{-1} , \\ {\rm Stokes} & = & 1\ {\rm (I\ polarization)} . \end{eqnarray*}$$

The CTYPE1 and CTYPE2 keywords denote a TAN (gnomonic) projection for which (x,y) are projection plane coordinates for the zenithal perspective projection with the source of the projection at the center of the sphere. Thus the native longitude and latitude are given by

$$\begin{eqnarray*}\phi & = & \arg~(-y,x) , \\ \theta & = & \tan^{-1} \left( \frac{180\hbox{$^\circ$ }}{\pi}\frac{1}{\sqrt{x^2+y^2}} \right) , \end{eqnarray*}$$

which, on substitution, become

$$\begin{eqnarray*}\phi & = & \arg~(p_2-257, p_1-256) + 180\hbox{$^\circ$ }, \\ ... ...$.\!\!^\circ$ }5932}{\sqrt{ (p_1-256)^2+(p_2-257)^2}} \right) , \end{eqnarray*}$$

where we have used Eqs. (15), (14), and (55).

The celestial coordinate system is equatorial since the CTYPE ia begin with RA and DEC and the RADESYS a and EQUINOX a cards denote that these are in the IAU 1984 system. Zenithal projections have $(\phi _0,\theta _0) = (0,90\hbox {$^\circ $ })$ for which the CRVAL i give equatorial coordinates $\alpha_{\rm p}=45\hbox{$.\!\!^\circ$ }83$ in right ascension and $\delta_{\rm p}=63\hbox{$.\!\!^\circ$ }57$ in declination. The equatorial north pole has a longitude of $180\hbox{$^\circ$ }$ in the native coordinate system from the LONPOLE a keyword. LATPOLE a is never required for zenithal projections and was not given. Thus, Eqs. (6) for the right ascension and declination become

$$\begin{eqnarray*}\begin{array}{l} \sin\delta = \sin\theta \sin(63\hbox{$.\!\!^\... ... \cos\theta\cos\phi\sin(63\hbox{$.\!\!^\circ$ }57) . \end{array}\end{eqnarray*}$$

Sample calculations for points on the diagonal near the three corners of the image are shown in Table 6.

If we define the projection non-linearity as the departure of $\theta$ from $r = \sqrt{x^2 + y^2}$, then in this $1\hbox{$.\!\!^\circ$ }5 \times 1\hbox{$.\!\!^\circ$ }5$ image it amounts to $\sim 0\hbox{$.\!\!^{\prime\prime}$ }5$ at the corners. However, in a $6\hbox{$^\circ$ }\times 6\hbox{$^\circ$ }$ image it quickly escalates to $\sim 0\hbox{$.\mkern-4mu^\prime$ }5$, sixty times larger. Comparison of $\alpha$ for the southeast and northeast corners indicates the significant effect of grid convergence even in this moderate sized image. The two differ by $22\hbox{$.\!\!^{\rm s}$ }16$, most of which is attributable to the $\cos\delta$ term in converting longitude offsets to true angular distances. The effect of grid convergence is small for $(\alpha _0,\delta _0)$ near the equator and large near the poles.

Figure 35 investigates the effect of projective non-linearities for moderate field sizes for the commonly used zenithal projections. It shows the difference in $R_{\theta}$ between the various projections and the SIN projection as a function of native latitude $\theta$. The SIN projection is used as reference since it always has $R_{\theta}$ less than the other zenithal projections. The difference for all projections exceeds 1 arcsec for values of $\theta$ less than $88\hbox{$^\circ$ }$ and the difference for the TAN projection exceeds one milliarcsec only 440 arcsec from the native pole. Such nonlinearities would be significant in optical, VLBI, and other high resolution observations.

   
7.3.2 Header interpretation example 2

While the previous header was a realistic example it overlooked many of the concepts introduced in this paper. Consider now the header given in Table 7. Although contrived to illustrate as much as possible in one example, this header could conceivably arise as a "tile'' from a conic equal area projection of a region of the southern galactic hemisphere centered at galactic coordinates $(\ell,b) = (90\hbox{$^\circ$ },-25\hbox{$^\circ$ })$. The tile size of $2048 \times 2048$ pixels is approximately $10\hbox{$^\circ$ }\times 10\hbox{$^\circ$ }$, and this particular tile is situated immediately to the north of the central tile.

 

 

Table 6: Sample calculations for example 1.

parameter	units	SE corner	NE corner	NW corner
(p1,p2)	pixel	(1, 2)	(1, 512)	(511, 512)
(p3,p4)	pixel	(1, 1)	(1, 1)	(196, 1)
x	deg	$0\hbox{$.\!\!^\circ$ }765000$	$0\hbox{$.\!\!^\circ$ }765000$	$-0\hbox{$.\!\!^\circ$ }765000$
y	deg	$-0\hbox{$.\!\!^\circ$ }765000$	$0\hbox{$.\!\!^\circ$ }765000$	$0\hbox{$.\!\!^\circ$ }765000$
$\phi$	deg	$45\hbox{$.\!\!^\circ$ }000000$	$135\hbox{$.\!\!^\circ$ }000000$	$225\hbox{$.\!\!^\circ$ }000000$
$\theta$	deg	$88\hbox{$.\!\!^\circ$ }918255$	$88\hbox{$.\!\!^\circ$ }918255$	$88\hbox{$.\!\!^\circ$ }918255$
$\alpha$	deg	$47\hbox{$.\!\!^\circ$ }503264$	$47\hbox{$.\!\!^\circ$ }595581$	$44\hbox{$.\!\!^\circ$ }064419$
$\delta$	deg	$62\hbox{$.\!\!^\circ$ }795111$	$64\hbox{$.\!\!^\circ$ }324332$	$64\hbox{$.\!\!^\circ$ }324332$
Velocity	ms-1	500000.00	500000.00	1890018.50
Stokes		$1.0 \equiv$ I	$1.0 \equiv$ I	$1.0 \equiv$ I

	Figure 35: $R_{\theta }(projection) - R_{\theta }(\hbox {{\tt SIN}})$ in arcsec plotted versus $\theta$ for various projections.
Open with DEXTER

The header defines ecliptic coordinates as an alternate coordinate system, perhaps to help define the distribution of zodiacal light. This has the same reference point and transformation matrix as the primary description and the reference values are translated according to the prescription given in Sect. 2. The RADESYSA card indicates that the ecliptic coordinates are mean coordinates in the IAU 1984 system, though the author of the header has been sloppy in omitting the EQUINOXA card which therefore defaults to J2000.0. Note that, in accord with Paper I, all keywords for the alternate description are reproduced, even those which do not differ from the primary description.

The problem will be to determine the galactic and ecliptic coordinates corresponding to a field point with pixel coordinates (1957.2,775.4). We begin as before by substituting in Eq. (1)

$$\begin{eqnarray*}\left( \begin{array}{c} x \\ y \end{array} \right) & = & ... ... - & 1024.5 \\ p_2 & + & 1023.5 \\ \end{array} \right) \cdot \end{eqnarray*}$$

Note that this transformation matrix describes a slight skewness and rotation; the x-, and y-axes correspond to the following two non-orthogonal lines in the pixel plane;

$$\begin{eqnarray*}p_2 & = & 0.002 \left( p_1 - 1024.5 \right) - 1023.5 , \\ * p_2 & = & \hphantom{0.} 250 \left( p_1 - 1024.5 \right) - 1023.5 . \end{eqnarray*}$$

The values of the (x,y) projection plane coordinates for the chosen point are shown in Table 8 together with the remaining calculations for this example.

The next step is to compute the native longitude and latitude. Like all standard conics, the conic equal area projection has two parameters, $\theta_{\rm a}$ and $\eta$, given by the keywords PV i_1a and PV i_2a attached to latitude coordinate i. In this example PV2_1 is given but the author has again been sloppy in omitting PV2_2, which thus defaults to 0. Explicit inclusion of this keyword in the header would obviate any suspicion that it had been accidentally omitted. The native coordinates, $(\phi,\theta)$ are obtained from (x,y) using Eqs. (117) and (129) via the intermediaries C, Y₀, and $R_{\theta}$ given by Eqs. (125), (127), and (116). These in turn are expressed in terms of $\gamma$, $\theta_1$, and $\theta_2$ given by Eqs. (128), (119), and (120). The calculations are shown in Table 8.

 

 

Table 7: Second example FITS header (blank lines have been inserted for clarity).

123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789

NAXIS = 2 / 2-dimensional image

NAXIS1 = 2048 / x axis (fastest)

NAXIS2 = 2048 / y axis (2nd fastest)

MJD-OBS = 44258.7845612 / MJD at start of observation

CRPIX1 = 1024.5 / Pixel coordinate of reference point

CRPIX2 = -1023.5 / Pixel coordinate of reference point

PC1_1 = 1 / Transformation matrix element

PC1_2 = -0.004 / Transformation matrix element

PC2_1 = -0.002 / Transformation matrix element

PC2_2 = 1 / Transformation matrix element

CDELT1 = -0.005 / x-scale

CDELT2 = 0.005 / y-scale

CTYPE1 = 'GLON-COE' / Conic equal area projection

CTYPE2 = 'GLAT-COE' / Conic equal area projection

PV2_1 = -25.0 / Conic mid-latitude

CRVAL1 = 90.0 / Galactic longitude at reference point

CRVAL2 = -25.0 / Galactic latitude at reference point

CRPIX1A = 1024.5 / Pixel coordinate of reference point

CRPIX2A = -1023.5 / Pixel coordinate of reference point

PC1_1A = 1 / Transformation matrix element

PC1_2A = -0.004 / Transformation matrix element

PC2_1A = -0.002 / Transformation matrix element

PC2_2A = 1 / Transformation matrix element

CDELT1A = -0.005 / x-scale

CDELT2A = 0.005 / y-scale

CTYPE1A = 'ELON-COE' / Conic equal area projection

CTYPE2A = 'ELAT-COE' / Conic equal area projection

PV2_1A = -25.0 / Conic mid-latitude

CRVAL1A = -7.0300934 / Ecliptic longitude at reference point

CRVAL2A = 34.8474143 / Ecliptic latitude at reference point

LONPOLEA= 6.3839706 / Native longitude of ecliptic pole

LATPOLEA= 29.8114400 / Ecliptic latitude of native pole

RADESYSA= 'FK5 ' / Mean IAU 1984 ecliptic coordinates

At this stage we have deprojected the field point. In this example the alternate coordinate description defines the same projection as the primary description. Indeed, it may seem odd that the formalism even admits the possibility that they may differ. However, this is a realistic possibility, for example in defining multiple optical plate solutions based on the TAN projection. It now remains to transform the native spherical coordinates into galactic coordinates, $(\ell,b)$, and ecliptic coordinates, $(\lambda,\beta)$. To do this we need to apply Eqs. (2) and in order to do that we need $(\ell_{\rm p},b_{\rm p})$ and $(\lambda_{\rm p},\beta_{\rm p})$.

Looking first at galactic coordinates, we have $(\ell_0,b_0) = (90\hbox{$^\circ$ },-25\hbox{$^\circ$ })$ and $\phi_{\rm p} = 0$. This conic projection has $(\phi_0,\theta_0) = (0\hbox{$^\circ$ },-25\hbox{$^\circ$ })$, and because $b_0 = \theta_0$ and $\phi_{\rm p} = 0\hbox{$^\circ$ }$, the native and galactic coordinate systems must coincide to within an offset in longitude. However, it is not obvious what $\ell_{\rm p}$ should be to produce this offset. Equation (8) has one valid solution, namely $b_{\rm p}=90\hbox{$^\circ$ }$. The special case in point 2 in the usage notes for Eqs. (9) and (10) must therefore be used to obtain $(\ell_{\rm p},b_{\rm p}) = (-90\hbox{$^\circ$ },90\hbox{$^\circ$ })$. The galactic coordinates of the field point listed in Table 8 are then readily obtained by application of Eqs. (2).

The header says that the ecliptic coordinates of the reference point are $(\lambda_0,\beta_0) = (-7\hbox{$.\!\!^\circ$ }0300934,34\hbox{$.\!\!^\circ$ }8474143)$ and the native longitude of the ecliptic pole is $\phi_{\rm p} = 6\hbox{$.\!\!^\circ$ }3839706$. It also specifies LATPOLEA as $29\hbox{$.\!\!^\circ$ }8114400$. In this case Eq. (8) has two valid solutions, $\beta_{\rm p} = -25\hbox{$.\!\!^\circ$ }1367794 \pm 54\hbox{$.\!\!^\circ$ }9482194$, and the one closest in value to LATPOLEA (in fact equal to it) is chosen. If LATPOLEA had been omitted from the header its default value of $+90\hbox{$^\circ$ }$ would have selected the northerly solution anyway, but of course it is good practice to make the choice clear. The value of $\lambda_{\rm p}$ may be obtained by a straightforward application of Eqs. (9) and (10), and the ecliptic coordinates of the field point computed via Eqs. (2) are listed in Table 8 as the final step of the calculation. The reader may verify the calculation by transforming the computed galactic coordinates of the field point to mean ecliptic coordinates.

7.3.3 Binary table representations of example 2

Table 9 shows the FITS header for a set of images stored in binary table image array column format. The images are similar to that described in header interpretation example 2, Sect. 7.3.2, with primary image header illustrated in Table 7.

In this example the images are stored as 2-D arrays in Col. 5 of the table and each row of the table contains a $2048 \times 2048$ pixel image of a different region on the sky. This might represent a set of smaller images extracted from a single larger image. In this case all coordinate system parameters except for the reference pixel coordinate are the same for each image and are given as header keywords. The reference pixel coordinate for the primary and secondary description are given in Cols. 1 to 4 of the binary table.

Table 10 shows the header for the same image given in pixel list format. There are 10 000 rows in this table, each one listing the pixel coordinates (XPOS, YPOS) of every detected "event'' or photon in the image. For illustration purposes, this table also contains an optional DATA_QUALITY column that could be used to flag the reliability or statistical significance of each event. A real image may be constructed from this virtual image as the 2-dimensional histogram of the number of events that occur within each pixel. The additional TLMIN n and TLMAX n keywords shown here are used to specify the minimum and maximum legal values in XPOS and YPOS columns and are useful for determining the range of each axis when constructing the image histogram.

 

 

Table 8: Calculations for example 2.

(p1,p2)	1957.2	775.4
(x,y)	$-4\hbox{$.\!\!^\circ$ }6275220$	$8\hbox{$.\!\!^\circ$ }9851730$
$\theta_{\rm a},\eta$	$-25\hbox{$.\!\!^\circ$ }0000000$	$0\hbox{$.\!\!^\circ$ }0000000$
$\theta_1,\theta_2$	$-25\hbox{$.\!\!^\circ$ }0000000$	$-25\hbox{$.\!\!^\circ$ }0000000$
$\gamma$	-0.8452365
C,Y0	-0.4226183	$-122\hbox{$.\!\!^\circ$ }8711957$
$(\phi,\theta)$	$-4\hbox{$.\!\!^\circ$ }7560186$	$-15\hbox{$.\!\!^\circ$ }8973800$
$(\ell_{\rm p},b_{\rm p})$	$-90\hbox{$.\!\!^\circ$ }0000000$	$90\hbox{$.\!\!^\circ$ }0000000$
$(\ell,b)$	$85\hbox{$.\!\!^\circ$ }2439814$	$-15\hbox{$.\!\!^\circ$ }8973800$
$(\lambda_{\rm p},\beta_{\rm p})$	$-179\hbox{$.\!\!^\circ$ }9767827$	$29\hbox{$.\!\!^\circ$ }8114400$
$(\lambda,\beta)$	$-14\hbox{$.\!\!^\circ$ }7066741$	$43\hbox{$.\!\!^\circ$ }0457292$

   
7.3.4 Header interpretation example 3

This example has been adapted from a real-life FITS data file. The simplicity of the header shown in Table 11 is deceptive; it is actually presented as an example of how not to write a FITS header, although the latent problem with its interpretation is quite subtle.

 

 

Table 9: Example binary table image array header (blank lines have been inserted for clarity).

123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789

XTENSION= 'BINTABLE' / Binary table extension

BITPIX = 8 / 8-bit bytes

NAXIS = 2 / 2-dimensional binary table

NAXIS1 = 16777232 / Width of table in bytes

NAXIS2 = 4 / Number of rows in table

PCOUNT = 0 / Size of special data area

GCOUNT = 1 / One data group (required keyword)

TFIELDS = 5 / number of fields in each row

TTYPE1 = '1CRP5 ' / Axis 1: reference pixel coordinate

TFORM1 = '1J ' / Data format of column 1: I*4 integer

TTYPE2 = '2CRP5 ' / Axis 2: reference pixel coordinate

TFORM2 = '1J ' / Data format of column 2: I*4 integer

TTYPE3 = '1CRP5A ' / Axis 1A: reference pixel coordinate

TFORM3 = '1J ' / Data format of column 3: I*4 integer

TTYPE4 = '2CRP5A ' / Axis 2A: reference pixel coordinate

TFORM4 = '1J ' / Data format of column 4: I*4 integer

TTYPE5 = 'Image ' / 2-D image array

TFORM5 = '4194304J' / Data format of column 5: I*4 vector

TDIM5 = '(2048,2048)' / Dimension of column 5 array

MJDOB5 = 44258.7845612 / MJD at start of observation

COMMENT The following keywords define the primary coordinate description

COMMENT of the images contained in Column 5 of the table.

11PC5 = 1 / Transformation matrix element

12PC5 = -0.004 / Transformation matrix element

21PC5 = -0.002 / Transformation matrix element

22PC5 = 1 / Transformation matrix element

1CDE5 = -0.005 / Axis 1: scale

2CDE5 = 0.005 / Axis 2: scale

1CTY5 = 'GLON-COE' / Axis 1: conic equal area projection

2CTY5 = 'GLAT-COE' / Axis 2: conic equal area projection

2PV5_1 = -25.0 / Conic mid-latitude

1CRV5 = 90.0 / Axis 1: galactic longitude at reference point

2CRV5 = -25.0 / Axis 2: galactic latitude at reference point

COMMENT The following keywords define the secondary coordinate description

COMMENT of the images contained in Column 5 of the table.

11PC5A = 1 / Transformation matrix element

12PC5A = -0.004 / Transformation matrix element

21PC5A = -0.002 / Transformation matrix element

22PC5A = 1 / Transformation matrix element

1CDE5A = -0.005 / Axis 1A: scale

2CDE5A = 0.005 / Axis 2A: scale

1CTY5A = 'ELON-COE' / Axis 1A: conic equal area projection

2CTY5A = 'ELAT-COE' / Axis 2A: conic equal area projection

2PV5_1A = -25.0 / Conic mid-latitude

1CRV5A = -7.0300934 / Axis 1A: ecliptic longitude at reference point

2CRV5A = 34.8474143 / Axis 2A: ecliptic latitude at reference point

LONP5A = 6.3839706 / Native longitude of ecliptic pole

LATP5A = 29.8114400 / Ecliptic latitude of native pole

RADE5A = 'FK5 ' / Mean IAU 1984 ecliptic coordinates

EQUI5A = 2000.0 / Coordinate epoch

END

 

 

Table 10: Example pixel list header (blank lines have been inserted for clarity).

123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789

XTENSION= 'BINTABLE' / Binary table extension

BITPIX = 8 / 8-bit bytes

NAXIS = 2 / 2-dimensional binary table

NAXIS1 = 5 / Width of table in bytes

NAXIS2 = 10000 / Number of rows in table

PCOUNT = 0 / Size of special data area

GCOUNT = 1 / One data group (required keyword)

TFIELDS = 3 / Number of fields in each row

TTYPE1 = 'DATA_QUALITY' / Quality flag value of the photon

TFORM1 = '1B ' / Data format of the field: 1-byte integer

TTYPE2 = 'XPOS ' / Axis 1: pixel coordinate of the photon

TFORM2 = '1I ' / Data format of column 2: I*2 integer

TLMIN2 = 1 / Lower limit of axis in column 2

TLMAX2 = 2048 / Upper limit of axis in column 2

TTYPE3 = 'YPOS ' / Axis 2: pixel coordinate of the photon

TFORM3 = '1I ' / Data format of column 3: I*2 integer

TLMIN3 = 1 / Lower limit of axis in column 3

TLMAX3 = 2048 / Upper limit of axis in column 3

MJDOB3 = 44258.7845612 / MJD at start of observation

COMMENT The following keywords define the primary coordinate description.

TCRP2 = 1024.5 / Axis 1: reference pixel coordinate

TCRP3 = -1023.5 / Axis 2: reference point pixel coordinate

TPC2_2 = 1 / Transformation matrix element

TPC2_3 = -0.004 / Transformation matrix element

TPC3_2 = -0.002 / Transformation matrix element

TPC3_3 = 1 / Transformation matrix element

TCDE2 = -0.005 / Axis 1: scale

TCDE3 = 0.005 / Axis 2: scale

TCTY2 = 'GLON-COE' / Axis 1: conic equal area projection

TCTY3 = 'GLAT-COE' / Axis 2: conic equal area projection

TPV3_1 = -25.0 / Conic mid-latitude

TCRV2 = 90.0 / Axis 1: galactic longitude at reference point

TCRV3 = -25.0 / Axis 2: galactic latitude at reference point

COMMENT The following keywords define the secondary coordinate description.

TCRP2A = 1024.5 / Axis 1A: reference pixel coordinate

TCRP3A = -1023.5 / Axis 2A: reference point pixel coordinate

TP2_2A = 1 / Transformation matrix element

TP2_3A = -0.004 / Transformation matrix element

TP3_2A = -0.002 / Transformation matrix element

TP3_3A = 1 / Transformation matrix element

TCDE2A = -0.005 / Axis 1A: scale

TCDE3A = 0.005 / Axis 2A: scale

TCTY2A = 'ELON-COE' / Axis 1A: conic equal area projection

TCTY3A = 'ELAT-COE' / Axis 2A: conic equal area projection

TV3_1A = -25.0 / Conic mid-latitude

TCRV2A = -7.0300934 / Axis 1A: ecliptic longitude at reference point

TCRV3A = 34.8474143 / Axis 2A: ecliptic latitude at reference point

LONP3A = 6.3839706 / Native longitude of ecliptic pole

LATP3A = 29.8114400 / Ecliptic latitude of native pole

RADE3A = 'FK5 ' / Mean IAU 1984 ecliptic coordinates

EQUI3A = 2000.0 / Coordinate epoch

END

 

 

Table 11: Third example FITS header.

123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789

NAXIS = 2 / 2-dimensional image

NAXIS1 = 181 / x axis (fastest)

NAXIS2 = 91 / y axis (2nd fastest)

CRPIX1 = 226.0 / Pixel coordinate of reference point

CRPIX2 = 46.0 / Pixel coordinate of reference point

CDELT1 = -1.0 / x-scale

CDELT2 = 1.0 / y-scale

CTYPE1 = 'GLON-CAR' / Plate carree projection

CTYPE2 = 'GLAT-CAR' / Plate carree projection

CRVAL1 = 30.0 / Galactic longitude at reference point

CRVAL2 = 35.0 / Galactic latitude at reference point

Observe that the image spans $180\hbox{$^\circ$ }$ in native longitude and $90\hbox{$^\circ$ }$ in native latitude and that the reference pixel lies outside the image. In fact, the reference pixel is located so that the native longitude runs from $45\hbox{$^\circ$ }$ to $225\hbox{$^\circ$ }$ and hence the image lies partly inside and partly outside the normal range of native longitude, [ $-180\hbox{$^\circ$ }$, $180\hbox{$^\circ$ }$].

In fact, as might be expected, this makes no difference to the computation of celestial coordinates. For example, in computing the celestial coordinates of pixel (1,1) we readily find from Eqs. (83) and (84) that the native coordinates are $(\phi,\theta) = (225\hbox{$^\circ$ },-45\hbox{$^\circ$ })$. The fact that $\phi$ exceeds $180\hbox{$^\circ$ }$ becomes irrelevant once Eqs. (2) are applied since the trigonometric functions do not distinguish between $\phi = 225\hbox{$^\circ$ }$ and $\phi = -135\hbox{$^\circ$ }$. The latter value is the appropriate one to use if the cylinder of projection is considered to be "rolled out'' over multiple cycles. Consequently the correct galactic coordinates are obtained.

The problem only surfaces when we come to draw a coordinate grid on the image. A meridian of longitude, for example, is traced by computing the pixel coordinates for each of a succession of latitudes along the segment of the meridian that crosses the image. As usual, in computing pixel coordinates, the celestial coordinates are first converted to native coordinates by applying Eqs. (5), and the native longitude will be returned in the normal range [ $-180\hbox{$^\circ$ },180\hbox{$^\circ$ }$]. Consequently, in those parts of the image where $\phi~>~180\hbox{$^\circ$ }$ the pixel coordinates computed will correspond to the point at $\phi~-~360\hbox{$^\circ$ }$, i.e. in the part of the principle cycle of the cylindrical projection outside the image.

In principle it is possible to account for this, at least in specific cases, particularly for the cylindrical projections which are somewhat unusual in this regard. In practice, however, it is difficult to devise a general solution, especially when similar problems may arise for projection types where it is not desirable to track $\phi$ outside its normal range. For example, consider the case where a Hammer-Aitoff projection is used to represent the whole sky; since its boundary is curved there will be out-of-bounds areas in the corner of the image. Normally a grid drawing routine can detect these by checking whether the inverse projection equations return a value for $\phi$ outside its normal range. It may thus determine the proper boundary of the projection and deal with the discontinuity that arises when a grid line passes through it.

How then should the header have been written? Note that the problem exists at the lowest level of the coordinate description, in the conversion between (x,y) and $(\phi,\theta)$, and the solution must be found at this level. The problem arose from a particular property of cylindrical projections in that they have $x~\propto~\phi$. We must use this same property, which we might call "$\phi$-translation similarity'', to recast the coordinate description into a more manageable form. $\phi$-translation similarity simply means that changing the origin of $\phi$ corresponds to shifting the image in the x-direction. In other words, we can transfer the reference point of the projection from its current location to another location along the native equator without having to regrid the image. The fact that the PC i_ja matrix is unity in this example makes this task a little simpler than otherwise.

Note first that because the image straddles $\phi = 180\hbox{$^\circ$ }$ we can't simply reset CRPIX1 so as to shift the reference point to $\phi = -360\hbox{$^\circ$ }$; the image would then straddle $\phi = -180\hbox{$^\circ$ }$, which is no improvement. In this example it is convenient and sufficient to shift the reference point to $(\phi,\theta) = (180\hbox{$^\circ$ },0\hbox{$^\circ$ })$, which corresponds to pixel coordinate (p₁,p₂) = (46.0,46.0). Hence we need to reset CRPIX1 to 46.0 and adjust the keywords which define the celestial coordinate system. The reader may readily verify that the galactic coordinates of the new reference point are $(\ell,b) = (210\hbox{$^\circ$ },-35\hbox{$^\circ$ })$ and whereas the old, implied value of LONPOLE was $0\hbox{$^\circ$ }$ when $\delta_0~>~0$, now that $\delta_0~<~0$ its new implied value is $180\hbox{$^\circ$ }$, and this is correct. However, we will set it explicitly anyway. The keywords to be changed are therefore

$$\begin{eqnarray*}& & \hbox{{\tt CRPIX1}} = 46.0 , \\ & & \hbox{{\tt CRVAL1}} =... ...$^\circ$ }, \\ & & \hbox{{\tt LONPOLE}} = 180\hbox{$^\circ$ }. \end{eqnarray*}$$

What if the PC i_ja matrix was not unity? The problem of determining the pixel coordinates where $(\phi,\theta) = (-180\hbox{$^\circ$ },0\hbox{$^\circ$ })$ would have presented little extra difficulty, although in general CRPIX2 would also need to be changed. On reflection it may come as a surprise that changing CRPIX j like this does not fundamentally alter the linear transformation. However it may readily be verified that the only effect of changing CRPIX j is to change the origin of the (x,y) coordinates.

   
7.4 Header construction examples

This section considers the more difficult problem faced by FITS writers; that of constructing world coordinate system headers.

Example 1 is contrived to illustrate the general methods used in constructing a spherical coordinate representation. Paying homage to Claudius Ptolemy, it actually constructs a terrestrial coordinate grid for the Mediterranean region as seen from space.

Example 2 constructs headers for the infra-red dust maps produced by Schlegel et al. (1998) who regridded data from the COBE/DIRBE and IRAS/ISSA surveys onto two zenithal equal area projections.

Example 3 considers the case of long-slit optical spectroscopy. It is concerned in particular with solving the problem of producing three world coordinate elements for a data array of only two dimensions.

Example 4 constructs a dual coordinate description for an image of the moon and considers the problem of producing consistent scales for each. It also suggests an extension to deal with the rings of the planet Saturn.

	Figure 36: The Earth from 2230 km above Cairo looking directly towards Athens.
Open with DEXTER

   
7.4.1 Header construction example 1

Our first example of FITS header construction concerns satellite photography of the Earth. Figure 36 shows a part of the world which probably would have been recognizable to Claudius Ptolemy.

A satellite 2230 km immediately above Cairo aims its digital camera directly towards Athens and adjusts the orientation and focal length to include Cairo in the field near the bottom edge of the frame. The $2048 \times 2048$ pixel CCD detector array employed by the camera is centered on its optical axis so Athens is at pixel coordinate (1024.5,1024.5). Cairo - at the satellite's nadir - is later found to be at (681.67,60.12). The task is to overlay a terrestrial coordinate grid on this image.

For the sake of simplicity we will assume pin-hole camera optics. Figure 37 identifies the geometry as that of a tilted, near-sided zenithal perspective projection with the nadir at the reference point. The point of projection P corresponds to the pin-hole of the camera and the tilted plane of projection is parallel to the camera's focal plane F. The image in the focal plane is simply scaled (and inverted) with respect to that on the plane of projection. Clearly a tilted AZP projection is a better match to the geometry than SZP in this instance.

The satellite altitude of 2230 km is 0.350 Earth radii and since the projection is near-sided we may immediately write $\mu = -1.350$.

Determination of $\gamma$ requires knowledge of the coordinates of Cairo $(31\hbox{$.\!\!^\circ$ }15~{\rm E}, 30\hbox{$.\!\!^\circ$ }03~{\rm N})$ and Athens $(23\hbox{$.\!\!^\circ$ }44~{\rm E}, 38\hbox{$.\!\!^\circ$ }00~{\rm N})$. From spherical trigonometry we may deduce that the angular separation between the two cities is $10\hbox{$.\!\!^\circ$ }2072$ and also that Athens is on a bearing $36\hbox{$.\!\!^\circ$ }6252$ west of north from Cairo.

Now since Cairo is at the native pole of the projection the angular separation between the two cities is just their difference in native latitude, $90\hbox{$^\circ$ }- \theta_{\rm A}$. Using the sine rule in triangle PAO in Fig. 37 we obtain

$$90\hbox{$^\circ$ }- \theta = -\gamma - \sin^{-1}(\mu \sin\gamma) .$$ (194)

This equation, which takes account of the fact that $\mu$ is negative, may be solved iteratively to obtain $\gamma = 25\hbox{$.\!\!^\circ$ }8458$.

	Figure 37: The geometry of Fig. 36. Cairo at C is at the reference point of the projection with Athens at A. The camera with aperture at P and focal plane F is not drawn to scale, though the rest of the diagram is. Note that $\mu < 0$ so $-\mu$, as marked, is a positive value.
Open with DEXTER

The obliqueness of this view of the Earth, occasioned by the fact that the image is not oriented along a meridian, is handled partly via LONPOLE a and partly as a bulk image rotation via PC i_ja. Note that these are not complementary; they produce distinct effects since the tilted AZP projection does not have point symmetry. Figure 37 shows the situation - the generating sphere with Cairo at the native pole must be oriented so that Athens is at native longitude $180\hbox{$^\circ$ }$. The native longitude of the terrestrial pole (substituting for the celestial pole in this example) is offset from this by the bearing of Athens from Cairo, i.e. $\phi_{\rm p} = 180\hbox{$^\circ$ }- 36\hbox{$.\!\!^\circ$ }6252 = 143\hbox{$.\!\!^\circ$ }3748$.

Since Athens is at native longitude $180\hbox{$^\circ$ }$ its x-coordinate, like Cairo's, must be zero. The inequality of the corresponding pixel coordinates must have arisen by the satellite rotating the camera about its optical axis thereby rotating the CCD detector in the focal plane. The angle of this bulk rotation is readily deduced from the pixel coordinates of Athens and Cairo, $\arg(j_{\rm A}-j_{\rm C}, i_{\rm A}-i_{\rm C}) = 19\hbox{$.\!\!^\circ$ }570$. This rotation is to be applied via PC i_ja. The direction of rotation is completely determined by the requirement that Athens' x-coordinate be zero, regardless of the sign of CDELT ia, and the resulting PC i_ja rotation matrix is shown below.

The reference pixel coordinates, CRPIX ja, were measured from the image and the reference coordinates for Cairo, CRVAL ia, were obtained from an atlas, so the last remaining unknowns are the scales, CDELT ia. From previous calculations we know that Athens is at $(\phi,\theta) = (180\hbox{$^\circ$ }, 79\hbox{$.\!\!^\circ$ }7928)$ and we may apply Eqs. (20) and (21) to obtain $(x,y) = (0,8\hbox{$.\!\!^\circ$ }7424)$. The distance in pixel coordinates between the two cities is readily found to be 1023.5 so the y-scale must be $8\hbox{$.\!\!^\circ$ }7424 / 1023.5 = 0\hbox{$.\!\!^\circ$ }008542$ per pixel.

The x-scale cannot be determined like this since both cities have the same x-coordinate. However, the x-, and y-scales must be equal because the focal plane is parallel to the plane of projection and ray-tracing through the pin-hole therefore results in an isotropic change of scale. Do not confuse this with the statement made in Sect. 5.1.1 that with $\gamma \neq 0$the projection is not scaled true at the reference point; this refers to the differential scale between (x,y) and $(\phi,\theta)$.

Though the image in the focal plane is inverted through the pin-hole, thus indicating a negative scale, we can assume that the camera compensated by reading out the CCD array in reverse order. Thus for the image of Fig. 36 we make both of the CDELT ia positive in order to have east to the right as befits a sphere seen from the outside. Putting this all together we have

$$\begin{eqnarray*}&& \hbox{{\tt NAXIS}} = 2 , \\ && \hbox{{\tt NAXIS1}} = 2048 ... ... \hbox{{\tt WCSNAME}} = \hbox{{\tt 'Terrestrial coordinates'}} . \end{eqnarray*}$$

This example was simplified by the fact that the nadir was known and was included in the field of view. This probably would not be the case in a more realistic example where there may also be some uncertainty in the target coordinates and the value of $\mu$. In such cases the mapping parameters would have to be determined by least squares via the identification of an adequate number of surface features with known coordinates.

   
7.4.2 Header construction example 2

The following example comes from the $4096 \times 4096$ pixel infrared dust maps produced by Schlegel et al. (SFD, 1998). The authors chose to regrid data from the COBE/DIRBE and IRAS/ISSA maps onto two zenithal equal area (ZEA) projections centered on the galactic poles. The projection formula given in their Appendix C expressed in terms of standard 1-relative FITS pixel coordinates (p₁,p₂) is

$$\begin{eqnarray*}p_1 - 1 & = & \hphantom{-n}2048\sqrt{1 - n\sin b}~\cos\ell + 20... ...\\ p_2 - 1 & = & -n 2048\sqrt{1 - n\sin b}~\sin\ell + 2047.5 , \end{eqnarray*}$$

where n = +1 for the north Galactic pole and n = -1 for the south Galactic pole maps. Now for ZEA from Eqs. (1), (12), (13), and (69)

$$\begin{eqnarray*}\hbox{{\tt CDELT1}}~ (p_1 - \hbox{{\tt CRPIX1}}) & = & \hphant... ...{2} \frac{180\hbox{$^\circ$ }}{\pi}\sqrt{1-\sin\theta}\cos\phi . \end{eqnarray*}$$

If we take the north Galactic pole case first, the SFD equations may be rewritten as

$$\begin{eqnarray*}p_1 - 2048.5 & = & -2048 \sqrt{1 - \sin b}~\sin(\ell-90\hbox{$^... ....5 & = & -2048 \sqrt{1 - \sin b}~\cos(\ell-90\hbox{$^\circ$ }) . \end{eqnarray*}$$

By inspection of the two sets of equations we must have

$$\begin{eqnarray*}& & \hbox{{\tt NAXIS}} = 2 , \\ & & \hbox{{\tt NAXIS1}} = 4096... ..., \\ & & \ell = \phi + 90\hbox{$^\circ$ }, \\ & & b = \theta . \end{eqnarray*}$$

Now, writing $\ell$ and b in place of $\alpha$ and $\delta$ in Eqs. (2), we have to determine $\ell_{\rm p}$, $b_{\rm p}$, and $\phi_{\rm p}$ to give $(\ell,b)$ in terms of $(\phi,\theta)$. This is easy because we know $b_{\rm p}=90\hbox{$^\circ$ }$ and the simple special case Eqs. (3) apply, so

$$\begin{eqnarray*}\ell = \phi + (\ell_{\rm p} - \phi_{\rm p} - 180\hbox{$^\circ$ }) . \end{eqnarray*}$$

We have a degree of freedom here since only $\ell_{\rm p}-\phi_{\rm p}$ is determined. It's best to let $\phi_{\rm p}$ take its default value of $0\hbox{$^\circ$ }$for zenithal projections with the celestial pole at the native pole (it's $180\hbox{$^\circ$ }$ otherwise), so we must have

$$\begin{eqnarray*}& &\hbox{{\tt CRVAL1}} = 270\hbox{$^\circ$ }, \\ & &\hbox{{\t... ... \\ & &\hbox{{\tt LONPOLE}} = \hphantom{-9} 0\hbox{$^\circ$ }. \end{eqnarray*}$$

Although LONPOLE assumes its default value here it would of course be wise to write it explicitly into the header. The procedure for the south Galactic pole case is similar. The SFD equations may be rewritten

$$\begin{eqnarray*}p_1 - 2048.5 & = & 2048 \sqrt{1 - \sin(-b)}~\sin(90\hbox{$^\cir... ...5 & = & 2048 \sqrt{1 - \sin(-b)}~\cos(90\hbox{$^\circ$ }-\ell) , \end{eqnarray*}$$

whence

$$\begin{eqnarray*}& &\hbox{{\tt CRPIX1}} = 2048.5 , \\ & &\hbox{{\tt CRPIX2}} =... ...\ & & \ell = 90\hbox{$^\circ$ }- \phi , \\ & & b = -\theta . \end{eqnarray*}$$

Equations (4) apply for $b_{\rm p}=-90\hbox{$^\circ$ }$ so

$$\begin{eqnarray*}\ell & = & (\ell_{\rm p} + \phi_{\rm p}) - \phi . \end{eqnarray*}$$

Again we let $\phi_{\rm p}$ take its default value of $180\hbox{$^\circ$ }$, so

$$\begin{eqnarray*}&& \hbox{{\tt CRVAL1}} = 270\hbox{$^\circ$ }, \\ && \hbox{{\t... ...{$^\circ$ }, \\ && \hbox{{\tt LONPOLE}} = 180\hbox{$^\circ$ }. \end{eqnarray*}$$

It's generally easier to interpret coordinate headers than to construct them so it's essential after formulating the header to test it at a few points and make sure that it works as expected. Note that this sort of translation exercise wouldn't be necessary if the formalism of this paper was used right from the start, i.e. in the regridding operation used to produce the maps.

   
7.4.3 Header construction example 3

Consider now the coordinate description for the two-dimensional image formed by a long slit spectrograph. We assume that the wavelength axis of length 1024 and dispersion $\Delta\lambda$ nm/pixel corresponds to the p₁ pixel coordinate, and the 2048 pixel spatial axis corresponds to p₂. The slit is centered on equatorial coordinates $(\alpha _0,\delta _0)$ and oriented at position angle $\rho$ measured such that when $\rho = 0$ the first spectrum is northwards. We will assume that the telescope and spectrograph optics are such that the distance along the slit is in proportion to true angular distance on the sky with a separation between pixels of $\sigma$ arcsec. We do not consider curvature of the slit here, that is a distortion of the sort to be handled by the methods of Paper IV.

Equiscaling along the length of the slit together with the one-dimensional nature of the spatial geometry indicate that any projection could be used that is equiscaled along a great circle projected as a straight line. Many projections satisfy this criterion. However, in practice the zenithal equidistant (ARC) projection is the most convenient choice.

The one-dimensional spatial geometry may at first seem problematical, with each point along the slit having a different  $(\alpha ,\delta )$. However, this is easily handled by introducing a third, degenerate axis and introducing a rotation. The rotation in position angle may be handled via the linear transformation matrix. Moreover, since we've opted for a zenithal projection, bearing in mind the discussion of Sect. 7.1, the rotation can also be handled via $\phi_{\rm p}$ (i.e. LONPOLE). We will demonstrate both methods.

Use of $\phi_{\rm p}$ is perhaps more straightforward, and the header may be written without further ado:

$$\begin{eqnarray*}& & \hbox{{\tt NAXIS}} = 3 , \\ & & \hbox{{\tt NAXIS1}} = 1024... ...ta_0 , \\ & & \hbox{{\tt LONPOLE}} = 90\hbox{$^\circ$ }+ \rho . \end{eqnarray*}$$

LONPOLE is the only card which requires explanation. Since no rotation is introduced by the linear transformation, the slit, which coincides with the p₂-axis, maps directly to the x-axis. However, because CDELT2 is negative, the two axes run in opposite directions and, given that when $\rho = 0$ the p₂-axis runs from north to south, x must run from south to north. Referring to the left side of Fig. 3 for zenithal projections, we see that when $\rho = 0$, the north celestial pole must be at $\phi = 90\hbox{$^\circ$ }$. Since position angle increases in an anticlockwise direction (i.e. north through east) in the plane of the sky, the celestial pole must rotate clockwise as $\rho$ increases. Thus, as $\rho$ increases so must $\phi_{\rm p}$, hence $\phi_{\rm p} = 90\hbox{$^\circ$ }+ \rho$.

To verify this representation we will test it with $(\alpha_0,\delta_0) = (150\hbox{$^\circ$ },-35\hbox{$^\circ$ })$, $\rho = 30\hbox{$^\circ$ }$, and $\sigma = 2''$ for pixel coordinate (1,1,1). The corresponding (s,x,y) coordinates are $(0,0\hbox{$.\!\!^\circ$ }5686111,0\hbox{$^\circ$ })$. The wavelength is therefore $\lambda_0$. From Eqs. (15), (14), and (67) the native longitude and latitude are $(\phi,\theta) = (90\hbox{$^\circ$ },89\hbox{$.\!\!^\circ$ }4313889)$, $\phi$is always $\pm 90\hbox{$^\circ$ }$, and the value of $\theta$ corresponds to an offset of 2047'' from the center of the slit as it should. From Eqs. (2) the right ascension and declination are $(\alpha,\delta) = (150\hbox{$.\!\!^\circ$ }3450039,-34\hbox{$.\!\!^\circ$ }5070794)$, which are in the correct quadrant.

Most optical telescopes are better described by the TAN projection. In this case the only difference in the header would be

$$\begin{eqnarray*}\hbox{{\tt CTYPE2}} & = & \hbox{{\tt 'RA---TAN'}} , \\ \hbox{{\tt CTYPE3}} & = & \hbox{{\tt 'DEC--TAN'}} . \end{eqnarray*}$$

Verifying it with the same values as before yields the same (s,x,y) and thus the same wavelength. From Eqs. (15), (14), and (54) we get $(\phi,\theta) = (90\hbox{$^\circ$ },89\hbox{$.\!\!^\circ$ }4314076)$, the offset of 2046 $.\!\!^{\prime\prime}$933 differing by only 67 mas. The right ascension and declination are $(\alpha,\delta) = (150\hbox{$.\!\!^\circ$ }3449926,-34\hbox{$.\!\!^\circ$ }5070956)$.

As an alternative, the original header could also have been written with the CTYPE i interchanged, so that the slit, still coincident with the p₂-axis, maps directly to the y-axis. The header would be as above but with the following changes:

$$\begin{eqnarray*}& & \hbox{{\tt CTYPE2}} = \hbox{{\tt 'DEC--ARC'}} , \\ & & \hb... ...a_0 , \\ & & \hbox{{\tt LONPOLE}} = 180\hbox{$^\circ$ }+ \rho . \end{eqnarray*}$$

LONPOLE differs because the slit now runs along the y-axis rather than the x-axis; the negative sign on CDELT2 is still needed to make y run counter to p₂.

There are several practical ways to rewrite the header using the PC i_j or CD i_j matrices. Looking first at the CD i_j form we have

$$\begin{eqnarray*}\left( \begin{array}{c} s \\ x \\ y \end{array} \right) & =... ...t CRPIX2}} \\ p_3 - \hbox{{\tt CRPIX3}} \end{array} \right) . \end{eqnarray*}$$

Since the third axis is degenerate, with $\hbox{{\tt CRPIX3}} = 1$, we have $p_3 - \hbox{{\tt CRPIX3}} = 0$, so the values of CD1_3, CD2_3, and CD3_3 are irrelevant. Moreover, CD1_2, CD2_1, and CD3_1 are all zero, hence

 

$$\hbox{{\tt CD1\_1}} \ (p_1 - \hbox{{\tt CRPIX1}}) ,$$

$$\hbox{{\tt CD2\_2}} \ (p_2 - \hbox{{\tt CRPIX2}}) ,$$ (195)

$$\hbox{{\tt CD3\_2}} \ (p_2 - \hbox{{\tt CRPIX2}}) .$$

Effectively this provides (x,y) coordinates along the slit via the parametric equations of a line, where $p_2 - \hbox{{\tt CRPIX2}}$ is the distance parameter.

As before, the main problem is to determine the correct angle of rotation. For zenithal projections $\phi _{\rm p} = 180\hbox {$^\circ $ }$ by default, and referring to the left side of Fig. 3 we see that this corresponds to the y-axis. Thus when $\rho = 0$ we require a rotation of $90\hbox{$^\circ$ }$ to transform the p₂-axis onto the y-axis. For $\rho > 0$ we need to rotate further so the angle of rotation is $90\hbox{$^\circ$ }+ \rho$.

Therefore, CDELT i and LONPOLE in the original header must be replaced with

$$\begin{eqnarray*}\hbox{{\tt CD1\_1}} & = & \Delta\lambda , \\ \hbox{{\tt CD2\_... ...\hbox{{\tt CD2\_3}} & = & 1 , \\ \hbox{{\tt CD3\_3}} & = & 1 . \end{eqnarray*}$$

The negative sign on CD3_2 is associated with the rotation, not the scale. Note that, although the value of CD3_3 is irrelevant, it must be set non-zero otherwise the zero defaults for CD i_j would make the third column zero, thereby producing a singular matrix. Likewise, in this example and similarly below, we set CD2_3 non-zero to prevent the second row of the matrix becoming zero for values of $\rho$ such that $\cos(90\hbox{$^\circ$ }+\rho) = 0$.

The PC i_j matrix formulation is similar but allows more flexibility in the way the scale is handled:

1.

Since CDELT i defaults to unity, omitting it allows PC i_j to duplicate the functionality of CD i_j. However, since PC i_j defaults to the unit matrix rather than zero, there is no need to set PC3_3 specifically, hence

$$\begin{eqnarray*}\hbox{{\tt PC1\_1}} & = & \Delta\lambda , \\ \hbox{{\tt PC2\_... ...irc$ }+\rho) (\sigma/3600) , \\ \hbox{{\tt PC2\_3}} & = & 1 . \end{eqnarray*}$$

2.

Equation (195), as a simple scaling relation, suggests the use of CDELT i. However, PC3_2 must be non-zero since x and y both depend on p₂. Setting it to unity and allowing the remaining PC i_j to default to the unit matrix we have

$$\begin{eqnarray*}& & \hbox{{\tt PC3\_2}} = 1 , \\ & & \hbox{{\tt CDELT1}} = \D... ...{\tt CDELT3}} = - \sin(90\hbox{$^\circ$ }+\rho) (\sigma/3600) . \end{eqnarray*}$$

3.

The "orthodox'' method is to encode the rotation and scale separately in PC i_j and CDELT i:

$$\begin{eqnarray*}& &\hbox{{\tt PC2\_2}} = \hphantom{-}\cos(90\hbox{$^\circ$ }+\r... ...2}} = \sigma/3600 , \\ & &\hbox{{\tt CDELT3}} = \sigma/3600 . \end{eqnarray*}$$

Arguably this is more amenable to human interpretation, especially if thoughtful comments are added.

The above six methods should all be regarded as equally legitimate. In fact, there are infinitely many ways to encode this header, though most would disguise the essential simplicity of the geometry.

   
7.4.4 Header construction example 4

Thompson (1999) has applied the methods of this paper to the definition of solar coordinates for a variety of coordinate systems. As the final header construction example we will consider a specific coordinate description for another solar system body, the Moon.

A short exposure plate taken at Sydney Observatory at 1:00 am on 15 February 1957 AEST (GMT +1000) of the full Moon near lunar perigee has been digitized and converted to a $4096 \times 4096$ pixel image in FITS format. The scale is $1\hbox{$^{\prime\prime}$ }$ per pixel with the moon centered in the image, and it is desired to construct a dual coordinate description. The first system is geocentric apparent equatorial coordinates in a gnomonic projection. The second is selenographic coordinates in a zenithal perspective projection attached to the surface of the Moon.

The ephemeris gives the geocentric apparent right ascension of the Moon at the time as $09^{\rm h}41^{\rm m}13\hbox{$.\!\!^{\rm s}$ }1$ and declination as $+08\hbox{$^\circ$ }34\hbox{$^\prime$ }26\hbox{$^{\prime\prime}$ }$. Diurnal parallax would have caused the Moon to appear slightly offset from this position as seen from Sydney Observatory, but to a good approximation the coordinate systems may be defined with the center of the Moon at the stated geocentric coordinates. The image was digitized in the normal orientation with north up and east to the left so the header for the primary description is straightforward:

$$\begin{eqnarray*}& & \hbox{{\tt NAXIS}} = 2 , \\ & & \hbox{{\tt NAXIS1}} = 409... ...= 180.0 , \\ & & \hbox{{\tt RADESYS}} = \hbox{{\tt 'GAPPT'}} . \end{eqnarray*}$$

If any rotation had been introduced it could have been corrected for in the linear transformation matrix.

	Figure 38: Geometry of header construction example 4 (not to scale). The observer is at O, a distance of $\vert\mu\vert$ Moon radii from the center of the Moon. The Moon is projected as a near-sided zenithal perspective projection ($\mu < -1$), and the celestial sphere as a gnomonic projection. An alternative plane of projection is shown.
Open with DEXTER

 

 

Table 12: Summary of new coordinate keywords introduced in this paper; see also Table 3 for alternate types used in binary tables.

Keyword	Type	Sect.	Use	Status	Comments
LONPOLE a	floating	2.2	coordinate rotation	new	Longitude in the native coordinate system of the
					celestial system's north pole.
					Default: $0\hbox{$^\circ$ }$ if $\delta_0\ge\theta_0$, $180\hbox{$^\circ$ }$ otherwise.
LATPOLE a	floating	2.4	coordinate rotation	new	Latitude in the native coordinate system of the
					celestial system's north pole, or equivalently, the
					latitude in the celestial coordinate system of the
					native system's north pole.
					Default: $90\hbox{$^\circ$ }$, unless $(\theta_0, \delta_0, \phi_{\rm p} - \phi_0) = (0,0, \pm 90\hbox{$^\circ$ })$
					in which case it has no default.
RADESYS a	character	3.1	frame of reference	new	Reference frame of equatorial and ecliptic
					coordinates; recognized values are given in Table 2.
					Default: FK4 if EQUINOX a < 1984.0, FK5 if
					$\ge\ 1984.0$, or ICRS if EQUINOX a is not given.
EQUINOX a	floating	3.1	coordinate epoch	new	Epoch of the mean equator and equinox in years;
					Besselian if RADESYS a is FK4 or FK4-NO-E,
					Julian if FK5; not applicable to ICRS or GAPPT.
					Default: EPOCH if given, else 1950.0 if RADESYS a is FK4
					or FK4-NO-E, or 2000.0 if FK5.
EPOCH	floating	3.1	coordinate epoch	deprecated	Replaced by EQUINOX a.
MJD-OBS	floating	3.1	time of observation	new	Modified Julian Date (JD - 2 400 000.5) of observation.
					Default: DATE-OBS if given, else no default.

The first step in constructing the secondary description is to determine the distance of the observer from the Moon. The ephemeris gives the equatorial horizontal parallax as $61\hbox{$^\prime$ }29\hbox{$.\!\!^{\prime\prime}$ }3$, corresponding to a distance between the centers of the Earth and Moon of 55.91 Earth radii. However, the observer may be closer or further away by up to one Earth radius depending on location and time of day and this 2% effect is deemed worthy of correction. At Sydney Observatory (longitude $10^{\rm h}04^{\rm m}49\hbox{$.\!\!^{\rm s}$ }2$, latitude $-33\hbox{$^\circ$ }51\hbox{$^\prime$ }41\hbox{$^{\prime\prime}$ }$) the Greenwich apparent sidereal time was $00^{\rm h}51^{\rm m}43\hbox{$.\!\!^{\rm s}$ }6$. Thus the apparent right ascension and declination of the zenith was $10^{\rm h}56^{\rm m}32\hbox{$.\!\!^{\rm s}$ }8$, $-33\hbox{$^\circ$ }51\hbox{$^\prime$ }41\hbox{$^{\prime\prime}$ }$. The vector dot product then gives the distance correction as 0.69 Earth radii. Using the ratio of the Earth and Moon radii of 3.670 the corrected distance of 55.22 Earth radii indicates that the distance parameter for the zenithal perspective projection is $\mu = 202.64$Moon radii.

We now need to determine the correct relative scale. Figure 38 shows the geometry of the two projections where we note, by analogy with Fig. 4, that $\mu < -1$. From the diagram we have

$$\beta = \tan^{-1} \left( \frac{\sin\gamma}{\vert\mu\vert - \cos\gamma} \right) \cdot$$ (196)

The figure shows that this equation is not influenced by the choice of the plane of projection. Evaluating the derivative  ${\rm d}\beta/{\rm d}\gamma$at $\gamma = 0$ we find the relative scale is $1/(\vert\mu\vert - 1)$.

The ephemeris records that the selenographic coordinates of the Earth at the time were $(\ell,b) = (+0\hbox{$.\!\!^\circ$ }26,+6\hbox{$.\!\!^\circ$ }45)$ and the position angle of the Moon's axis was $C = 19\hbox{$.\!\!^\circ$ }03$. However, since the Earth subtends an angle of over $2\hbox{$^\circ$ }$ in the lunar sky, topocentric optical libration - the correction for the observatory's location - is significant. Application of the correction formulæ derived by Atkinson (1951) gives the selenographic coordinates of Sydney Observatory as $(\ell',b') = (-0\hbox{$.\!\!^\circ$ }23,5\hbox{$.\!\!^\circ$ }89)$ and $C' = 18\hbox{$.\!\!^\circ$ }94$. Since the image was taken in the normal orientation and we have a zenithal projection it is convenient to account for the position angle by setting $\phi_{\rm p} = 180\hbox{$^\circ$ }- C'$. Adopting keyword values of SELN and SELT for selenographic coordinates we may write

$$\begin{eqnarray*}& & \hbox{{\tt CRPIX1S}} = 2048.5 , \\ & & \hbox{{\tt CRPIX2S... ...box{{\tt WCSNAMES}} = \hbox{{\tt 'SELENOGRAPHIC COORDINATES'}} . \end{eqnarray*}$$

We have used the letter S to denote the alternate coordinate system simply to demonstrate that there is no requirement to start with A. A WCSNAME a keyword, defined in Paper I, is used to identify the coordinate system. Note that selenographic coordinates are right-handed so that selenographic longitude increases towards the west on the celestial sphere as seen from the Earth and this is handled by setting CDELT1S positive.

As an extension of this example, a FITS header with three coordinate systems might be constructed for an image of Saturn; a celestial grid for the background, a saturnographic system for the surface of the planet, and a third system for its rings. The rings might be described by a zenithal equidistant (ARC) projection with associated linear transformation matrix set to match the oblique viewing angle.

 

 

Table 13: Summary of projection codes, default values of $\phi _0$ and $\theta _0$, special properties, full name, and required parameters. Where a projection has a special property other than being conformal, equiareal, or equidistant a reference is given to the section where this is discussed.

FITS			Special		Projection parameters associated with latitude$^\dag$ axis i
Code	$\phi _0$	$\theta _0$	properties	Projection name	PV i_0a	PV i_1a	PV i_2a	PV i_3a	PV i_ma
AZP	$0\hbox{$^\circ$ }$	$90\hbox{$^\circ$ }$	Sect. 5.1.1	Zenithal perspective		$\mu$	$\gamma$
SZP	$0\hbox{$^\circ$ }$	$90\hbox{$^\circ$ }$	Sect. 5.1.2	Slant zenithal perspective		$\mu$	$\phi_{\rm c}$	$\theta_{\rm c}$
TAN	$0\hbox{$^\circ$ }$	$90\hbox{$^\circ$ }$	Sect. 5.1.3	Gnomonic
STG	$0\hbox{$^\circ$ }$	$90\hbox{$^\circ$ }$	Sect. 5.1.4, Conformal	Stereographic
SIN	$0\hbox{$^\circ$ }$	$90\hbox{$^\circ$ }$	Sect. 5.1.5	Slant orthographic		$\xi$	$\eta$
ARC	$0\hbox{$^\circ$ }$	$90\hbox{$^\circ$ }$	Equidistant	Zenithal equidistant
ZPN	$0\hbox{$^\circ$ }$	$90\hbox{$^\circ$ }$		Zenithal polynomial	P0	P1	P2	P3	$\ldots P_{20}$
ZEA	$0\hbox{$^\circ$ }$	$90\hbox{$^\circ$ }$	Equiareal	Zenithal equal-area
AIR	$0\hbox{$^\circ$ }$	$90\hbox{$^\circ$ }$	Sect. 5.1.9	Airy		$\theta_{\rm b}$
CYP	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$		Cylindrical perspective		$\mu$	$\lambda$
CEA	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$	Equiareal	Cylindrical equal area		$\lambda$
CAR	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$	Equidistant	Plate carrée
MER	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$	Conformal	Mercator
SFL	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$	Equiareal	Sanson-Flamsteed
PAR	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$	Equiareal	Parabolic
MOL	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$	Equiareal	Mollweide
AIT	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$	Equiareal	Hammer-Aitoff
COP	$0\hbox{$^\circ$ }$	$\theta_{\rm a}$		Conic perspective		$\theta_{\rm a}$	$\eta$
COE	$0\hbox{$^\circ$ }$	$\theta_{\rm a}$	Equiareal	Conic equal-area		$\theta_{\rm a}$	$\eta$
COD	$0\hbox{$^\circ$ }$	$\theta_{\rm a}$	Equidistant	Conic equidistant		$\theta_{\rm a}$	$\eta$
COO	$0\hbox{$^\circ$ }$	$\theta_{\rm a}$	Conformal	Conic orthomorphic		$\theta_{\rm a}$	$\eta$
BON	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$	Equiareal	Bonne's equal area		$\theta_1$
PCO	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$		Polyconic
TSC	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$	Sect. 5.6.1	Tangential Spherical Cube
CSC	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$	Sect. 5.6	COBE Quadrilateralized Spherical Cube
QSC	$0\hbox{$^\circ$ }$	$0\hbox{$^\circ$ }$	Sect. 5.6	Quadrilateralized Spherical Cube

$^\dag$ Parameters PV i_0a, PV i_1a, and PV i_2a associated with longitude axis i implement user-specified values of $(\phi _0,\theta _0)$ as discussed in Sect. 2.5, and PV i_3a, and PV i_4a encapsulate the values of LONPOLE a and LATPOLE a respectively, as discussed in Sect. 2.6.

   
7.5 Realization

Calabretta (1995) has written, and made available under a GNU license, a package of routines, WCSLIB, which implements all projections and coordinate conversions defined here. It contains independent C and Fortran libraries.

The Fortran library includes a routine, PGSBOX, which is based on PGPLOT and uses WCSLIB to draw general curvilinear coordinate axes. It was used to generate Fig. 36.

   
8 Summary

We have developed a flexible method for associating celestial coordinates with a FITS image and implemented it for all projections likely to be of use in astronomy. It should be a relatively simple matter to add new projections should the need ever arise.

The FITS-header keywords defined in this paper are summarized in Table 12 with the 3-letter projection codes listed in Table 13. The column labeled $\theta _0$ gives the native latitude of the reference point in degrees ( $\phi_0 = 0$ always) and the parameters associated with each projection are listed in the nomenclature of Sect. 5.

Acknowledgements

The authors would particularly like to thank Patrick Wallace (U.K. Starlink) for his helpful comments which included the text for some of the paragraphs of Sect. 3.1 and the suggestion of the zenithal polynomial projection. Also, Steve Allen (UCO/Lick Observatory) who provided valuable feedback and suggestions over the long period of this work's development. We are grateful for the assistance provided by Bob Hanisch (Space Telescope Science Institute) in bringing the manuscript to its final form. We also thank William Pence, Arnold Rots, and Lorella Angelini (NASA Goddard Space Flight Center) for contributing the text of Sect. 4.

The authors thank the following for comments and suggestions: Rick Balsano, David Berry, Emmanuel Bertin, Jeff Bloch, Peter Bunclark, Julian Daniels, Lindsey Davis, Rick Ebert, Douglas Finkbeiner, Immanuel Freedman, Brian Glendenning, Bill Gray, Dave Green, Mike Kesteven, Doug Mink, Jonathan McDowell, François Ochsenbein, Fred Patt, Tim Pearson, William Thompson, Doug Tody, Stephen Walton, and Don Wells.

Many others contributed indirectly via requests for clarification of particular points leading potentially to improvements in the wording or formalism.

Mark Calabretta would also like to thank the following for feedback on WCSLIB: Klaus Banse, David Barnes, David Berry, Wim Brouw, Lindsey Davis, Rick Ebert, Ken Ebisawa, Brian Glendenning, Neil Killeen, Thomas A. McGlynn, Doug Mink, Clive Page, Ray Plante, Keith A. Scollick, and Peter Teuben.

Coastline data for Fig. 36 are from the "CIA World DataBank II'' map database, which is in the public domain. Dave Pape's plain-text version was obtained from http://www.evl.uic.edu/pape/data/WDB/ and this site also provides references to the original, highly compressed binary encodings.

The Australia Telescope is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO.

The National Radio Astronomy Observatory is a facility of the (U.S.) National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

   
Appendix A: Projection aliases

Table A.1 provides a list of aliases which have been used by cartographers for special cases of the projections described in Sect. 5.

  

Table A.1: Projection aliases.

   
Appendix B: Mathematical methods

B.1 Coordinate rotation with matrices

The coordinate rotations represented in Eqs. (2) or (5) may be represented by a matrix multiplication of a vector of direction cosines. The matrix and its inverse (which is simply the transpose) may be precomputed and applied repetitively to a variety of coordinates, improving performance. Thus, we have

$$\left( \begin{array}{c} l \\ m \\ n \end{array} \right) = \l... ... l^{\prime} \\ m^{\prime} \\ n^{\prime} \end{array} \right) ,$$ (B.1)

where

$$\begin{eqnarray*}(l',m',n') & = & (\cos\delta\cos\alpha,\cos\delta\sin\alpha, \s... ...m p} \cos\delta_{\rm p} , \\ r_{33} & = & \sin\delta_{\rm p} . \end{eqnarray*}$$

The inverse equation is

$$\left( \begin{array}{c} l^{\prime} \\ m^{\prime} \\ n^{\prime... ...ht) \left( \begin{array}{c} l \\ m \\ n \end{array} \right) .$$ (B.2)

B.2 Iterative solution

Iterative methods are required for the inversion of several of the projections described in this paper. One, Mollweide's, even requires solution of a transcendental equation for the forward equations. However, these do not give rise to any particular difficulties.

On the other hand, it sometimes happens that one pixel and one celestial coordinate element is known and it is required to find the others; this typically arises when plotting graticules on image displays. Although analytical solutions exist for a few special cases, iterative methods must be used in the general case. If, say, p₂ and $\alpha$ are known, one would compute pixel coordinate as a function of $\delta$ and determine the value which gave p₂. The unknown pixel coordinate elements would be obtained in the process.

This prescription glosses over many complications, however. All bounded projections may give rise to discontinuities in the graph of p₂ versus $\delta$ (to continue the above example), for example where the meridian through $\alpha$ crosses the $\phi = \pm180\hbox{$^\circ$ }$ boundary in cylindrical, conic and other projections. Even worse, if the meridian traverses a pole represented as a finite line segment then p₂ may become multivalued at a particular value of $\delta$. The derivative $\partial p_2/\partial\delta$will also usually be discontinuous at the point of discontinuity, and it should be remembered that some projections such as the quad-cubes may also have discontinuous derivatives at points within their boundaries.

We will not attempt to resolve these difficulties here but simply note that WCSLIB (Calabretta, 1995) implements a solution.

   
Appendix C: The slant orthographic projection

The slant orthographic or generalized SIN projection derives from the basic interferometer equation (e.g. Thompson et al. 1986). The phase term in the Fourier exponent is

$$\wp = \vec{(e - e_0) \cdot B} ,$$ (C.1)

where $\vec{e_0}$ and $\vec{e}$ are unit vectors pointing towards the field center and a point in the field, $\vec{B}$ is the baseline vector, and we measure the phase $\wp$ in rotations so that we don't need to carry factors of $2\pi$. We can write

$$\wp = p_u u + p_v v + p_w w ,$$ (C.2)

where (u,v,w) are components of the baseline vector in a right-handed coordinate system with the w-axis pointing from the geocenter towards the source and the u-axis lying in the equatorial plane, and

$$\begin{array}{rccl} p_u & = & & \cos \theta \sin \phi , \\ ... ...eta \cos \phi , \\ p_w & = & & \sin \theta - 1 , \end{array}$$ (C.3)

are the coordinates of $\vec{(e - e_0)}$, where $(\phi,\theta)$ are the longitude and latitude of $\vec{e}$ in the spherical coordinate system with the pole towards $\vec{e_0}$ and origin of longitude towards negative v, as required by Fig. 3. Now, for a planar array we may write

n_u u + n_v v + n_w w = 0 (C.4)

where (n_u,n_v,n_w) are the direction cosines of the normal to the plane. Using this to eliminate w from Eq. (C.2) we have

$$\wp = \left[p_u - \frac{n_u}{n_w} p_w\right] u + \left[p_v - \frac{n_v}{n_w} p_w\right] v\cdot$$ (C.5)

Being the Fourier conjugate variables, the quantities in brackets become the Cartesian coordinates, in radians, in the plane of the synthesized map. Writing

$$\begin{array}{rcl} \xi & = & n_u / n_w , \\ \eta & = & n_v / n_w , \end{array}$$ (C.6)

Eqs. (61) and (62) are then readily derived from Eqs. (C.3) and (C.5). For 12-hour synthesis by an east-west interferometer the baselines all lie in the Earth's equatorial plane whence $(\xi,\eta) = (0, \cot\delta_0)$, where $\delta_0$ is the declination of the field center. For a "snapshot'' observation by an array such as the VLA, Cornwell & Perley (1992) give $(\xi,\eta) = (-\tan Z\sin\chi, \tan Z\cos\chi)$, where Z is the zenith angle and $\chi$ is the parallactic angle of the field center at the time of the observation.

In synthesizing a map a phase shift may be applied to the visibility data in order to translate the field center. If the shift applied is

$$\Delta\wp = q_u u + q_v v + q_w w$$ (C.7)

where (q_u,q_v,q_w) is constant then Eq. (C.2) becomes

$$\wp = (p_u - q_u) u + (p_v - q_v) v + (p_w - q_w) w ,$$ (C.8)

whence Eq. (C.5) becomes

$$\begin{array}{rcl} \wp & = & [(p_u - q_u) - \xi (p_w - q_w)] u ~ \\ & & +[(p_v - q_v) - \eta(p_w - q_w)] v . \end{array}$$ (C.9)

Equations (61) and (62) become

$$\hphantom{-}\frac{180\hbox{$^\circ$ }}{\pi}\left[~ \cos\theta \si... ...ta - 1)~ \right] - \frac{180\hbox{$^\circ$ }}{\pi}\left[q_u - \xi q_w \right] ,$$

$$- \frac{180\hbox{$^\circ$ }}{\pi}\left[~ \cos\theta \cos\phi - \e... ...a - 1)~ \right] - \frac{180\hbox{$^\circ$ }}{\pi}\left[q_v - \eta q_w \right] ,$$ (C.10)

from which on comparison with Eqs. (61) and (62) we see that the field center is shifted by

$$(\Delta x, \Delta y) = \frac{180\hbox{$^\circ$ }}{\pi}(q_u - \xi q_w, ~ q_v - \eta q_w) .$$ (C.11)

The shift is applied to the coordinate reference pixel.

References

 
• Atkinson, R. d'E. 1951, MNRAS, 111, 448 In the text NASA ADS  
• Blaauw, A., Gum, C. S., Pawsey, J. L., & Westerhout, G. 1960, MNRAS, 121, 123 NASA ADS  
• Brouw, W. N. 1971, Ph.D. Thesis, Sterrenwacht, Leiden In the text In the text  
• Bunclark, P., & Rots, A. 1996, Precise Re-definition of DATE-OBS Keyword Encompassing the Millennium, located at http://www.cv.nrao.edu/fits/documents/standards/ as year2000.txt  
• Calabretta, M. R. 1992, Spherical Map Projections, AIPS++ Memo No. 107, The AIPS++ Consortium, c/o National Radio Astronomy Observatory, Charlottesville, Virginia, Dec 15. In the text http://aips2.nrao.edu/docs/memos/107.ps  
• Calabretta, M. R. 1995, WCSLIB, version 2, located at In the text ftp://ftp.atnf.csiro.au/pub/software/wcslib/wcs- lib.tar.gz  
• Calabretta, M. R., et al. 2002, Representations of distortions in FITS world coordinate systems, A&A, in preparation In the text  
• Chan, F. K., & O'Neill, E. M. 1975, Feasibility Study of a Quadrilateralized Spherical Cube Earth Data Base, EPRF Technical Report 2 - 75 (CSC), Computer Sciences Corporation (Silver Springs, Maryland) In the text  
• Cornwell, T. J., & Perley, R. A. 1992, A&A, 261, 353 In the text NASA ADS  
• Cotton, W. D., Tody, D., & Pence, W. D. 1995, A&AS, 113, 159 In the text NASA ADS  
• Evenden, G. I. 1991, Cartographic Projection Procedures for the UNIX Environment - A User's Manual, Open-File Report 90-284, United States Department of the Interior Geological Survey. Available as OF90-284.ps.gz at In the text ftp://kai.er.usgs.gov/pub/PROJ.4/ with accompanying code and related documents. The documents are also cached at http://www.cv.nrao.edu/fits/wcs/USGS/  
• Górski, K. M., & Hivon, E. F. 1999, HEALPix homepage, located at In the text http://www.eso.org/science/healpix  
• Geldzahler, B., & Schlesinger, B. M. 1997, A User's Guide for the Flexible Image Transport System (FITS), NASA/OSSA Office of Standards and Technology, Greenbelt, Maryland, version 4.0 (April 1997). Located via In the text http://www.gsfc.nasa.gov/astro/fits/fits_home.html  
• Greisen, E. W. 1983, Non-linear Coordinate Systems in AIPS, AIPS Memo No. 27, National Radio Astronomy Observatory, Charlottesville, Virginia. Revised January 1993, located at In the text http://www.cv.nrao.edu/fits/wcs/ as aips27.ps.Z  
• Greisen, E. W. 1986, Additional Non-linear Coordinates, AIPS Memo No. 46, National Radio Astronomy Observatory, Charlottesville, Virginia, May 20, available as aips46.ps.Z from In the text http://www.cv.nrao.edu/fits/wcs/  
• Greisen, E. W., Valdes, F. G., Calabretta, M. R., & Allen, S. L. 2002, Representations of spectral coordinates in FITS, A&A, in preparation In the text  
• Greisen, E. W., & Calabretta, M. R. 2002, A&A, 395, 1061 (Paper I) In the text  
• Greisen, E. W., & Harten, R. H. 1981, A&AS, 44, 371 In the text NASA ADS  
• Hanisch, R. J., & Wells, D. C. 1988, World Coordinate Systems Representations Within the FITS Format, draft based on notes from a meeting sponsored by the National Aeronautics and Space Administration, Code EZ, held in Charlottesville, Virginia in January 1988, located at In the text http://www.cv.nrao.edu/fits/wcs/ as wcs88.ps.Z  
• Hanisch, R. J., Farris, A., Greisen, E. W., et al. 2001, A&A, 376, 359 In the text NASA ADS  
• Hohenkerk, C. Y., Yallop, B. D., Smith, C. A., & Sinclair, A. T. 1992, Explanatory Supplement to the Astronomical Almanac, ed. P. K. Seidelmann, Chap. 3 (University Science Books, Mill Valley, California) In the text  
• IAU, 1998, Proceedings of the Twenty-Third General Assembly, Kyoto 1997, Transactions of the IAU, vol. XXIII-B In the text  
• Jones, M. H. 1993, Observatory, 113, 213 NASA ADS  
• Kellaway, G. P. 1946, Map Projections (Methuen & Co. Ltd., London, England) In the text  
• O'Neill, E. M., & Laubscher, R. E. 1976, Extended Studies of the Quadrilateralized Spherical Cube Earth Data Base, NEPRF Technical Report 3 - 76 (CSC) Computer Sciences Corporation (Silver Springs, Maryland) In the text  
• Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525 In the text NASA ADS  
• Snyder, J. P. 1993, Flattening the Earth (University of Chicago Press, Chicago and London) In the text  
• Tegmark, M. 1996, ApJ, 470, L81 In the text NASA ADS  
• Thompson, A. R., Moran, J. W., & Swenson, G. W. 1986, Interferometry and Synthesis in Radio Astronomy (J. Wiley, New York, USA) In the text  
• Thompson, W. 1999, Coordinate Systems for Solar Image Data (proposed). Located at In the text http://orpheus.nascom.na- sa.gov/~thompson/papers/coordinates.ps  
• Wells, D. C., Greisen, E. W., & Harten, R. H. 1981, A&AS, 44, 363 In the text NASA ADS