Comments on Draft WCS Standard
Fred Patt (301)286-4569
rather than to "fits...@nrao.edu", and so it came to me, Don Wells,
the maintainer of the exploder, rather than being exploded as
intended.
------- Start of forwarded message -------
From: fr...@manono.span (Fred Patt (301)286-4569)
To: fitsbits...@fits.CV.NRAO.EDU
Subject: Comments on Draft WCS Standard
Date: Thu, 18 Feb 93 11:12:32 EST
Dr. Greisen,
I was offered a chance by Richard White to review your draft World
Coordinate System standard for FITS, due to my interest in map projections
in general and the quad sphere in particular. My recent involvement has
been with Earth science data sets but I was also involved with the COBE
mission for several years and am familiar with their skymap implementation.
I have advocated the adoption of the quad sphere by local Earth science
projects and was responsible for the use of this scheme in the Advanced
Very High Resolution Radiometer (AVHRR) Pathfinder Land prototype
processing system here at Goddard. I provided details of the projection
equations to the AVHRR group in the form of pseudocode, which has been
implemented in C by the developers. The data system manager, Mary James of
Goddard Code 902, has given me encouragement to publish my work on this
project.
My point here is that an exact form of the projection was derived
in the second of the two references (O'Neill and Laubscher, 1976) in terms
of polar coordinates, and this is the form which was implemented by AVHRR.
The reference only derived the projection in the forward direction, from
sky to cube coordinates, but I found that the derivation of the reverse
projection was very straightforward. I also found that the equations could
be significantly simplified by direct use of cartesian coordinates, thereby
eliminating the polar angles as intermediaries and, as a side benefit,
reducing the number of trig function evaluations required.
The information regarding the COBE implementation by Immanuel
Freedman is completely accurate. COBE has accepted the non-exact form of
the transformation, but I suggest that the description of the exact form is
at least as appropriate for your work. The combination of the overall
complexity of the projection equations and the errors in the projection
described for COBE would seem to discourage most potential users.
I have attached a set of summary notes which I provided for the
AVHRR group and the derivation of the reverse projection and cartesian form
equations for the exact transformation (in computer-language syntax; it
would be reasonable for me to put these into a more appropriate notation if
desired). I hope that this information may be of use to you.
Fred Patt
SUMMARY
The quadrilaterlized spherical cube, or quad sphere for short, is
an equal-area mapping and binning scheme for data collected on a spherical
surface (either Earth data or the celestial sphere). It was first proposed
by Chan and O'Neill in 1975 for the Naval Environmental Prediction Research
Facility (Reference 1).
There are two key elements to the quad sphere:
o The mapping consists of projecting the sphere onto the faces of an
inscribed cube using a curvilinear projection which preserves area. The
sphere is divided into six equal areas which correspond to the faces of the
cube. The vertices of the cube correspond to the cartesian coordinates
defined by |x|=|y|=|z| on the unit sphere. For an Earth projection, the
cube is normally oriented with one face normal to the North Pole and one
face centered on the Greenwich meridian (although any definition of pole
and meridian could be used). The faces of the cube are divided into square
bins, where the number of bins along each edge is a power of 2, selected to
produce the desired bin size. Thus the number of bins on each face is
2**(2*N), where N is the binning level, and the total number of bins is
6*2**(2*N). For example, a level of 10 gives 1024x1024 bins on each face
and 6291456 (6*2**20) total bins, which are 23.605 square arcminutes
(1.99737E-6 steradians) in size.
o The bins are numbered serially, rather than being rastered as for an
image. The bin numbers are determined as follows. The total number of
bits required for the bin numbers at level N is 2*N+3, where the 3 MSBs are
used for the face numbers and the remaining bits are used to number the
bins within each face. The faces are numbered 0-5 with 0 being the North
face, 1 through 4 being equatorial with 1 corresponding to Greenwich, and 5
being South. Thus at level 10, face 0 has bin numbers 0-1048577, face 1
has numbers 1048576-2097151, etc. Within each face the bins are numbered
serially from one corner (the convention is to start at the "lower left")
to the opposite corner, with the ordering such that each pair of bits
corresponds to a level of bin resolution. This ordering in effect is a
two-dimensional binary tree, which is referred to as the quad-tree The
conversion between bin numbers and coordinates is straightforward. If
4-byte integers are used for the bin numbers the maximum practical bin
level is 14, which uses 31 of the 32 bits and results in a bin size of
0.0922 square arcminutes (7.80223E-9 steradians). The advantages of this
numbering approach are stated below.
In principal the mapping and numbering schemes are separable; the
projection onto the cube could be used with another bin numbering scheme,
and the numbering scheme itself could be used with any arrangement of bins
which can be partitioned as a set of square arrays. Used together, they
comprise a flexible and efficient system for storing map data.
ADVANTAGES
The quad sphere projection does not produce singularities at the
poles or elsewhere, as do some other equal-area mapping schemes. The
distortion is moderate over the entire sphere, so that at no point are
shapes distorted beyond recognition.
Individual faces can be used to display parts of the sphere without
remapping by converting from pixel order to raster order; coordinate grids
can be overlaid on the faces to assist in registration.
The bin numbering allows the resolution to be changed (by factors
of two) simply by adding or deleting LSBs. This is particularly useful if
it is desired to increase the bin size, for example to compare maps with
different resolutions. This is performed simply by dividing the bin
numbers by four for each factor of two increase in bin size.
The bin numbers tend be closer together in sequence for neighboring
bins.
The quadtree bin numbering allows maps to be stored as FITS
BINTABLES with no wasted storage needed for blank areas (unfilled bins).
Only bins containing valid data need be included in a file, as the bin
numbers themselves identify the locations of the bins.
The serial numbering frees up one array dimension for another use
(e.g., time, a third coordinate, multichannel data, etc.).
DERIVATIONS OF REVERSE AND CARTESIAN TRANSFORMATIONS
The equations for the polar form of the transformation from the sphere to
the cube are given in reference 2 as equations (3-21) and (3-38). I have
taken this approach further in two areas: 1) inversion of these equations
to derive the cube-to-sphere transformation; and 2) simplification of these
equations to perform the transformation directly in terms of cartesian
coordinates. Both of these are straightforward as shown below.
1. DERIVATION OF THE CUBE-TO-SPHERE TRANSFORMATION
The exact form of the sphere-to-cube tranformation in reference 2 is
derived in terms of the co-elevation and azimuth angles (theta, phi) on the
sphere and an equivalent set of angles (mu,nu) on the cube. The equal-area
tranformation is performed using the following equations:
TAN(mu) = (12/pi)*(theta+ACOS(SIN(theta)*SIN(pi/4))-pi/2) (3-21)
TAN(nu) = SQRT((SEC(mu))**2*(1-COS(theta))/(1-COS(ATAN(SEC(theta))))) (3-38)
The inversion of these equations is as follows. Using the trig identity
ASIN(x)+ACOS(x)=pi/2, and the knowledge that SIN(pi/4) = 1/SQRT(2),
equation (3-21) can be written as
TAN(mu) = (12/pi)*(theta-ASIN(SIN(theta)/SQRT(2))) (3-21a)
and then rearranged as
theta - (pi/12)*TAN(mu) = ASIN(SIN(theta)/SQRT(2))
Taking the sine of both sides gives
SIN(theta)*COS((pi/12)*TAN(mu)) - COS(theta)*SIN((pi/12)*TAN(mu))
= SIN(theta)/SQRT(2)
Dividing through by COS(theta) and solving for TAN(theta) yields the
desired inverse relationship:
TAN(theta) = SIN((pi/12)*TAN(mu))/(COS((pi/12)*TAN(mu))-1/SQRT(2)) (1)
The inversion of equation (3-38) is performed by squaring both sides and
solving for COS(phi):
COS(phi) = 1 - (TAN(nu))**2*(1-COS(ATAN(SEC(theta))))/(SEC(mu)**2 (2)
As in the original derivation, this equation used the prior determination
of theta as a function of mu.
2. EXPRESSION OF THE TRANSFORMATIONS IN CARTESIAN COORDINATES
Both sets of transformation equations can be expressed in cartesian
coordinates (three-dimensional coordinates of the unit vector on the sphere
and two-dimensional coordinates on the cube face). Assume the unit vector
has coordinates (q,r,s) and the coordinates on the cube face are (x,y),
where q is normal to the plane of the cube face and (r,s) cooresponds to
(x,y). Then the functions of angles can be expressed as follows:
TAN(theta) = s/r
SIN(theta) = s/SQRT(r*r+s*s)
COS(phi) = q
TAN(mu) = y/x
TAN(nu) = SQRT(x*x+y*y)
Substituting these into the transformation equations results in the
following revised forms for the equations.
Equation (3-21a):
y/x = (12/pi)*(ATAN(s/r)- ASIN(s/SQRT(2*r*r+2*s*s))) (3)
Equation (3-38) can first be simplified by the following substitution:
COS(ATAN(SEC(theta))) = 1/SQRT(1+(SEC(theta))**2)
Substituting this and the appropriate cartesian coordinate expressions
gives:
x*x+y*y = (1+(y/x)**2)*(1-q)/(1-1/SQRT(2+(s/r)**2))
which can be re-written to eliminate y:
x = sqrt((1-q)/(1-1/SQRT(2+(s/r)**2))) (4)
Note that there is no ambiquity in the sign of the square root since the
derivation of (3-38) assumes that x is positive. This and equation (3) for
y/x above give the solution for both y and x, noting also that (3-21)
assumes that x is not less than y.
Equation (1) above can be rewritten also by direct substitution:
s/r = SIN((pi/12)*(y/x))/(COS((pi/12)*(y/x))-1/SQRT(2)) (5)
Equation (2) above can be rewritten by solving (4) for q:
q = 1 - x*x*(1-1/SQRT(2+(s/r)**2))) (6)
Thus equations (3) through (6) can be used to convert directly from
coordinates on the sphere to the cube and back without the need to convert
to angles as intermediate variables, and the number of trigonometric
function evaluations required is greatly reduced.
REFERENCES
1. Chan, F.K. and O'Neill (1975), Feasibility Study of a Quadrilateralized
Spherical Cube Earth Data Base, Computer Sciences Corporation, EPRF
Technical Report 2-75 (CSC). Prepared for the Environmental Prediction
Research Facility, Monterey, California.
2. O'Neill, I.M. and Laubscher, R.E. (1976), Extended Studies of a
Quadrilateralized Spherical Cube Earth Data Base, Computer Sciences
Corporation, NEPRF Technical Report 3-76 (CSC). Prepared for the Naval
Environmental Prediction Research Facility, Monterey, California.
------- End of forwarded message -------