7 Discussion

   
7.1 Oblique projections

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

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

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

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

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

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

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

   
7.2 Choice of projection

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

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

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

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

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

   
7.3 Header interpretation examples

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

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

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

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

   
7.3.1 Header interpretation example 1

 

 

Table 5: FITS header for example 1.

123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789

NAXIS = 4 / 4-dimensional cube

NAXIS1 = 512 / x axis (fastest)

NAXIS2 = 512 / y axis (2nd fastest)

NAXIS3 = 196 / z axis (planes)

NAXIS4 = 1 / Dummy to give a coordinate

CRPIX1 = 256 / Pixel coordinate of reference point

CDELT1 = -0.003 / 10.8 arcsec per pixel

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

CRVAL1 = 45.83 / RA at reference point

CUNIT1 = 'deg ' / Angles are degrees always

CRPIX2 = 257 / Pixel coordinate of reference point

CDELT2 = 0.003 / 10.8 arcsec per pixel

CTYPE2 = 'DEC--TAN' / Gnomonic projection

CRVAL2 = 63.57 / Dec at reference point

CUNIT2 = 'deg ' / Angles are degrees always

CRPIX3 = 1 / Pixel coordinate of reference point

CDELT3 = 7128.3 / Velocity increment

CTYPE3 = 'VELOCITY' / Each plane at a velocity

CRVAL3 = 500000.0 / Velocity in m/s

CUNIT3 = 'm/s ' / metres per second

CRPIX4 = 1 / Pixel coordinate of reference point

CDELT4 = 1 / Required here.

CTYPE4 = 'STOKES ' / Polarization

CRVAL4 = 1 / Unpolarized

CUNIT4 = ' ' / Conventional unitless = I pol

LONPOLE = 180 / Native longitude of celestial pole

RADESYS = 'FK5 ' / Mean IAU 1984 equatorial coordinates

EQUINOX = 2000.0 / Equator and equinox of J2000.0

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

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

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

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

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

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

which, on substitution, become

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

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

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

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

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

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

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

   
7.3.2 Header interpretation example 2

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

 

 

Table 6: Sample calculations for example 1.

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

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

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

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

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

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

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

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

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

 

 

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

123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789

NAXIS = 2 / 2-dimensional image

NAXIS1 = 2048 / x axis (fastest)

NAXIS2 = 2048 / y axis (2nd fastest)

MJD-OBS = 44258.7845612 / MJD at start of observation

CRPIX1 = 1024.5 / Pixel coordinate of reference point

CRPIX2 = -1023.5 / Pixel coordinate of reference point

PC1_1 = 1 / Transformation matrix element

PC1_2 = -0.004 / Transformation matrix element

PC2_1 = -0.002 / Transformation matrix element

PC2_2 = 1 / Transformation matrix element

CDELT1 = -0.005 / x-scale

CDELT2 = 0.005 / y-scale

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

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

PV2_1 = -25.0 / Conic mid-latitude

CRVAL1 = 90.0 / Galactic longitude at reference point

CRVAL2 = -25.0 / Galactic latitude at reference point

CRPIX1A = 1024.5 / Pixel coordinate of reference point

CRPIX2A = -1023.5 / Pixel coordinate of reference point

PC1_1A = 1 / Transformation matrix element

PC1_2A = -0.004 / Transformation matrix element

PC2_1A = -0.002 / Transformation matrix element

PC2_2A = 1 / Transformation matrix element

CDELT1A = -0.005 / x-scale

CDELT2A = 0.005 / y-scale

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

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

PV2_1A = -25.0 / Conic mid-latitude

CRVAL1A = -7.0300934 / Ecliptic longitude at reference point

CRVAL2A = 34.8474143 / Ecliptic latitude at reference point

LONPOLEA= 6.3839706 / Native longitude of ecliptic pole

LATPOLEA= 29.8114400 / Ecliptic latitude of native pole

RADESYSA= 'FK5 ' / Mean IAU 1984 ecliptic coordinates

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

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

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

7.3.3 Binary table representations of example 2

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

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

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

 

 

Table 8: Calculations for example 2.

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

   
7.3.4 Header interpretation example 3

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

 

 

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

123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789

XTENSION= 'BINTABLE' / Binary table extension

BITPIX = 8 / 8-bit bytes

NAXIS = 2 / 2-dimensional binary table

NAXIS1 = 16777232 / Width of table in bytes

NAXIS2 = 4 / Number of rows in table

PCOUNT = 0 / Size of special data area

GCOUNT = 1 / One data group (required keyword)

TFIELDS = 5 / number of fields in each row

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

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

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

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

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

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

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

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

TTYPE5 = 'Image ' / 2-D image array

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

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

MJDOB5 = 44258.7845612 / MJD at start of observation

