[fitswcs] Status of WCS negotiations
Patrick Wallace
ptw at star.rl.ac.uk
Tue Jul 28 12:36:33 EDT 1998
Warning - this message is long and boring.
BACKGROUND
----------
On July 20th, Don Wells mentioned my having drawn attention to an
apparent change in the meaning of the established CDELTn keyword. I've
known Don long enough to be confident that "the water has been muddied"
was not meant to sound pejorative. But Mark Calabretta's response
(July 21st), "Technically invalid arguments shouldn't be accepted by
scientists" is harder to ignore!
Accordingly, for anyone who's interested, I'll try to outline the history
of this issue, then I'll give a technical example where I feel the G+C
proposals lead us into trouble, and finally suggest some ways forward.
BASIC ASSUMPTIONS
-----------------
To keep things really simple, I'll deal only with sky coordinates and
stick to easy cases where, although there is scaling and rotation, there
are no zero offsets and no skew.
Coordinates (i,j) are pixels. Coordinates (x,y) are sky coordinates,
for example "standard coordinates" (xi,eta) for some celestial reference
frame.
I'm going to ignore the historical problems with the signs of various
angles and matrix elements. If it turns out that my definition of a
CROTA angle or CD element has the wrong sign, it shouldn't matter.
HISTORY
-------
There are two existing conventions for dealing with world coordinates in FITS:
1) The AIPS conventions have scaling factors CDELT1 and CDELT2 which you
apply to pixel coordinates (i,j) prior to rotation through the angle
CROTA2.
x = CDELT1*i*cos(CROTA2) - CDELT2*j*sin(CROTA2)
y = CDELT1*i*sin(CROTA2) + CDELT2*j*cos(CROTA2)
(Remember I've left out the zero offsets CRPIXn, to avoid clutter.)
2) The Hanisch & Wells conventions employ a matrix, the elements of which
describe both the scaling and orientation (and any skew):
x = CD001001*i + CD001002*j
y = CD002001*i + CD002002*j
Drafts of Greisen & Calabretta c1993 used the CD matrix, albeit with a
different subscript naming convention. No CDELTs were required, because
the scaling was defined by the CD matrix. Sometime between 1993 and the
1996 draft, the CD matrix was dropped. The PC matrix had appeared in its
place, and the CDELTs had returned...except the order in which they were
applied had changed, giving this result:
not 2 as in AIPS
v
x = CDELT1*i*PC001001 + CDELT1*j*PC001002
y = CDELT2*i*PC002001 + CDELT2*j*PC002002
^
not 1 as in AIPS
This change was first brought to my attention by my colleague David Berry,
in March 1996. He pointed out that, as a consequence of the change in
meaning, G+C's formulae for interpreting old, AIPS-style FITS data in
the new scheme were wrong. Told of this by various people, G+C responded
by changing the formulae (Eqn 130 in the old drafts, Eqn 132 in recent
drafts), introducing adjustments to the PC matrix to undo the effects of
applying the CDELTs after rather than before the "rotation".
I waded into the debate about a year later. I had written a library of WCS
functions for the Gemini project, one of which generated FITS headers for
Gemini imagers. I had, naturally, embraced the latest and best FITS
conventions, namely the G+C proposals, and I did succeed in generating FITS
headers which expressed the astrometry properly and which used the G+C
conventions. However, there were certain awkwardnesses about the process,
and I began to suspect that the G+C proposals, though perfectly workable in
a technical sense, would in practice catch people out. This matters; I
fully concur with G+C when they say (p3) "Good form in FITS always extends
to include the human as well as the software readers".
My suspicions were strengthened when, as a check that I'd implemented the
new design properly, I gave Bill Pence one of my sample headers and asked
him to tell me the (RA,Dec) of the corner pixels. Because his software
supported the existing de facto methods rather than the new G+C proposals,
he needed a CROTA2 angle, and so he used the obvious interpretation of the
PC matrix to derive one. He then used this CROTA2 estimate in conjunction
with the CDELT1 and CDELT2 values that I had supplied, and duly got wrong
answers.
So I took it up with G+C...and was unable to persuade them that there was
any problem at all. The essence of the disagreement was that whereas I
maintained that the old and new CDELTs mean something different, Mark
adamantly denied this (see the appendix to his July 21st message). What
he has claimed all along is that the new CDELTs are simply more general
than the old ones, incorporating existing uses of CDELTs while adding
new capabilities. He points out the following:
(i) When CDELT1=CDELT2 (the case I refer to as "square pixels"),
it doesn't matter which order you do things in and hence the
old and new CDELT values are the same.
(ii) When CROTA2=0, you can use the identity matrix for PC so that
it doesn't matter which order you do things in and hence the
old and new CDELT values are the same.
(iii) When there isn't a PC matrix you can contrive one that makes
the CDELTs have the same numerical values in the old and the
new schemes.
However, (i) and (ii) are special cases and not representative of the
sort of real examples I had in mind. Point (iii) is valid, though
for my money a bit kludgey. But the important case that Mark has not
dealt with is where (a) the pixels aren't square, (b) the rotation isn't
zero, and (c) you want to write a new-style FITS file *that can be
correctly interpreted by existing readers*. In the next section, I'll
give an example of this sort and demonstrate that the G+C precepts lead
to incorrect results for an image generated in an obvious and legal way.
But first, at the risk of labouring the point, let me prove that CDELTn
has changed in meaning by writing down the two definitions.
AIPS convention:
CDELT1 is the distance on the sky between the centres of two
pixels whose detector coordinates differ by (1,0). CDELT2
is the same for two pixels whose detector coordinates differ
by (0,1).
G+C proposal:
CDELT1 is the distance on the sky between the centres of two
points whose "true plane pixel coordinates" differ by (1,0).
CDELT2 is the same for two points whose a "true plane pixel
coordinates" differ by (0,1).
Obviously, "detector coordinates" and "true plane pixel coordinates"
are not, in general, the same thing. Hence the two definitions *are*
different, which I feel on its own is enough to justify using a different
nomenclature. Yes, you can devise a PC matrix that is compatible with
the CROTA and CDELTn of an existing AIPS-style FITS image and preserves
the old CDELTn numbers, so for this application the difference in
meaning is academic. But what about the reverse - when you want older
FITS-reading software (or a human reader) to interpret a new-style FITS
image?
AN EXAMPLE
----------
Our detector has 2:1 oblong pixels, as a result of binning. The pixel
size, projected onto the sky, is:
0.0001 deg in i
0.0002 deg in j
The camera can rotate to any angle, CROTA. The relationship between
pixel coordinates (i,j) and sky coordinates (x,y) is:
x = +0.0001*i*cos(CROTA)-0.0002*j*sin(CROTA)
y = +0.0001*i*sin(CROTA)+0.0002*j*cos(CROTA)
On this particular occasion the angle is 90 deg. Therefore:
x = -0.0002*j
y = +0.0001*i ....(1)
Our FITS-writing software, we will suppose, uses CROTA directly when it
generates the headers. (It does not treat CROTA = 90, 180, 270 as
special cases - why should it have to? - so there's going to be no
swapping around of axes or anything like that. In any case, the only
reason I picked a cardinal point for this example was to make the mental
arithmetic easier.) The obvious way to generate the PC matrix (and I
assume our FITS-writing program is entitled to do this, because G+C
doesn't say otherwise) is as follows:
| PC001001 PC001002 | | +cos(CROTA) -sin(CROTA) | | 0 -1 |
| | = | | = | |
| PC002001 PC002002 | | +sin(CROTA) +cos(CROTA) | | +1 0 |
To deliver the required transformation, equations (1) tell us that we
require a CDELT matrix as follows:
| CDELT1 0 | | +0.0002 0 |
| | = |
| 0 CDELT2 | | 0 +0.0001 |
So far, so good. We have used the G+C precepts to generate this set
of keywords, which validly describe the astrometry for our case:
PC001001 0
PC001002 -1
PC002001 +1
PC002002 0
CDELT1 +0.0002
CDELT2 +0.0001
Now we turn to Section 5 of G+C, where we are asked to help old FITS
readers by writing a CROTA2 value (I know from private correspondence
that Mark is against this, but that's what the present draft recommends
we do). Equation 133 shows us how to do this:
CROTA2 = (arg(0,1)+arg(0,1))/2 = 90 deg
So our old FITS reader, which doesn't know anything about PCiiijjj,
sees these keywords:
CROTA2 = 90
CDELT1 = +0.0002
CDELT2 = +0.0001
When it transforms pixel (i,j) to sky (x,y) it ends up using:
x = -0.0001*j
y = +0.0002*i ....(2)
Equations (1) and (2) disagree. QED.
WAYS FORWARD
------------
PLAN A:
* Drop PC/nuCDELT and revert to a CD matrix approach.
* Provide for multiple CD matrices, eliminating the need for separate
scaling factors and supporting alternative astrometric solutions
(e.g. observed place for definitive accuracy, ICRS for practical
applications).
* Recommend including CROTA2, CDELT1 and CDELT2 in headers, generated
from the CD matrix, for backwards compatibility.
PLAN B (to retain post-rotation scaling and hence secondary axis
descriptions):
* Use a different name for the scale factors, not CDELTn.
* In Section 5, recommend helping existing readers by expressing PC +
nuCDELTn as CDELTn + CROTA (in low-skew cases), and using a PC matrix
that requires unit nuCDELTn (in other words the CD matrix).
PLAN C (minimal changes to fix up G+C):
. Include a definition of the term "TRUE PLANE PIXEL COORDS" and spell
out the old and new meanings of CDELTn.
. Either (i) omit Section 5 altogether, or (ii) lay down restrictions
on how the PC matrix should be constructed so that the rules of
Section 5 work properly in the case where we want old AIPS-compatible
readers to understand new FITS files, or (iii) say that we aren't
going to address this case, and eliminate the recommendation to write
CROTA2.
Finally, I should stress that I am prepared to vote for whatever consensus
emerges in this and all other aspects of the WCS discussion, and support
the Chairman's efforts to reach such a consensus.
Patrick Wallace
____________________________________________________________________________
Starlink Project Manager Internet: ptw at star.rl.ac.uk
Rutherford Appleton Laboratory Tel: +44-1235-445372
Chilton, Didcot, Oxon OX11 0QX, UK Fax: +44-1235-446667
____________________________________________________________________________
More information about the fitswcs
mailing list