A&A 395, 10771122 (2002)
DOI: 10.1051/00046361:20021327
M. R. Calabretta^{1}  E. W. Greisen^{2}
1  Australia Telescope National Facility,
PO Box 76,
Epping, NSW 1710, Australia
2  National Radio Astronomy Observatory,
PO Box O,
Socorro, NM 878010387, 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
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 twodimensional 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 multistep 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 substeps, 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
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
to intermediate world coordinates
becomes
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 IAUendorsed 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 nonlinear 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 FITSreading programs. Section 7 discusses the concepts presented here, including worked examples of header interpretation and construction.
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, . The equations for the transformation depend on the particular projection and this will be specified via the CTYPE ia keyword. Paper I defined "43'' form for such purposes; the rightmost threecharacters 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 , 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 threeletter 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.
The last step in the chain of transformations shown in Fig. 1 is the spherical rotation from native coordinates, , to celestial^{} coordinates . 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 . 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.
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  
p_{j}  Pixel coordinates  
r_{j}  Reference pixel coordinates  CRPIX ja 
m_{ij}  Linear transformation matrix  CD i_ja or PC i_ja 
s_{i}  Coordinate scales  CDELT ia 
x_{i}  Intermediate world coordinates (in general)  
(x,y)  Projection plane coordinates  
Native longitude and latitude  
Celestial longitude and latitude  
Native longitude and latitude of the fiducial point  PV i_1a, PV i_2a  
Celestial longitude and latitude of the fiducial point  CRVAL ia  
Native longitude and latitude of the celestial pole  LONPOLE a ( = PV i_3a), LATPOLE a, ( = PV i_4a)  
Celestial longitude and latitude of the native pole  
Inverse tangent function that returns the correct quadrant 
Accordingly we specify only that a fiducial celestial coordinate pair given by the CRVAL ia will be associated with a fiducial native coordinate pair defined explicitly for each projection. For example, zenithal projections all have , while cylindricals have . The AIPS convention has been honored here as far as practicable by constructing the projection equations so that 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 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 chosen at some midlatitude 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, , specified by the new FITS keyword
LONPOLE a (floatingvalued). 
The default value of LONPOLE a will be 0 for or for . 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 (unless ) since . In cylindrical projections, where , the default value for LONPOLE a is 0 for , but it is for .
Since differs for different projections it is apparent that the relationship between and the required Euler angles also differs.
For zenithal projections,
so the CRVAL ia
specify the celestial coordinates of the native pole,
i.e.
.
There is a
simple relationship between the Euler angles for consecutive rotations about
the Z, X, and Zaxes and
,
and
;
the ZXZ Euler angles are
.
Given this close correspondence it is convenient to
write the Euler angle transformation formulæ directly in terms of
,
and
:
Projections such as the cylindricals and conics for which are handled by providing formulae to compute from whence the above equations may be used.
Given that
are the celestial coordinates of the point
with native coordinates
,
Eqs. (6) and
(7) may be inverted to obtain
Note, however, that if , as is usual, then when LONPOLE a ( ) takes its default value of 0 or (depending on ) then any combination of and produces a valid argument to the in Eq. (8), and at least one of the solutions is valid.
LATPOLE a (floatingvalued) 
is chosen. It is acceptable to set LATPOLE a to a number greater than to choose the northerly solution (the default if LATPOLE a is omitted), or a number less than to select the southern solution.
These rules governing the application of Eqs. (810) are certainly the most complex of this formalism. FITS writers are well advised to check the values of , , , , and against them to ensure their validity.
In Sect. 2.2 we formally decoupled from the reference point and associated it with . One implication of this is that it should be possible to allow to be userspecifiable. 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 and 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 by introducing an implied offset (x_{0},y_{0}) which is computed for 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 projectionspecific defaults given in Sect. 5.
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 for each is at right ascension , declination (marked). The two sets of FITS header cards differ only in their CTYPE ia keyword values. The nonoblique graticule could be obtained in the current system by setting .  
Open with DEXTER 
A change of coordinate system may be effected in a straightforward way if the transformation from the original system, , to the new system, , and its inverse are known. The new coordinates of the , namely , are obtained simply by transforming . To obtain , 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 . As a check, compute via Eq. (8) and verify that .
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 , and illustrated in Fig. 18. Figure 1 shows that while the transformation from (x,y) to may be linear (in fact identical), there still remains the nonlinear transformation from to . 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 RACAR, DECCAR 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 right ascension and declination and with . It is evident that since the plate carrée has , a nonoblique graticule may only be obtained by setting with . 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.
As mentioned in Sect. 2.1, Paper I defined "43'' 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, , the left half simply identifies the celestial system with which 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 righthanded and that the pole be at latitude . 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), , and . This leads us to associate x, , and on the one hand, and y, , and on the other. This simply means, for example, that if CTYPE3 = 'GLONAIT', then the third element of the intermediate world coordinate calculated via Eq. (1) corresponds to what we have been calling the xcoordinate in the plane of projection, the association being between 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.
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 (charactervalued) 
to specify the particular system. Recognized values are given in Table 2. Apart from FK4NOE these keywords are applicable to ecliptic as well as equatorial coordinates.
RADESYS a  Definition 
ICRS  International Celestial Reference System 
FK5  mean place, 
new (IAU 1984) system  
FK4  mean place, 
old (BessellNewcomb) system  
FK4NOE  mean place, 
old system but without eterms  
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 (floatingvalued), 
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 FK4NOE, 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 . 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 DATEOBS keyword may be used for this purpose. However, to provide a more convenient specification we here introduce the new keyword
MJDOBS (floatingvalued), 
that provides DATEOBS as a Modified Julian Date ( ) but is identical to it in all other respects. MJDOBS does not have a version code since there can only be one time of observation. Following the year2000 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 socalled "eterms'' which amount to 343 milliarcsec. Strictly speaking, therefore, a map obtained from, say, a radio synthesis telescope, should be regarded as FK4NOE 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 eterms are effectively corrected to first order only.
The randomgroups 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 onepixel 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, MJDOBS, 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.
Paper I defined the coordinate keywords used to describe a multidimensional 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.
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:
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  MJDOBS  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_{1} and p_{2} 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.
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, SansonFlamsteed, HammerAitoff 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 infinitevalued 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 specialcase 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 SansonFlamsteed 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. , and (right) with the intersection of the equator and prime meridian at the reference point i.e. .  
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,
Conformality is a property which applies to points in the plane of projection which are locally distortionfree. 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. . 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_{1},p_{2}) pixel coordinates.
Figure 4: (Left) geometry of the zenithal perspective projections, the point of projection at P is 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_{1}pixel coordinate increasing to the right with p_{2} increasing upwards, i.e. a righthanded system. This means that the (x,y) coordinates must be lefthanded as shown in Fig. 3. Note, however, that the approximation cannot hold unless 1) , and hence , do actually map to the reference point (Sect. 2.2), 2) 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_{0}, known as the "radius of the generating sphere''. However, in this work r_{0} is set explicitly to in order to maintain backwards compatibility with the AIPS convention. This effectively sets the circumference of the generating sphere to 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_{0} 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, , celestial coordinates, , 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.
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
(11) 
Figure 5: Alternate geometries of slant zenithal perspective projections with and : (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 and 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
spherical radii from the
center of the sphere with
increasing in the direction away from the
plane of projection, then it is straightforward to show that
(17) 
(18) 
A nearsided perspective projection may be obtained with . This correctly represents the image of a sphere, such as a planet, when viewed from a distance times the planetary radius. The coordinates of the reference point may be expressed in planetary longitude and latitude, . 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 nonzero lookangle, , 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:
(22) 
=  (23) 
=  (24) 
=  (25) 
=  (26) 
=  (27) 
R=  (28) 
Figure 6: Slant zenithal perspective ( AZP) projection with and which corresponds to the lefthand side of Fig. 5. The reference point of the corresponding SZP projection is marked at .  
Open with DEXTER 
With
the projection is not scaled true at the reference point.
In fact the x scale is correct but the y scale is magnified by
,
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 ,
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
(29) 
(30) 
(31) 
Definition of the perimeter of the projection is more complicated for the
slant projection than the orthogonal case. As before, for
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
(32) 
In general, for
,
the projection is bounded, otherwise
it is unbounded. However, the latitude of divergence is now a function of
:
=  (34) 
=  (35) 
The FITS keywords PV i_1a and PV i_2a, attached to latitude coordinate i, will be used to specify, respectively, in spherical radii and in degrees, both with default value 0.
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 and , which corresponds to the righthand side of Fig. 5. The reference point of the corresponding AZP projection is marked at .  
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 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 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
,
then
a = X'^{2} + Y'^{2} + 1 ,  (39) 
b = X'(XX') + Y'(YY') ,  (40) 
c = (X  X')^{2} + (Y  Y')^{2}  1 ,  (41) 
(42)  
(43) 
=  (45) 
=  (46) 
=  (47) 
For
the sphere is divided into two unequal segments that are
projected in superposition. The limb is defined by computing the native
latitude
as a function of
(49)  
(50)  
(51)  
(52) 
(53) 
Figure 8: Gnomonic ( TAN) projection; diverges at .  
Open with DEXTER 
The zenithal perspective projection with ,
the gnomonic
projection^{}, is widely
used in optical astronomy and was given its own code within the AIPS
convention, namely TAN^{}. For ,
Eq. (16) reduces to
Figure 9: Stereographic ( STG) projection; diverges at .  
Open with DEXTER 
The stereographic projection, the second important special case of the
zenithal perspective projection defined by the AIPS convention, has ,
for which Eq. (16) becomes
=  (56) 
= 
(57) 
(58) 
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.
The zenithal perspective projection with , 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, , north and south sides begin to overlap at ; (bottom left) , i.e. ; (bottom right) projection appropriate for an eastwest array observing at , , .  
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
(59) 
(60) 
=  (63) 
=  (64) 
It can be shown that the slant orthographic projection is equivalent to an
orthographic projection centered at
which has been
stretched in the
direction by a factor of
.
The projection equations may be inverted using
Eqs. (38) and (44) except that
Eq. (43) is replaced with
(65) 
(66) 
Figure 11: Zenithal equidistant ( ARC) projection; no limits.  
Open with DEXTER 
Some nonperspective 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
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
(68) 
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 nthdegree expansion of one of the standard zenithal projections.
If P_{0} is nonzero the native pole is mapped to an open circle centered on the reference point as illustrated in Fig. 12. In other words, is not at (x,y) = (0,0), which in fact lies outside the boundary of the projection. However, we do not dismiss 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 (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, , PV i_20a, attached to latitude coordinate i, all of which have default values of zero.
Figure 13: Zenithal equal area projection ( ZEA); no limits.  
Open with DEXTER 
Lambert's zenithal equalarea projection illustrated in Fig. 13 is
constructed by defining
so that the area enclosed by the native parallel
at latitude
in the plane of projection is equal to the area of the
corresponding spherical cap. It may be generated using
(70) 
The Airy projection^{} minimizes the error for the region within latitude
(Evenden 1991). It is defined by
= 
= 
Figure 14: Airy projection ( AIR) with ; diverges at .  
Open with DEXTER 
The FITS keyword PV i_1a, attached to latitude coordinate i, will be used to specify in degrees with a default of . This projection is illustrated in Fig. 14.
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
(72) 
(73) 
Figure 15 illustrates the geometry for the construction of
cylindrical perspective projections. The sphere is projected onto a cylinder
of radius
spherical radii from points in the equatorial plane of the
native system at a distance
spherical radii measured from the center of
the sphere in the direction opposite the projected surface. The cylinder
intersects the sphere at latitudes
.
It is
straightforward to show that
x=  (74) 
y=  (75) 
=  (76) 
=  (77) 
(78) 
The case with
is covered by the class of cylindrical equal area
projections. No other specialcases 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
,
giving
Figure 16: Gall's stereographic projection, CYP with , ; no limits.  
Open with DEXTER 
The cylindrical equal area projection is the special case of the cylindrical
perspective projection with
.
It is conformal at latitudes
where
.
The formulæ are
x=  (79) 
y=  (80) 
=x ,  (81) 
=  (82) 
Figure 17: Lambert's equal area projection, CEA with ; no limits.  
Open with DEXTER 
Lambert's^{} equal area projection, the case with , 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.
The equator and all meridians are correctly scaled in the plate carrée
projection^{}, whose main virtue is that of simplicity. Its formulæ are
Figure 18: The plate carrée projection ( CAR); no limits.  
Open with DEXTER 
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
(85) 
x=  (86) 
y=  (87) 
(88) 
Refer to Sect. 6.1.4 for a discussion of the usage of MER in AIPS.
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
(89) 
Figure 19: Mercator's projection ( MER); diverges at .  
Open with DEXTER 
Bonne's projection (Sect. 5.5.1) reduces to the pseudocylindrical
SansonFlamsteed^{} projection
when
.
Parallels are equispaced and projected at their true
length which makes it an equal area projection. The formulæ are
x=  (90) 
y=  (91) 
=  (92) 
=y .  (93) 
Figure 20: SansonFlamsteed 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: HammerAitoff projection ( AIT); no limits.  
Open with DEXTER 
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
x=  (94) 
y=  (95) 
=  (96) 
=  (97) 
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
x=  (98) 
y=  (99) 
(100) 
=  (101) 
=  
(102) 
The HammerAitoff^{} 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 allsky projections.
Figure 24: Construction of conic perspective projections ( COP) and the resulting graticules; (left) twostandard projection with , ; (right) onestandard projection with . Both projections have and this accounts for their similarity. Both diverge at .  
Open with DEXTER 
The formulæ for the projection and its inverse are derived in Greisen
(1986) and Calabretta (1992) among others. They are
x=  (103) 
y=  (104) 
(105) 
=  (106) 
=  (107) 
Z=  (108) 
=  (109) 
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.
Twostandard conic projections are characterized by two latitudes,
and ,
whose parallels are projected at their true length.
In the conic perspective projection these are the latitudes at which the cone
intersects the sphere. Onestandard conic projections have
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 midway between the two standard parallels maps
to the reference point so we set
(110) 
=  (112) 
x=  (114) 
y=  (115) 
Figure 25: Conic equal area projection ( COE) with , and , also illustrating the inversion of southern hemisphere conics; no limits.  
Open with DEXTER 
The conics will be parameterized in FITS by
(given by
Eq. (111)) and
where
(118) 
As noted in Sect. 5, the zenithal projections are special cases of the conics with . Likewise, the cylindrical projections are conics with . However, we strongly advise against using conics in these cases for the reasons given previously. Nevertheless, the only formal requirement on and in the equations presented below is that .
Development of Colles' conic perspective projection is shown in
Fig. 24. The projection formulæ are
C=  (121) 
=  (122) 
Y_{0}=  (123) 
(124) 
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
Figure 26: Conic equidistant projection ( COD) with and ; 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
C=  (130) 
=  (131) 
Y_{0}=  (132) 
(133) 
(134)  
(135)  
(136) 
(137) 
The requirement for conformality of conic projections is
(138) 
C=  (139) 
=  (140) 
Y_{0}=  (141) 
=  (142) 
=  (143) 
(144) 
Figure 27: Conic orthomorphic projection ( COO) with and ; diverges at .  
Open with DEXTER 
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 subclass 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 ,
as shown in Fig. 24, we have
.
Since both are constructed with the native coordinate system origin at the
reference point we set
(145) 
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
for which
.
The projection is conformal at this
latitude and along the central meridian. The equations for Bonne's projection
become divergent for
and this special case is handled as the
SansonFlamsteed projection. The projection formulæ are
x=  (146) 
y=  (147) 
=  (148) 
=  (149) 
Y_{0}=  (150) 
=  (151) 
=  (152) 
=  (153) 
=  (154) 
Figure 28: Bonne's projection ( BON) with ; no limits.  
Open with DEXTER 
Each parallel in Hassler's polyconic projection is projected as an arc of a
circle of radius
at its true length,
,
and
correctly divided. The scale along the central meridian is true and
consequently the parallels are not concentric. The projection formulæ are
x=  (155) 
y=  (156) 
(158) 
Quadrilateralized spherical cube (quadcube) 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 quadcubes 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 closepacked arrangement rather than the simple rectangular arrangement supported by FITS.
The six faces of quadcube projections are numbered and laid out as
Face  
0  l  
1  
2  l  
3  m  l  
4  m  
5  n 
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
(159) 
=  (161) 
=  (162) 
While perspective quadcube 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
and
as
described above, then
x=  (163) 
y=  (164) 
=  (165) 
=  (166) 
=  (167) 
=  (168) 
=  (169) 
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
x=  (170) 
y=  (171) 
Figure 31: COBE quadrilateralized spherical cube projection ( CSC); no limits.  
Open with DEXTER 
=  (173) 
=  (174) 
X= 
Y= 
Given the face number, , and , the native coordinates 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
to
with
Eq. (175) and then transforming back to
with Eq. (172), i.e.
Measures of equalarea 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 % within the inner 64% of the face.
Figure 32: Quadrilateralized spherical cube projection ( QSC); no limits.  
Open with DEXTER 
O'Neill & Laubscher (1976) derived an exact expression for an equalarea 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
and must be reflected into the other
quadrants. This has the effect of making the projection nondifferentiable
along the diagonals as is evident in Fig. 32. To compute the
forward projection first identify the face and find
and
from Table 4. Then
(176) 
u=  (177) 
v=  (178) 
=  
S= 
(179) 
(180) 
(181) 
=  (182) 
=  (183) 
=  (184) 
=  (185) 
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.
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
(188) 
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:
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.
The "northcelestialpole'' projection defined by Greisen (1983) is
a special case of the new generalized SIN projection. The old header
cards
The TAN, ARC and STG projections defined by Greisen (1983, 1986) are directly equivalent to those defined here and no translation is required.
Special care is required in interpreting the AIT (HammerAitoff), GLS (SansonFlamsteed), 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 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 nonzero 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 within the AIPS convention. For AIT and MER (but not GLS), CDELT i must also be corrected for the scaling factors and incorporated into the AIPS projection equations.
Of the three projections only GLS is known to have been used with nonzero CRVAL i. Consequently we have renamed it as SFL as a warning that translation is required.
As mentioned in Sect. 6, FITS interpreters will need to recognize the AIPS convention virtually forever. It stands to reason, therefore, that if modern FITSwriters 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 FITSwriters 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
=  
=  (191) 
(192) 
Solutions for CROTA2 come in pairs separated by . The above formulæ give the solution which falls in the halfopen interval . The other solution is obtained by subtracting from CROTA2 and negating CDELT1 and CDELT2. While each solution is equally valid, if one makes and 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.
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 northsouth 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 nonoblique graticules are often plotted together on the same map. It wouldn't make sense to describe such a projection as being simultaneously oblique and nonoblique; 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, HammerAitoff, 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, , together with by setting as described in Sect. 2.5.
The first graticule, A, when compared to the nonoblique native coordinate graticules presented earlier for each projection, illustrates the effect of changing (and hence ). Comparison of graticules A and B shows that changing (and hence ) results in a simple change in origin of longitude. Graticules A, C, and D show the more interesting effect of varying (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 graticules plotted for the zenithal equal area ( ZEA) and conic equidistant ( COD) projections with parameters (A) , , ; (B) , , ; (C) , , ; (D) , , .  
Open with DEXTER 
Figure 34: Oblique ( ) graticules plotted for the HammerAitoff ( AIT) and COBE quadrilateralized spherical cube ( CSC) projections using the same parameters as Fig. 33.  
Open with DEXTER 
The projected meridians and parallels in Figs. 33 and 34 also serve to illustrate the distortions introduced by the various projections. In particular, the quadcube 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 . The same is true for spacecraft generated images of distant moons and planets. Allsky 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 HammerAitoff (AIT) projection is probably the most widely used for allsky maps. However, the SansonFlamsteed (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 . 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 "Australiasized'' regions of the sphere and at this they excel. Oblique forms of the zenithal equalarea 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 then a number of projections will serve provided also that 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 ( ), ZPN ( P_{0}=0, P_{1}=1), AIR ( ), CYP (), CEA (), COP (), COD (), COE (), and COO (). 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_{1},p_{2}). Thus, for example, the plate carrée projection (CAR) also serves as the more general equirectangular projection.
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.
123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789 
NAXIS = 4 / 4dimensional 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 = 'RATAN' / 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 = 'DECTAN' / 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 
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
for which the CRVAL i give equatorial
coordinates
in right ascension and
in declination. The equatorial north pole has a
longitude of
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
If we define the projection nonlinearity as the departure of from , then in this image it amounts to at the corners. However, in a image it quickly escalates to , sixty times larger. Comparison of for the southeast and northeast corners indicates the significant effect of grid convergence even in this moderate sized image. The two differ by , most of which is attributable to the term in converting longitude offsets to true angular distances. The effect of grid convergence is small for near the equator and large near the poles.
Figure 35 investigates the effect of projective nonlinearities for moderate field sizes for the commonly used zenithal projections. It shows the difference in between the various projections and the SIN projection as a function of native latitude . The SIN projection is used as reference since it always has less than the other zenithal projections. The difference for all projections exceeds 1 arcsec for values of less than 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.
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
.
The tile size of
pixels is approximately
,
and this
particular tile is situated immediately to the north of the central tile.
parameter  units  SE corner  NE corner  NW corner 
(p_{1},p_{2})  pixel  (1, 2)  (1, 512)  (511, 512) 
(p_{3},p_{4})  pixel  (1, 1)  (1, 1)  (196, 1) 
x  deg  
y  deg  
deg  
deg  
deg  
deg  
Velocity  ms^{1}  500000.00  500000.00  1890018.50 
Stokes  I  I  I 
Figure 35: in arcsec plotted versus 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)
The next step is to compute the native longitude and latitude. Like all standard conics, the conic equal area projection has two parameters, and , 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, are obtained from (x,y) using Eqs. (117) and (129) via the intermediaries C, Y_{0}, and given by Eqs. (125), (127), and (116). These in turn are expressed in terms of , , and given by Eqs. (128), (119), and (120). The calculations are shown in Table 8.
123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789 
NAXIS = 2 / 2dimensional image 
NAXIS1 = 2048 / x axis (fastest) 
NAXIS2 = 2048 / y axis (2nd fastest) 
MJDOBS = 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 / xscale 
CDELT2 = 0.005 / yscale 
CTYPE1 = 'GLONCOE' / Conic equal area projection 
CTYPE2 = 'GLATCOE' / Conic equal area projection 
PV2_1 = 25.0 / Conic midlatitude 
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 / xscale 
CDELT2A = 0.005 / yscale 
CTYPE1A = 'ELONCOE' / Conic equal area projection 
CTYPE2A = 'ELATCOE' / Conic equal area projection 
PV2_1A = 25.0 / Conic midlatitude 
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, , and ecliptic coordinates, . To do this we need to apply Eqs. (2) and in order to do that we need and .
Looking first at galactic coordinates, we have and . This conic projection has , and because and , the native and galactic coordinate systems must coincide to within an offset in longitude. However, it is not obvious what should be to produce this offset. Equation (8) has one valid solution, namely . The special case in point 2 in the usage notes for Eqs. (9) and (10) must therefore be used to obtain . 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 and the native longitude of the ecliptic pole is . It also specifies LATPOLEA as . In this case Eq. (8) has two valid solutions, , 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 would have selected the northerly solution anyway, but of course it is good practice to make the choice clear. The value of 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.
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 2D arrays in Col. 5 of the table and each row of the table contains a 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 2dimensional 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.
This example has been adapted from a reallife 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.
123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789 
XTENSION= 'BINTABLE' / Binary table extension 
BITPIX = 8 / 8bit bytes 
NAXIS = 2 / 2dimensional 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 ' / 2D 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 = 'GLONCOE' / Axis 1: conic equal area projection 
2CTY5 = 'GLATCOE' / Axis 2: conic equal area projection 
2PV5_1 = 25.0 / Conic midlatitude 
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 = 'ELONCOE' / Axis 1A: conic equal area projection 
2CTY5A = 'ELATCOE' / Axis 2A: conic equal area projection 
2PV5_1A = 25.0 / Conic midlatitude 
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 
123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789 
XTENSION= 'BINTABLE' / Binary table extension 
BITPIX = 8 / 8bit bytes 
NAXIS = 2 / 2dimensional 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: 1byte 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 = 'GLONCOE' / Axis 1: conic equal area projection 
TCTY3 = 'GLATCOE' / Axis 2: conic equal area projection 
TPV3_1 = 25.0 / Conic midlatitude 
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 = 'ELONCOE' / Axis 1A: conic equal area projection 
TCTY3A = 'ELATCOE' / Axis 2A: conic equal area projection 
TV3_1A = 25.0 / Conic midlatitude 
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 
123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789 
NAXIS = 2 / 2dimensional 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 / xscale 
CDELT2 = 1.0 / yscale 
CTYPE1 = 'GLONCAR' / Plate carree projection 
CTYPE2 = 'GLATCAR' / Plate carree projection 
CRVAL1 = 30.0 / Galactic longitude at reference point 
CRVAL2 = 35.0 / Galactic latitude at reference point 
Observe that the image spans in native longitude and 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 to and hence the image lies partly inside and partly outside the normal range of native longitude, [ , ].
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 . The fact that exceeds becomes irrelevant once Eqs. (2) are applied since the trigonometric functions do not distinguish between and . 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 [ ]. Consequently, in those parts of the image where the pixel coordinates computed will correspond to the point at , 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 outside its normal range. For example, consider the case where a HammerAitoff projection is used to represent the whole sky; since its boundary is curved there will be outofbounds 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 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 , and the solution must be found at this level. The problem arose from a particular property of cylindrical projections in that they have . We must use this same property, which we might call "translation similarity'', to recast the coordinate description into a more manageable form. translation similarity simply means that changing the origin of corresponds to shifting the image in the xdirection. 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
we can't simply
reset CRPIX1 so as to shift the reference point to
;
the image would then straddle
,
which is no improvement. In
this example it is convenient and sufficient to shift the reference point to
,
which corresponds to pixel coordinate
(p_{1},p_{2}) = (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
and whereas the old, implied value of
LONPOLE was
when
,
now that
its
new implied value is
,
and this is correct. However, we will set it
explicitly anyway. The keywords to be changed are therefore
What if the PC i_ja matrix was not unity? The problem of determining the pixel coordinates where 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.
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 infrared 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 longslit 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 
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 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 pinhole camera optics. Figure 37 identifies the geometry as that of a tilted, nearsided zenithal perspective projection with the nadir at the reference point. The point of projection P corresponds to the pinhole 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 nearsided we may immediately write .
Determination of requires knowledge of the coordinates of Cairo and Athens . From spherical trigonometry we may deduce that the angular separation between the two cities is and also that Athens is on a bearing 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,
.
Using the sine rule in triangle PAO in
Fig. 37 we obtain
(194) 
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 so , 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 . 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. .
Since Athens is at native longitude its xcoordinate, 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, . This rotation is to be applied via PC i_ja. The direction of rotation is completely determined by the requirement that Athens' xcoordinate 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 and we may apply Eqs. (20) and (21) to obtain . The distance in pixel coordinates between the two cities is readily found to be 1023.5 so the yscale must be per pixel.
The xscale cannot be determined like this since both cities have the same xcoordinate. However, the x, and yscales must be equal because the focal plane is parallel to the plane of projection and raytracing through the pinhole therefore results in an isotropic change of scale. Do not confuse this with the statement made in Sect. 5.1.1 that with the projection is not scaled true at the reference point; this refers to the differential scale between (x,y) and .
Though the image in the focal plane is inverted through the pinhole, 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
The following example comes from the
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 1relative
FITS pixel coordinates (p_{1},p_{2}) is
Consider now the coordinate description for the twodimensional image formed by a long slit spectrograph. We assume that the wavelength axis of length 1024 and dispersion nm/pixel corresponds to the p_{1} pixel coordinate, and the 2048 pixel spatial axis corresponds to p_{2}. The slit is centered on equatorial coordinates and oriented at position angle measured such that when 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 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 onedimensional 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 onedimensional spatial geometry may at first seem problematical, with each point along the slit having a different . 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 (i.e. LONPOLE). We will demonstrate both methods.
Use of
is perhaps more straightforward, and the header may be
written without further ado:
To verify this representation we will test it with , , and for pixel coordinate (1,1,1). The corresponding (s,x,y) coordinates are . The wavelength is therefore . From Eqs. (15), (14), and (67) the native longitude and latitude are , is always , and the value of corresponds to an offset of 2047'' from the center of the slit as it should. From Eqs. (2) the right ascension and declination are , 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
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_{2}axis, maps directly to the yaxis. The header would be as above but
with the following changes:
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
As before, the main problem is to determine the correct angle of rotation. For zenithal projections by default, and referring to the left side of Fig. 3 we see that this corresponds to the yaxis. Thus when we require a rotation of to transform the p_{2}axis onto the yaxis. For we need to rotate further so the angle of rotation is .
Therefore, CDELT i and LONPOLE in the original header must be
replaced with
The PC i_j matrix formulation is similar but allows more flexibility in the way the scale is handled:
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.
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 pixel image in FITS format. The scale is 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
and declination as
.
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:
Figure 38: Geometry of header construction example 4 (not to scale). The observer is at O, a distance of Moon radii from the center of the Moon. The Moon is projected as a nearsided zenithal perspective projection (), and the celestial sphere as a gnomonic projection. An alternative plane of projection is shown.  
Open with DEXTER 
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: if , 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: , unless  
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  
, 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 FK4NOE,  
Julian if FK5; not applicable to ICRS or GAPPT.  
Default: EPOCH if given, else 1950.0 if RADESYS a is FK4  
or FK4NOE, or 2000.0 if FK5.  
EPOCH  floating  3.1  coordinate epoch  deprecated  Replaced by EQUINOX a. 
MJDOBS  floating  3.1  time of observation  new  Modified Julian Date (JD  2 400 000.5) of observation. 
Default: DATEOBS 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 , 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 , latitude ) the Greenwich apparent sidereal time was . Thus the apparent right ascension and declination of the zenith was , . 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 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 .
From the diagram we have
(196) 
The ephemeris records that the selenographic coordinates of the Earth at the
time were
and the position angle of the
Moon's axis was
.
However, since the Earth subtends an angle of
over
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
and
.
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
.
Adopting keyword values of SELN and SELT for selenographic
coordinates we may write
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.
FITS  Special  Projection parameters associated with latitude axis i  
Code  properties  Projection name  PV i_0a  PV i_1a  PV i_2a  PV i_3a  PV i_ma  
AZP  Sect. 5.1.1  Zenithal perspective  
SZP  Sect. 5.1.2  Slant zenithal perspective  
TAN  Sect. 5.1.3  Gnomonic  
STG  Sect. 5.1.4, Conformal  Stereographic  
SIN  Sect. 5.1.5  Slant orthographic  
ARC  Equidistant  Zenithal equidistant  
ZPN  Zenithal polynomial  P_{0}  P_{1}  P_{2}  P_{3}  
ZEA  Equiareal  Zenithal equalarea  
AIR  Sect. 5.1.9  Airy  

Cylindrical perspective  
CEA  Equiareal  Cylindrical equal area  
CAR  Equidistant  Plate carrée  
MER  Conformal  Mercator  

Equiareal  SansonFlamsteed  
PAR  Equiareal  Parabolic  
MOL  Equiareal  Mollweide  
AIT  Equiareal  HammerAitoff  

Conic perspective  
COE  Equiareal  Conic equalarea  
COD  Equidistant  Conic equidistant  
COO  Conformal  Conic orthomorphic  

Equiareal  Bonne's equal area  
PCO  Polyconic  

Sect. 5.6.1  Tangential Spherical Cube  
CSC  Sect. 5.6  COBE Quadrilateralized Spherical Cube  
QSC  Sect. 5.6  Quadrilateralized Spherical Cube 
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.
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 FITSheader keywords defined in this paper are summarized in Table 12 with the 3letter projection codes listed in Table 13. The column labeled gives the native latitude of the reference point in degrees ( 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 plaintext 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.
Table A.1 provides a list of aliases which have been used by cartographers for special cases of the projections described in Sect. 5.
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
(B.1) 
(B.2) 
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_{2} and are known, one would compute pixel coordinate as a function of and determine the value which gave p_{2}. 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_{2} versus (to continue the above example), for example where the meridian through crosses the boundary in cylindrical, conic and other projections. Even worse, if the meridian traverses a pole represented as a finite line segment then p_{2} may become multivalued at a particular value of . The derivative will also usually be discontinuous at the point of discontinuity, and it should be remembered that some projections such as the quadcubes 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.
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
(C.1) 
n_{u} u + n_{v} v + n_{w} w = 0  (C.4) 
(C.6) 
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
(C.7) 
(C.8) 
(C.9) 
x=  
y=  (C.10) 
(C.11) 