COMMENT The following keywords define the primary coordinate description

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

11PC5 = 1 / Transformation matrix element

12PC5 = -0.004 / Transformation matrix element

21PC5 = -0.002 / Transformation matrix element

22PC5 = 1 / Transformation matrix element

1CDE5 = -0.005 / Axis 1: scale

2CDE5 = 0.005 / Axis 2: scale

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

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

2PV5_1 = -25.0 / Conic mid-latitude

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

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

COMMENT The following keywords define the secondary coordinate description

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

11PC5A = 1 / Transformation matrix element

12PC5A = -0.004 / Transformation matrix element

21PC5A = -0.002 / Transformation matrix element

22PC5A = 1 / Transformation matrix element

1CDE5A = -0.005 / Axis 1A: scale

2CDE5A = 0.005 / Axis 2A: scale

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

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

2PV5_1A = -25.0 / Conic mid-latitude

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

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

LONP5A = 6.3839706 / Native longitude of ecliptic pole

LATP5A = 29.8114400 / Ecliptic latitude of native pole

RADE5A = 'FK5 ' / Mean IAU 1984 ecliptic coordinates

EQUI5A = 2000.0 / Coordinate epoch

END

 

 

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

123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789

XTENSION= 'BINTABLE' / Binary table extension

BITPIX = 8 / 8-bit bytes

NAXIS = 2 / 2-dimensional binary table

NAXIS1 = 5 / Width of table in bytes

NAXIS2 = 10000 / Number of rows in table

PCOUNT = 0 / Size of special data area

GCOUNT = 1 / One data group (required keyword)

TFIELDS = 3 / Number of fields in each row

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

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

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

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

TLMIN2 = 1 / Lower limit of axis in column 2

TLMAX2 = 2048 / Upper limit of axis in column 2

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

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

TLMIN3 = 1 / Lower limit of axis in column 3

TLMAX3 = 2048 / Upper limit of axis in column 3

MJDOB3 = 44258.7845612 / MJD at start of observation

COMMENT The following keywords define the primary coordinate description.

TCRP2 = 1024.5 / Axis 1: reference pixel coordinate

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

TPC2_2 = 1 / Transformation matrix element

TPC2_3 = -0.004 / Transformation matrix element

TPC3_2 = -0.002 / Transformation matrix element

TPC3_3 = 1 / Transformation matrix element

TCDE2 = -0.005 / Axis 1: scale

TCDE3 = 0.005 / Axis 2: scale

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

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

TPV3_1 = -25.0 / Conic mid-latitude

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

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

COMMENT The following keywords define the secondary coordinate description.

TCRP2A = 1024.5 / Axis 1A: reference pixel coordinate

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

TP2_2A = 1 / Transformation matrix element

TP2_3A = -0.004 / Transformation matrix element

TP3_2A = -0.002 / Transformation matrix element

TP3_3A = 1 / Transformation matrix element

TCDE2A = -0.005 / Axis 1A: scale

TCDE3A = 0.005 / Axis 2A: scale

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

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

TV3_1A = -25.0 / Conic mid-latitude

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

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

LONP3A = 6.3839706 / Native longitude of ecliptic pole

LATP3A = 29.8114400 / Ecliptic latitude of native pole

RADE3A = 'FK5 ' / Mean IAU 1984 ecliptic coordinates

EQUI3A = 2000.0 / Coordinate epoch

END

 

 

Table 11: Third example FITS header.

123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789

NAXIS = 2 / 2-dimensional image

NAXIS1 = 181 / x axis (fastest)

NAXIS2 = 91 / y axis (2nd fastest)

CRPIX1 = 226.0 / Pixel coordinate of reference point

CRPIX2 = 46.0 / Pixel coordinate of reference point

CDELT1 = -1.0 / x-scale

CDELT2 = 1.0 / y-scale

CTYPE1 = 'GLON-CAR' / Plate carree projection

CTYPE2 = 'GLAT-CAR' / Plate carree projection

CRVAL1 = 30.0 / Galactic longitude at reference point

CRVAL2 = 35.0 / Galactic latitude at reference point

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

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

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

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

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

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

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

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

   
7.4 Header construction examples

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

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

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

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

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

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

   
7.4.1 Header construction example 1

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

   
7.4.2 Header construction example 2

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

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

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

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

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

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

By inspection of the two sets of equations we must have

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

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

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

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

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

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

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

whence

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

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

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

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

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

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

   
7.4.3 Header construction example 3

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

 

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

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

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

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

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

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

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

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

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

1.

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

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

2.

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

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

3.

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

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

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

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

   
7.4.4 Header construction example 4

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

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

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

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

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

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

 

 

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

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

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

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

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

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

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

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

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

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

 

 

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

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

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

   
7.5 Realization

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

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