FITS WCS Discussions and Email

NOAO conducted an internal review of the FITS WCS proposal over the period December 2000 through June 2001. The purpose of this review was to do in-depth studies of the applicability of FITS WCS to certain areas of high priority to NOAO, and then provide feedback to the community regarding possible further revision of the WCS proposal. The areas of application studied included development of a WCS formulation for spectroscopic data, application of FITS WCS to astrometric calibration of wide field cameras such as the NOAO Mosaic, and integration of FITS WCS into data analysis systems such as IRAF.

This section provides a record of the various emails and reports generated during the FITS WCS review at NOAO through June 2001 (as it turned out Don Wells also got heavily involved and his email is included as well). We are making this record available to provide additional in-depth material documenting the thinking leading up to our proposed revisions. Some of the text has been edited slightly for inclusion here.

Date: Mon, 11 Dec 2000 15:13:17 -0700 (MST)
From: Doug Tody <tody@tucana.tuc.noao.edu>
To: Don Wells <dwells@nrao.edu>
Subject: WCS issues (fwd)

Don -  Lindsey says you called to ask about WCS issues.  Our working plan
here regarding WCS has been to proceed as follows:

    o   Revise FITS MEF format for Mosaic CCD data so that it is fully
	self contained and self descriptive (this involves some WCS issues
	and related layered conventions).  Finish implementing new Ximtool
	and related MSCRED changes (the new version of Ximtool fully
	supports real time WCS display, including multiple images each
	with multiple WCS).

    o	Finish implementation in IRAF of image compression, pixel mask,
	and variance plane extensions (a partial implementation already
	exists).

    o	Design new IRAF (FITS-based) spectral data format.  Initially this
	will be just the data format, with the resolution of the spectral
	WCS issues to follow second.

	<interaction with outside groups on WCS [spectral data] format
	(STScI, UK, Lick)>

    o	With the finer points of the above in mind, do a careful review
	of the latest FITS WCS drafts and propose any changes needed to
	address general image WCS concerns (not just the sky projections).

	<interaction with standards committees re general WCS)

    o	Consider the spectral WCS issues and propose any changes.  Test
	implementation of spectral WCS in IRAF.

	<interaction with standards committees re spectral WCS)

	We are also working on a related convention for describing
	multiaperture spectral data using FITS keywords.  I don't think
	this needs to be part of the basic spectral WCS representation,
	but it does use WCS and is related.

[...] there will be very little time for this until at least mid-January.

	- Doug


Date: Mon, 11 Dec 2000 17:41:11 -0500
From: Don Wells <dwells@cv.nrao.edu>
To: Doug Tody <tody@noao.edu>
Subject: Re: WCS issues (fwd)

Doug,

At the moment my main concern and highest priority is the question of
whether the negotiations for Paper-I(General) and Paper-II(Celestial)
are completed.  If these two drafts are acceptable to everybody I want
to begin moving toward putting them to a vote in the Committees.

Do you expect that, as a result of the review you are planning, you
will propose any change to the content of Paper-I or Paper-II?

Doug Tody writes:
 > o Design new IRAF (FITS-based) spectral data format. Initially this
 >   will be just the data format, with the resolution of the spectral
 >   WCS issues to follow second.
 >   <interaction with outside groups on [spectral] format (STScI, UK, Lick)>

It appears to me that Paper-III(Spectroscopic) needs to add support
for distorted spectra. Bill Cotton and Walter Jaffe (Leiden) are
probably going to propose a polynomial correction function; they are
motivated by the case of optical interferometers, in which the fringes
are spectrally resolved.  Your interest in a new spectral design is an
opportunity for the whole community to finally address this issue that
we have put off for twelve years. 

I think that it would be entirely reasonable for the General(I) and
Celestial(II) papers to proceed to a vote, and for the Spectral(III)
paper to continue in negotiation for months more.

[...]

-Don


Date: Mon, 11 Dec 2000 15:53:52 -0700 (MST)
From: Doug Tody <tody@tucana.tuc.noao.edu>
To: Don Wells <dwells@cv.nrao.edu>
Subject: Re: WCS issues (fwd)

> At the moment my main concern and highest priority is the question of
> whether the negotiations for Paper-I(General) and Paper-II(Celestial)
> are completed.  If these two drafts are acceptable to everybody I want
> to begin moving toward putting them to a vote in the Committees.
> 
> Do you expect that, as a result of the review you are planning, you
> will propose any change to the content of Paper-I or Paper-II?

Don - As we discussed at the ADASS, we think we may have several minor
issues yet to resolve with Paper-I and II (referred to as "general WCS"
in my previous email).  The best time (for us at least) to address this
would be in January.  The ximtool/MSCRED WCS-related work currently in
progress here, for example, is a complex WCS application that is not just
a classical celestial coordinate system WCS application.

I hope you weren't planning to try to stir things up just before the
holidays?  That would likely get a very mixed response.  What with the
holidays and the AAS, the schedule is pretty well packed already through
mid-January.

I think papers I and II are very close, but we should allow time for one
last careful review before proceeding.  If we do this and everyone is in
good agreement informally, I would guess that the subsequent voting review
process could be sped up greatly.

	- Doug


Date: Mon, 2 Apr 2001 10:13:48 -0700 (MST)
From: Lindsey Davis <davis@noao.edu>
To: tody@perseus.tuc.noao.edu
Subject: FITS WCS


                Comments on the FITS WCS Papers
		        To Doug Tody
		     From Lindsey Davis
		        April 2001

	     [Revised (minor edits only) June 2001]

1. Paper I

1.1. Default Values of the CD matrix (FITS convention / IRAF convention)

    Section 2.1 and Table 1 define the default values of the CDj_is matrix
to be 0.0 if i != j and 1.0 if i = j. This means that if no CD matrix
elements are defined in the image header then the CD matrix defaults to 
the identity matrix, which is the same as the convention adopted by IRAF.
If only the diagonal elements of the CD matrix are defined then the
off-diagonal elements default to 0.0, which is also the convention adopted
by IRAF. However if the off-diagonal elements are defined then the diagonal
elements default to 1.0 which is not the convention adopted by IRAF.

    Basically IRAF assumes that if the CD matrix is undefined it defaults
to the identity matrix but, if any single element of the CD matrix is defined
all the undefined elements default to 0.0 even the diagonal ones. This is
economical, but is in contradiction to the proposed FITS standard. This
produces a portability problem for us. See the example below.

If the IRAF header contains the following  CD matrix elements

    CD1_2 = 2.0E-05
    CD2_1 = 2.0E-05

then IRAF assumes

    CD1_1 = 0.0
    CD2_2 = 0.0

whereas the proposed FITS standard assumes

    CD1_1 = 1.0
    CD2_2 = 1.0

Note that under the new standard existing IRAF headers may be wrong. In many
cases precision in the fitting process e.g. the fact that the CD elements are
not quite zero due to roundoff may save us from the effects of the proposed
FITS standard, but in some cases where the user forces a perfectly aligned
WCS there may be a problem.

    I did not realize the implications of this until quite recently
when a user claimed to have trouble reading data with AIPS. I thought is
was the usual CDELT/CROTA2 versus CD matrix issue, but it may actually have
been the above problem. 

    The FITS people will reply that while writing transposed headers as above
is legal, see section 2.1, it is strongly discouraged, and in that case
we should have explictly set the diagonal elements to 0.0. However we have
been using the CD matrix notation for quite a while with our current
convention. I took a look at the original CD matrix paper by Hanisch and
Wells. There is brief mention that missing off-diagonal terms should assume
default values of 0.0 but no formal mention of what missing diagonal terms
should default to.

    Basically the problem boils down to different conventions and we need
to decide which is better or at least how this affects the data archives.
The new convention causes portability problems for us.


1.2. Default values for CRPIX and CRVAL (FITS convention ?)

    As far as I can tell default values for CRPIX and CRVAL are not defined
in section 2.1 and are not listed in Table 1. I assume that this is just
an oversight. IRAF has adopted the convention that these quantities
default to  0.0 if undefined, and I assume that the unstated FITS convention
is the same ?


1.3 String Parameter Values (Better Interface / Compatibility with IRAF)

    There is no support for string valued functon parameters. The PVj_ms are
by definition floating point numbers. To date in IRAF we have used builtin
string parameters for coordinate type (axtype), projection type (wtype),
units(units), label(label), and format (format) quantities, and function 
specific string parameters for encoding a variety of arbitrary functions
including polynomials, splines, etc. The coordinate (axtype) and projection
types (wtype) have always been aliased to the FITS CTYPEjs parameter, and the
units quantity can know be encoded in the proposed FITS CUNITjs parameter.
However there is no analog to the label and formats quantities. People will
argue that the label quantity is not that important (the important info is
in the units parameter and I agree with this), and that encoding formats
is tricky anyway as they are system and language dependent (and I also agree
with this comment).

    However I think the standard would be significantly enhanced if
string parameters in which any quantity, characters, integers, reals,
keyword=value pairs could be encoded were supported. What the quantities
mean would be determined solely by the function type defined by the
CTYPEjs parameter, which need not be one of the standard functions but
may instead be a personal or user function. We already do this in IRAF
and have found it to be a very powerful and flexible way to repreresent
complicated and arbitrary functions, including splines, sampled data,
arbitrary not fixed form polynomials, etc. In my opininin it is also a
more compact and readable way of transmitting parameter information. One use
I can see for it is to carry along information about the the function actually
used to compute the plate solution, dispersion function etc, which may be
different the than the form allowed by one of the standard functions.

    One possible way to implement this that might be more acceptable to the
current FITS community than the IRAF scheme, is to adopt a set of string
parameters WVj_ns where j is the usual axis number (1-99) and n is a sequence
number which is (1-99) or (1-999) if we allow the sequence number to grow, if
the number of axes < 10 which is the usual case. The actual value of the
string parameter would be the concatenation of all the SVj_ns cards which
have the same value of j and s. Using this scheme we immediately remove the
builtin limitation of 99 parameters and acquire instead a limitation of 99
FITS cards, or 999 if we don't get too restrictive on j.

    There are some currently supported IRAF spectral WCS's that cannot be fit
into the current FITS WCS proposal scheme but could be supported legally
in FITS without change if string parameters were supported. If adopted this
scheme woul means that the current TNX and ZPX IRAF WCSs could be considered
as engineering data and accomodated without change. This does not mean that
these WCSs' should be adopted as standard FITS WCS functions. We will adopt
one of the standards for export in the future. However we will still need to
support these systems and it would be nice to do so inside the FITS WCS
framework.

    Adopting an WVj_n parameter does not affect the standard or standard
functions in any way but it does provide more flexibilty to FITS representation.


1.4 Dimensional Reductions and the WCS (FITS conventions / IRAF conventions)

    One issue that is not discussed explictly in the first paper is the
issue of dimensional reduction, i.e what happens to the WCS if you extract
a line from a 2D image, a 2D plane from a 3D cube etc? This is something we 
do all the time in IRAF.

    My impression is that in FITS you either just lose the information on
the reduced axis, which is ok if the axis coordinate transformation are
independent of one another and you don't care about the image history, or
that dimensional reduction is not allowed, i.e. you can only extract a 2D
line with the length of the second axis set to 1 from a 2D image. Is this
true ? Does this mean one a can never represent a true 1D image in FITS by
extracting a line from an image with a celestial coordinate WCS where the
coordinate axes are dependent and information for both coordinates is
required ? or is this just a convention ?

    I would like to propose an alternative convention that could be equally
acceptable and that is to permit the dimension of the WCS to be higher than
that of the image encoded in NDIM. The FITS interpreter would look
to see if there are any CTYPE/CRPIX/CRVAL/CD/ etc keywords with j
indices greater than NDIM and interpret the WCS accordingly. I see no reason
why this convention can not be as valid a convention as the current one
of maintaining image axes of length 1.
    
    In IRAF the WCS may have indeed have a dimension greater than the image
if it was derived from a parent image. We use the WCSDIM keyword which
defaults to the dimension of the image if absent. There is an accompanying
axis map keyword WAXMAP which defines how the WCS axes of the parent
image map to the current image axes. This parameter is a string parameter
whose default value is "" for a one to one mapping. 

    I think there should be some discussion of this issue in the first paper.
I personally think the current convention of using extra axes of length 1 to
represent the extra WCS axes is too restrictive, and that the alternative
technique of permitting the WCS to have a higher dimension than the image
should also be permitted.  IRAF can support either one.  Although the
WCSDIM keyword may not be strictly necessary, and the the axis mapping
idea  beyond the scope of the current paper, it should not be foreclosed on. 
A string is a very tidy way to encode an axis map.


1.5 Pseudo N-Dimensional Images

    I do not like the FITS convention where pseudo 3D and 4D images are
created solely for the purpose of dealing with wavelength, polarization, etc
definitions. This convention appears to be entrenched in the radio astronomy
community and we can support it in IRAF, but I think it should be viewed
as a convention not necessaritly the way to encode position, wavelength, or
even time dependent quantities which are really parameters not coordinates.
The optical people have in general not adopted this convention to define an
effective wavelength for direct images, the coordinates of a spectrum, etc.

    Examples of this type of data mighgt be 1) the spectrum of a star
(which has a fixed position on the sky) and a perfectly simple linear
spectrum. I think that the current FITS convention suggests that this data set
should be encoded as a 3D image with the lengths of the second and third
axes set to 1. However there is no projection of scaling involved here and I
doubt most optical astronomers would appreciate their image becoming 3D in
pixel space just in order to define the RA and DEC, or even 4D data if we
want to start including time as another dimensions.

    I don't really have a good alternative suggestion for dealing with this
since in the case case of equatorial coordinates it is desirable to
know things like RADESYS, EQUINOX, MJD-OBS, DATE-OBS etc even though the
CD matrix is meaningless. In fact all you really want is the CRVAL ...
One possiblity is to set all the CD matrix elements to 0 in this
case and use the CTYPE, and CRVAL to define the extra coordinates. This
would be ok as long as the dimension of the image does not have to be
the same as the WCS, but rows of 0's produce singular matrices when you
try to invert the WCS, and this condition would have to be tested for.

     In any case I think this is another of those conventions that needs
to be discussed more.



1.6 The IRAF Lterm (Compatibility with FITS)

In IRAF we support a builtin coordinate system which transforms from physical
pixels to logical pixels. For most astronomical images this is a 1 to 1
transformation. However for images with a history, i.e. they are subsections of
larger parent images, or are linear transformation from some larger image,
then the physical coordinates are pixel coordinates in the parent image and
logical coordinates are pixel coordinates in the parent image.


In IRAF this information is represented in the following way

l = LTM * p + LTV

where l are the logical coordinates, p are the physical coordinates ,
LTM is a linear transformation analogous to the CD matrix, and LTV is
the offset.

We can fit this scheme into the FITS framework by adopting one of the
multiple WCS, say p, defining the WCSNAME keyword value to be "physical",
and inverting the above transformation, e.g.

p = LTM(-1) * (l - LTV)

where LTM(-1) is the inverse of the LTM matrix. In this case the appropriate
FITS values are

CTYPEjp = 'linear  '
CRVALjp = 0
CRPIXjp = LTVj
CDj_ip = LTMj_i(-1)

This would be a IRAF convention layered on top of the FITS convention and
the key to it would be seeing the WCSNAME = 'physical' keyword = value
pair.



1.7  Function / Projection Parameters (What is and is not a parameter ?)

When I look at all 3 WCS papers taken as a whole I am concerned that there
seems to be a false distinction between what are called function parameters
the PVj_Ms and other parameters, including such things as LONPOLEs,
LATPOLEs, RESTFRQs, RESTWAVs, FPRADIUS, etc. I take parameters to mean
anything that is required to actually compute the coordinates of a particular
pixel in the image using the function defined by CTYPEjs. By that definition
LONPOLEs and LATPOLEs are parameters since their default and/or actual
values depend on which projection function is chosen. Similarly RESTFRQ and
RESTWAV are used by some of the spectral functions. Similarly with the
FPRADIUS parameters. To me this is going outside the FITS framework and
introducing yet more arbitrary parameters. I think this tendency should
be discouraged.


1.8 Alternate FITS Image Representations

Appendix C discusses alternate FITS image representations including
a multi-dimensional vector in a single element of a FITS binary
table, and a tabulated list of pixel coordinates in a FITS ASCII or
binary table.

I have not looked at this in detail and I have not worked much with
FITS tables except to say that it is messy to have 3 conventions for
representing essentially the same thing. Given the nature of FITS
there may be no way around it though.




		    Comments on the FITS WCS Papers
			     To Doug Tody
			  From Lindsey Davis
			      April 2001

	       [Revised (Minor Revisions Only) June 2001]


2. Paper II


2.1. LONPOLEs and LATPOLEs (Are they really different from the PVj_ms ?).

These parameters are used by all the celestial coordinate functions to rotate
native to celestial spherical coordinates, so in a way it does make sense to
separate them from the PVj_m functions.  However it is also true that LONPOLEs
and LATPOLEs are are associated with a longitude or latitude axis or both. 
It is also true that the default values assumed by these parameters
are projection, i.e. CTYPEjs dependent. This means that if the parameters
are not defined, then the function driver must assign or compute the required
value, just as it reads in the PVj_n parameters.  I don't see any essential
difference between LONPOLEs / LATPOLEs and the other PVj_ms parameters which
are also function specific. 

2.2 Default values for LONPOLEs and LATPOLEs (Documentation suggestion)

These quantities are well described in section 2. As the default values
of these quantities depend on the projection class, e.g. zenithal, cylindrical,
or conic it might be useful to list the default values where appropriate in
Table 7.


2.3 CTYPEjs Changes (Minor syntax changes that may matter in future)

Section 3 suggests that in order to avoid exceeding the 26 character limit on
longitude and latitude pair names the following syntax

CTYPEjs = 'xyLN-TAN'
CTYPEjs = 'xyLT-TAN'

should be supported as well as the currently form 

CTYPEjs = 'xLON-TAN'
CTYPEjs = 'xLAT-TAN'

This should not cause any problem for IRAF and is just a heads up so that we
don't overlook this point when when upgrading the IRAF / MWCS code.


2.4 LONPOLEs and RADESYSs (LONGPOLE and RADECSYS already in use)

Some FITS writers are already using LONGPOLE and / or RADECSYS so FITS
readers may need to be aware of LONPOLE/LONGPOLE and RADECSYS/RADESYS.
In IRAF we have used RADECSYS (the keyword defined in the Hanisch and Wells
paper) but have not used LONGPOLE. I believe the ROSAT people have also
used RADECSYS.

This is not a problem for IRAF, merely a heads up so that we don't forget this
point when upgrading the WCS code.


2.5 RADESYSs Applied to Ecliptic Coordinates (Currently used only for
    equatorial coordinates) 

The original Hanisch and Wells paper specified that RADE(C)SYS applied
only to equatorial coordinates. This is actually implied by the name.
They also implied that the ecliptic coordinates were specified by
MJD-OBS or DATE-OBS not EQUINOX. The current standard extends the RADESYSs
and EQUINOX definitions to the ecliptic system as well. 

This note is just a heads up that there was a change between the two papers.
This change does not affect the IRAF WCS code but does imply that the package
of sky coordinates routines we support will have to be modified. In restrospect
maybe we should have thought about this ourselves but ...


2.6 Clarify the Handedness Issue (Paper is a bit confusing)

Section 3.0 paragraph 2 implies that new celestial coordinate systems
must have the pole at +90 and be right-handed. Is this intended to apply only
to planetary coordinates, i.e. you are outside the sphere looking down or
is it more general? Section 4.  paragraph 5 seems to imply that pixel
coordinates are right-handed and and the (x,y) coordinates are left-handed.
Doesn't this mean that the celestial coordinates are left-handed as well? Is
there a real discrepancy or is the paper different issues in the two places.


2.7 Description of TAN projection

In paragraph 11 I think the  expression

(xi,eta) = (x*r**n, y*r**n)

should be

(xi,eta) = (x*r**n-1, y*r**n-1)

It might be worth mentioning that although the (x,y) coordinates are
zero at the reference point the (xi,eta) coordinates may not be because
the (a0,b0) coefficients are not necessarily 0. Therefore the CRVAL are
not the coordinates of the CRPIX. This makes the TAN projection different
from all the other celestial projections (I think). Of course it may be
possible to define the projection so that the a0,b0 are 0.0,0.0 but this
is not required. There is nothing wrong with this, but a note to the
reader might be useful.

I do object to the idea that the units of (x,y) may be anything
other than degrees. This breaks the celestial coordinates model and means
that the units of the CD matrix are not defined by CUNIT. The higher order
solution should be part of the projection just as it is for the ZPN projection,
not part of the CD matrix. In the proposed scheme how are we supposed to
know what the units of x,y are except via more special keywords ?

If "facilitate the translation of existing plate solutions" means the DSS,
it appears to me that the DSS plate solution does fit into the existing scheme
of (x,y) in degrees. It should be easy to preserve the (x,y) in millimeters
options by assigning that transformation to one of the multiple WCS
where units etc can be assigned properly and the model is not violated.

However it is possible to appropriate the idea of the lambda parameter
for another purpose, i.e. applying an optional normalization (this
is based on an suggested by Frank) . This might make evaluating the
distortion term better behaved in some better behaved in some circumstances.
If adopted there would be two lamda parameters one for x and on for y,
with default values of 1.0. The lamda parameters would then become
PVj_ms quantities.

After reading through the spectral coordinate paper and noting the use
of the same polynomial distortion function for all the spectral functions,
I would like to add one further comment / suggestion. If the first say 12
PVj_ms quantities are reserved for the 'projection parameters' then the
remaining PVj_ms could be used to support the same DSS style polynomials
for all the projections. Even doesn't wish to add this option to Paper II
defining the PVj_ms this way would leave the possibility of supporting
a more general polynomial term easier and still leave around ~50 parameters
at least to play with.


2.8 Description of ZPN Projection

I think the paper meant to say that PVj_1s has a default value of 1.0
and all the other PVj_ms take on a default value 0.0.


2.9 New FLS Projection (was GLS)

This is just an IRAF heads up that we should support the old projection and the
new one. We may not have code for the old GLS projection anymore?


2.10 Description of Quad-cube Projections

I have a comment about the 3D image and CTYPEjs = 'CUBEFACE' parameter
convention that is proposed to suport the quad-cube projection. The authors
state that this format was adopted for storage efficiency. Would it not be just
as easy to store this kind of projection in an MEF file using EXTENSION
= 'IMAGE   ' which produces the same gains in pixel storage, removes the
need for the CUBEFACE CTYPE value, while preserving the origins of each
cubeface in a more natural way then the current scheme?
We may have to support the 3D version for historic reasons but I think the
alternative is more natural. If this is true then this section goes beyond
defining the function and actually defines a convention for representing
this type of data.


2.11 Supporting Old Interpreters (Mixing CDELT/CROTA and CD elements)

Section 5.2 states that that FITS writers must not mix CDELT/CROTA
with the CD matrix. We follow this convention except in the 2D case that
CROTA2 = 0.0 where both CD and CDELT may be written, and in the 1D case
where for historic reasons CD1_1 and CDELT1 may be written. I prefer the
"don't mix conventions" option myself so I don't have a problem with
changing this. I don't know whether this affects any old spectrocopy code
me have though.

2.12 Pseudo N-Dimensional Images (Header Interpretation Example 1)

Section 6.3 shows how to construct a header for a simple 3D image which has
been made into a pseudo 4D image to support the Stokes parameter. Although we
can support this type of image fine in IRAF I think the community should be
discouraged from using coordinate axes to represent what are really
parameters for the whole image. An analogy would be to make every direct
image that came off an optical telescope a pseudo 3D image with the 3rd
axis being the effective wavelength of the filter, or making every spectrum
a 3D image with ra and dec of the object observed the 2nd and 3rd axis.
The optical community has in general not used the extra axes conventions
and with NVO issues coming up I think using this convention to represent
parameters should not necessarily be addopted at the only one. One can start
extending it indefinitely to all kinds of other quantities, temperature, etc.

I don't mean to imply that there is anything wrong with this convention but
I think it should be viewed as a convention not a standard.


2.13 Suggestion for New Header Interpretation Section (Quad-Cubes)

I would really like to see a whole sky quad-cube example demonstrating the
use of CTYPEjs='cubeface' and the path to evaluating a point on a arbitrary
face shown.  A parallel example showing how to do the same thing using
EXTENSION = 'IMAGE   ' would also be useful. I think the second option
is somewhat simpler but I may be misuderstanding the function. These
projections seem to me to be the most complicated  internally, so a
concrete example would be useful. 


2.14 DSS Header Construction Example

With respect to section 1.7 above I think evaluating x,y in anything other
than degrees should be discouraged for celestial coordinate projections.


2.15 Example 6.7

I am glad this example is here because what it implies it that it is not
possible to fully support dimensional reduction of images in FITS, i.e.
FITS does not support cases where the dimensionality of the WCS is greater than
that of the image and the coordinate axes are not orthogonal. This issue was
raised in the first review paper because in IRAF we can support this, i.e.
extraction a line from a 2D tan projected image, and maintain knowledge of the
coordinates, or extract a 2D image from a 3D cube and maintain a knowledge of
the coordinate system of the original image. This is done via the keywords
WCSDIM which defines the dimensionality of the  WCS and an axis map keyword
which describes which axis in the current image mapped to which axis in a
parent image.


2.16 Suggestion for Header Construction Section (Simple Coordinate Conversion)

I think it might be useful to include a simple example of a TAN projection
celestial coordinate conversion from equatorial to galactic coordinates that
is performed using 1) a rotation of the CD matrix and 2) the LONPOLEs
parameter which fully illustrates the coordinate conversion recipe at the
end of section 2.



                      Comment on the FITS WCS Papers
			       To Doug Tody
			     From Lindsey Davis
			        April 2001 

		   [Revised June 2001 (minor edits only)]

Paper III


3.1 Higher Order Corrections (Absence of cross-terms)

I think the higher order correction representation given in section 2.2
is too limited. First and most important it does not permit cross-terms or
terms dependent on other axes which I think is a problem. How are you
supposed to represent S-distortion and functions like that. As long as support
for such an option with the current framework is not foreclosed on, we can
defer that discussion. However it appears to me that the comments in section 2
imply that the authors feel that these issues are dealt with in sections 4
and 5. After reading these sections I am not at all sure this is the case.
(more comments later)

Along these lines I note there is no model schematic for the spectral WCS
as there is in papers I and II (Figure 1 in both cases).


3.2 Higher Order Functions / Corrections (no provision for alternate forms ?)

It appears that only one higher order correction function, 1D polynomials,
is defined and this function is assumed to be the same for all
values of CTYPEjs. Aside from the 1D limitation this may be OK as long
as other functional forms can still be represented. In a way I  find the use
of polynomial distortion functions for all the spectral functions "interesting"
as when I originally suggested doing something similar for the celestial
coordinates I did not get much support. I still think that supporting a
DSS style polynomial (although not necesssarily in detail) for all the
celestial projections is a good idea.

One question I do have in connection with this issue is how we are going to
support alternate functional forms. Let us say we a have a spectrum
that is in wavelength units that is almost linear but with a sinusoidal
distortion term. Can I now define a new spectral functions which has
CTYPE='WAVE-TRG' where TRG stands for trigometric and some of the
PVj_s can be defined to be the trig function parameters, e.g. amplitude,
periodicity etc. If so the polynomials must be no longer valid for this
function or in a set of reserved PVj_ms keywords. Or is it intended that
the functions defined in section 3 somehow be a complete set?


3.3 Coordinate Parameters

To me the RESTFREQs and RESTWAVs parameters defined in section 3.3
are analogous to the LONPOLEs and LATPOLEs parameters for celestial
coordinates. They are required to compute spectral coordinates and they
are "projection" dependent in that they are required for some values
of CTYPEjs but no for others. I don't see why they are not represented
by PVj_s keywords in both cases.


3.3 Coordinate Reference Frames

As a non-spectroscopist I found this section somewhat confusing. After
reading it a couple of times I finally (think I) got the general idea that
SPECSYSs unlike RADESYSs is not a passive parameter defining the current system,
but an active parameter describing a transformation which may be but is
not necessarily the identity transform. If the two parts of SPECSYS are the
"same" then the spectral coordinates are correct over the image plane and
require no modification. If they are different then there are systematic
errors over the image plane that must be corrected for and this correction
is outside the purview of the usual FITS interpreter, i.e. it is not encoded
in the FITS header explicitly but instead requires specialized spectroscopy
knowledge and the parameters in Table 2. Is this really the case ? If so it
would be like requiring simple programs which list coordinates, label simple
graphs, etc to know all about spectral coordinates to do simple things
correctly. Again maybe this is a misunderstanding on my part but I think
evaluating correct coordinates and interpreting there meaning should be
separable functions as is the case for the celestial coordinates ...

Again this may be a naive question but is there a position on the sky to go
with the proposed VSOURCE / ZSOURCE values ...

3.4 Spectral Coordinate Variation with other Coordinates

I don't see why a there needs to be yet another special parameter in
this case FPRADIUS. Can this not be defined via a PVj_ms parameter? It is
a "projection" parameter and that is what these parameters were intended
to be used for.


Date: Fri, 30 Mar 2001 13:08:29 -0700 (MST)
From: Frank Valdes <valdes>
Message-Id: <200103302008.NAA01892@tucana.tuc.noao.edu>
To: davis
Subject: Spectral WCS comments

Lindsey,  look this over and incorporate or pass it on as necessary.
If something is unclear or unconvincing I can refine it.  If you disagree
we can talk about it further.  Frank

PS Note that I have dropped my concerns about physical coordinates, etc.,
since the CD transformation provides this and subraster extractions,
binning, etc can be applied to the CD matrix without requiring the
PV coefficients to be recomputed.  Also note that by defining CRPIX
to be the middle of the fitted range and PVN to be 1/(CD*2*(xmax-xmin))
gives the same (-1,1) normalization we use in the curfit package.



Comments on "Representation of spectral coordinates in FITS"

The "Representation of spectral coordinates in FITS" is very good for
describing spectral coordinates along dispersion axes which are "decoupled"
from other axes after application of the linear CD transformation.  The
usage for this type of decoupled dispersion axis has wide applicability
such as for extracted one dimensional spectra and geometrically rectified
long-slit and spatial/dispersion data cubes.  The question of non-linear
axis coupling is not really worked out.  Being able to attach a coupled
system to optical data is a low priority.  Therefore my first suggestion is
to simply state that the paper deals with decoupled dispersion data axes
and defer situations with more complex coupled spatial and dispersion axes
to a later stage.

I like the organization of the CTYPE definitions.  It has nice features
for converting between the various dispersion coordinate types without
regridding.

In my experience with optical spectral data and spectral data reductions,
one first determines a 1D non-linear dispersion function (the analog to an
astrometric solution in celestial coordinates) and then optionally regrids
to a simple sampling; primarily WAVE-WAV/WAVE-LOG or FREQ-FRQ/FREQ-LOG.
The velocity types may be used for display but not too commonly as a
transportable/archive format.

It is in the area of non-linear 1D dispersion functions that I have some
suggestions.  The concept of dealing with this in the PV correction is good
and I like this better than trying to define things in the CTYPEs.  There
are two problems I see with the PV prescription.  One is the lack of a
normalization for the intermediate coordinates and the other is the lack of
different functional forms.

One would like the CD transformation to produce approximate linear
dispersion coordinates to which small corrections are applied.  But these
intermediate coordinates can be in a wide variety of units such that values
at the extremes of the data from the pixel reference might have values from
10E-10 to 10E+10.  Applying any polynmial function to such values can be
numerically difficult.  For example if you wanted a function up to a power
of 10 then you would be computing numbers between 10E-100 and 10E+100.  The
way to deal with this is to allow a normalization of the intermediate
coordinate.  Of course one could do this with the CD terms but then the
uncorrected intermediate coordinate would no longer be the desired linear
approximation to the final coordinates.

So I suggest defining a numeric normalization keyword such as PVNks.  Then
the independent value for the functional evaluation would be (w/PVNks)
instead of just w in the notations of the paper.

Using a normalization would make the power law functional form of equation
(9) be reasonable.  But why limit the formalism to one functional form?
Polynomials are good but other functions can have better physical meaning
or stability.  For example one might want to use a form appropriate for
a grating (a trigonometric function) or use locally lower order representation
as provided by splines.

What I suggest is defining a string-valued function keyword such as PVFks.
Some functions can be defined for this paper and others can be added latter
by the usual augmentation methods by which the FITS standard evolves.  The
function types that I would be interested in are "POWER" (the current
equation 9), "CHEBYSHEV" or "LEGENDRE", and "SPLINE3" (cubic spline).
Dispersion functions derived for optical/IR spectra in IRAF are about
equally CHEBYSHEV/LEGENDRE and SPLINE3.

One minor comment is that I don't see why the PVk_ms coefficients have
to restrict m to the range 0-99.  Most data have axes numbers below 10
and I have never encountered real data with dimensionality above 4.  If
the requirement is made that the combination of k and m have less
than 4 digits then for most practical situations the functional forms
could have up to 999 coefficients.

For your consideration,
Frank Valdes



Date: Wed, 11 Apr 2001 18:15:46 -0400
From: Don Wells <dwells@cv.nrao.edu>
To: Doug Tody <tody@tucana.tuc.noao.edu>
Subject: Re: FW: WCS issues

Doug, (and colleagues)

Doug Tody writes:
 > ..I have gotten feedback now from both Lindsey and Frank..  We
 > still need to distill these down to the minimim we think is
 > necessary and pose things in such a way that we hopefully won't
 > inflame the issue..

Thanks for forwarding the drafts.

 > ..relationship of IRAF WCS and the FITS WCS, .. an internal issue
 > ..but a useful test case..

Right.  This is exactly why I have waited for your Group to reach this
point. 

 > 1. Paper I
 > 1.1. Default Values of the CD matrix (FITS convention / IRAF convention)
 > Section 2.1 and Table 1 define the default values of the CDj_is
 > matrix to be 0.0 if i != j and 1.0 if i = j.. if the off-diagonal
 > elements are defined then the diagonal elements default to 1.0
 > which is not the convention adopted by IRAF.  .. the problem boils
 > down to different conventions and we need to decide which is
 > better. The new convention causes portability problems for us..

This question can be raised in the committees if the IRAF Group thinks
that the proposed convention should be changed.

 > 1.2. Default values for CRPIX and CRVAL (FITS convention ?)  
 > As far as I can tell default values for CRPIX and CRVAL are not
 > defined in section 2.1 and are not listed in Table 1..  IRAF
 > .. convention [is] that these quantities default to 0.0 if undefined..

I expect that Eric can add the appropriate sentence(s).

 > 1.3 String Parameter Values (Better Interface / Compatibility with IRAF)
 > There is no support for string valued functon parameters. The PVj_ms are
 > by definition floating point numbers..

Although I am unable to find a statement of the rule in NOST_2.0, it
appears to me that there is a consensus in the FITS community that the
value fields of all instances of each reserved keyword should have the
same type.  Therefore I do not expect that community will ever approve
using strings in addition to numeric values with the PVj_ms keywords.

 > .. in IRAF we have used builtin string parameters for coordinate
 > type (axtype), projection type (wtype), units(units), label(label),
 > and format (format) quantities, and function specific string
 > parameters for encoding a variety of arbitrary functions including
 > polynomials, splines,.. I think the standard would be significantly
 > enhanced if string parameters in which any quantity, characters,
 > integers, reals, keyword=value pairs could be encoded were
 > supported. What the quantities mean would be determined solely by
 > the function type defined by the CTYPEjs parameter which need not
 > be one of the standard functions but may be a personal or user
 > function. We already do this in IRAF..

Obviously you need some semi-standard mechanism to interchange your
proprietary info between IRAF installations.

 > One possible way to implement this that might be more acceptable to
 > the current FITS community than the IRAF scheme, is to adopt a set
 > of string parameters WVj_ns ..  Adopting an WVj_n parameter does
 > not affect the standard or standard functions in any way but it
 > does provide more flexibilty to the scheme..

This is a reasonable proposal.  The string parameters could be used
for any purpose.  In order to protect your own usages, you should
probably prefix your string values with something like '<IRAF>'.  It
is possible that the community would accept this generic string
keyword as an addition to the reserved keywords, but it might also
turn out that there would be reluctance to make it an official part of
the interchange Agreement.  This sort of thing has happened several
times before in our negotiations, and we have developed a solution:
the unofficial appendix!  Indeed, even though such appendices always
carry disclaimers that they are not part of the formal standard and
are included "for informational purposes only", it is a fact that all
of the appendices devised so far are reproduced in the NOST document
(see Appendices B and D).  I recommend that the IRAF Group draft such
an appendix to be added to WCS-I.

 > 1.4 Dimensional Reductions and the WCS (FITS conventions / IRAF conventions)
 > ..what happens to the WCS if you extract a line from a 2D image, a
 > 2D plane from a 3D cube etc? ..impression.. that in
 > FITS.. dimensional reduction is not allowed, i.e. you can only
 > extract a 2D line with the length of the second axis set to 1 from
 > a 2D image. Is this true?

This is how AIPS has always done it, so it is no surprise that Eric
composed the draft standard in that fashion.  

 > Does this mean one a can never represent a true 1D image in FITS by
 > extracting a line from an image with a celestial coordinate wcs
 > where the coordinate axes are dependent and information for both
 > coordinates is required?

Probably I don't understand the issue that you are raising.  If by
'1D' and 'line' you mean that the 1D is something like a crossection
of objects in the scene at some position angle and with some sample
spacing (a 'slice' in AIPS-speak), why wouldn't NAXIS=2 be a
completely feasible convention?  I will be surprised if you can
suggest any NAXIS=1 scheme for this case that will be better than
NAXIS=2 in terms of generality and economy of notation (as measured by
negotiation effort).

 > I would like to suggest an alternative convention that can be
 > equally acceptable and that is to permit the dimension of the WCS
 > to be higher than that of the image encoded in NDIM.. I see no
 > reason why this convention can not be as valid a convention as the
 > current one of maintaining image axes of length 1.

Occam's Razor.   The actual bits in the data array will be identical,
the information content of the keywords will be identical, the only
difference is the addition of one more keyword in your notation.

 > ..In IRAF..  axis map keyword WAXMAP.. defines how the WCS axes of
 > the parent image map to the current image axes..

Although the draft proposal discourages the usage, it is a fact that
the CD keywords can be used to permute and flip axes with complete
flexibility. 

 > .. I think there should be some discussion of this issue in the
 > first paper.. I think the technique of using extra axes of length 1
 > to represent the extra WCS axes is too restrictive and the
 > alternative technique of permitting the WCS to have a higher
 > dimension than the image should also be permitted..

Maybe I am missing something here -- how can the extra-WCS-axes
approach add any functionality that is not also in the
axes-of-length-1 method?

 > 1.5 Pseudo N-Dimensional Images
 > 
 >     I do not like the FITS convention where pseudo 3D and 4D images
 > are created solely for the purpose of dealing with wavelength,
 > polarization, etc definitions.. I think it should be viewed as a
 > convention not as THE way to encode position, wavelength or even
 > time depedent quantities which are really parameters not
 > coordinates..

For radio astronomy the 3D/4D coding is a compact generalized notation
that tells the truth about the data. The VLA and other synthesis
imaging telescope are inherently four-Stokes-parameter multi-spectral
systems.  Continuum imaging is regarded as a special case of spectral
line imaging, and is treated as such by much of the data processing
and analysis software.  Unpolarized imaging (Stokes-I only) is a
special case of polarized imaging.  It is true that these special
cases are _very_ popular (involving the [vast] majority of the
datasets), and that that fact could be used to justify some special
notation in our interchange format, but our designers think that it is
better to use the same generalized notation for all radio observations
whereever it is possible.  An interchange standard for all of
astronomy should also use a generalized notation.

 >  The optical people have in general not adopted this convention to
 > define an effective wavelength for direct images, the coordinates
 > of a spectrum, etc.

but I say that they should..

 > Examples of this type of data mighgt be 1) a spectrum of a star
 > (which has a fixed position on the sky) and a perfectly simple
 > linear spectrum.  ..the current FITS convention suggests that this
 > be encoded as a 3D image with the lengths of the second and third
 > axes set to 1..

Occam's Razor again: negotiation of more keyword agreements and
composition of more descriptive text will be required so that a 1D
spectrum of a point source can be interchanged with a different coding
convention than is used for interchange of a 2D long-slit spectrogram.
The WCS-I architecture does both cases with the same notation, without
any special rules to distinguish the 1D and 2D cases.

 > In any case I think this is another of those conventions that needs
 > to be discussed more.

I agree that the generalized axes concept needs to be stressed more,
as an antidote to 'photonic provincialism'.

 > 1.7  Function / Projection Parameters (What is and is not a parameter) ?
 > 
 > When I look at all 3 WCS papers taken as a whole I am concerned
 > that there seems to be a false distinction between what are called
 > function parameters the PVj_Ms and other parameters, including such
 > things as LONPOLEs, LATPOLEs, RESTFRQs, RESTWAVs, FPRADIUS, etc.
 > ..LONPOLEs and LATPOLEs are parameters since their default and/or
 > actual values depend on which projection function is
 > chosen. Similarly RESTFRQ and RESTWAV are used by some of the
 > spectral functions..

I too am somewhat uncomfortable about the proliferation of these
'special' keywords. However, it appears to me that all of them are
specified to have the 's' suffix.  This means that a FITS header can
contain two or more WCS specifications with differing values of these
'parameters'.  The only limitation seems to be that any single WCS
specification cannot have two values of LATPOLEs or two values of
SPECSYSs, but I am unable to imagine any use for two LATPOLEs or two
SPECSYSs values in a WCS, so I guess it is OK.

 > 3.1 Higher Order Corrections (No provision for alternate forms ?)
 > 
 > I think the higher order correction representation given in section
 > 2.2 is too limited. First of all it does not permit cross-terms or
 > dependent axes which I think is a problem. However as long as
 > support for such an option with the current framework is not
 > foreclosed on we can defer that discussion.

I agree. Section_2.2 needs to be generalized to cover more cases. It
should be completely capable of representing long-slit spectra, and if
possible it should represent the imaging Fabry-Perot/Michelson data as
well. 

 > 3.3 Coordinate Reference Frames
 > 
 > .. I think it is true that SPECSYS is intended to be like RADECSYS
 > a "passive" parameter.  If the two parts of SPECSYS agree then the
 > spectral coordinates are correct over the image plane and require
 > no modification. If they are different then there are errors over
 > the image that must be corrected for, and this correction is
 > outside the scopy of the FITS interpreter. Is this the case ?

I think so.  The SPECSYSs definition is a new thing, probably not yet
implemented in any software system, so we all have some learning to
do.

 > .. is there a position on the sky to go with the VSOURCE / ZSOURCE
 > values ...

The text should discuss this issue, which is important for proper
interpretation of interchanged data.  For example, if optical and radio
spectra are compared, it makes sense to do the comparison of line
profiles in _velocity_ units with the same reference frame. 

 > 3.4 Spectral Coordinate Variation with other Coordinates
 > 
 > I don't see why a there needs to be a special FPRADIUS
 > parameter. Can this not be defined via a PVj_ms parameters. It is a
 > projection parameter of a sort.

Good question.   I am uncomfortable with the current version of
section_5 of WCS-III.  It needs more work.  The FPRADIUS question
should be reviewed by Bob Schommer and his colleagues in the imaging
Fabry-Perot community. Steve Ridgway could probably comment on the
imaging Michelson case. We need some sort of new generalization...

 > .. My biggest concern is that there does not seem to be
 > well-defined model for spectra like there is for celestial
 > coordinates, and there is a tendency to introduce ever more "adhoc"
 > keywords and tricks ...

Right.  WCS-III needs some more work, and the relevant astronomical
specialists need to take it seriously.  My concern with this is that
it might lead to requirements for change in WCS-I.  If so, we need to
know about that _now_!  Our community could afford to defer agreement
on WCS-III itself for quite a while, but we need agreement on WCS-I
and WCS-II ASAP.

 > .. represent S-distortion and functions like that..

Are there still any examples of instruments that exhibit S-distortion?
As far as I know, the only production instruments of that type in
recent years were on HST, probably only the FOC.  Do you know of
any others?  

 > .. spectrum.. in wavelength units that is almost linear but with a
 > sinusoidal distortion term. Can I now define a new spectral
 > functions which has CTYPE='WAVE-TRG' where TRG stands for
 > trigometric and some of the PVj_s can be defined to be the trig
 > function parameters, e.g. amplitude, periodicity etc...

Obviously the architecture permits such a definition, but why bother?
The idea is to fit the data in your datasystem with any formula that
suits your taste, but then to convert your formula to the equivalent
polynomial for interchange.  If this procedure is mathematically
satisfactory, then in my guise as Chair I will argue strongly that we
adopt it in order to get an interchange agreement with minimum
negotiating effort. 

 > 3.3 Coordinate Reference Frames
 > 
 > .. SPECSYSs.. 

Probably Eric should add a numeric example to illustrate why SPECSYSs
is needed for imaging spectroscopy systems.  The basic point is that
the various 'systems' involve adjusting the spectral axis by the
velocity vectors of the observer wrt the systems projected into the
direction of the target.  An imaging spectroscopy system observes
targets in multiple directions, so its data need different adjustments
in different parts of the field.  Eric's draft may be the first public
discussion of this detail...

 > From: Frank Valdes <valdes>
 > Comments on "Representation of spectral coordinates in FITS"
 > 
 > The "Representation of spectral coordinates in FITS" is very good
 > for describing spectral coordinates along dispersion axes which are
 > "decoupled" from other axes after application of the linear CD
 > transformation.. my first suggestion is to simply state that the
 > paper deals with decoupled dispersion data axes and defer
 > situations with more complex coupled spatial and dispersion axes to
 > a later stage.

We could do that.  I would prefer to attempt to get a more general
agreement, but I would settle for solving the decoupled interchange
problem. 

 > ..[optical astronomers regrid] to a simple sampling; primarily
 > WAVE-WAV/WAVE-LOG or FREQ-FRQ/FREQ-LOG.  The velocity types may be
 > used for display but not too commonly as a transportable/archive
 > format..

Ditto for radio people, but when the astrophysical implications are
being discussed it usually has to be done using a velocity scale.
Currently as a practical matter astronomers cannot interchange spectra
in the velocity form, so the comparisons are somewhat ad hoc.  I hope
that the WCS agreement, by enabling us to convey two or more WCS
representations in headers and by codifying conventions for the
velocity form, will promote and encourage more frequent combining of
spectral data from the different observing disciplines.

 > ..intermediate coordinates can be in a wide variety of units such
 > that values at the extremes of the data from the pixel reference
 > might have values from 10E-10 to 10E+10.  Applying any polynmial
 > function to such values can be numerically difficult.  For example
 > if you wanted a function up to a power of 10 then you would be
 > computing numbers between 10E-100 and 10E+100.  The way to deal
 > with this is to allow a normalization of the intermediate
 > coordinate..  I suggest defining a numeric normalization keyword
 > such as PVNks..

Several years ago (at the Sonthofen BOF I think), I argued for
interchange of Chebyshev polynomial coefficients plus a normalization
dimension, but got no support from other negotiators.  Obviously one
should never do a high-order _fit_ with unnormalized polynomials, and
orthogonal polynomials (Legendre, Chebyshev) are mathematically
superior to ordinary polynomials for _fitting_, but _interchange_ of
equivalent ordinary polynomial coefficients with double precision
E-notation value fields will probably work just fine. Do you have any
reason to think that _interchange_ of ordinary polynomial coefficients
is mathematically unsound?

 > ..I would like to add one further comment / suggestion. If the
 > first say 12 PVj_ms quantities are reserved for the 'projection
 > parameters' then the remaining PVj_ms could be used to support the
 > same DSS style polynomials for all the projections. Even I you
 > don't wish to add this option to Paper II defining the PVj_ms this
 > way would leave the possibility of supporting a general polynomial
 > term easier and still leave around ~50 parameters at least to play
 > with.

an interesting idea.

				 -=-

I want to emphasize the question of whether WCS-I (overall
architecture of the notation) needs change.  If study of WCS-III
suggests a need for revision of WCS-I, please bring this to the
attention of the community ASAP.  If you think that WCS-I is OK,
please annouce this so that the community can proceed with the formal
review and voting process.

Regards,
Don


Date: Wed, 11 Apr 2001 16:33:22 -0700
From: Doug Tody <tody@noao.edu>
To: Don Wells <dwells@cv.nrao.edu>
Subject: RE: FW: WCS issues

    [ The following text is in the "iso-8859-1" character set. ]
    [ Your display is set for the "US-ASCII" character set.  ]
    [ Some characters may be displayed incorrectly. ]

Hi Don -

Thanks for the detailed comments.  I was not expecting such detailed
comments at this stage, before we have our internal discussions, but
they are useful and should help speed things along.
 
> I want to emphasize the question of whether WCS-I (overall
> architecture of the notation) needs change.  If study of WCS-III
> suggests a need for revision of WCS-I, please bring this to the
> attention of the community ASAP.  If you think that WCS-I is OK,
> please annouce this so that the community can proceed with the formal
> review and voting process.

I understand...  I think it is necessary to look at the complex test
cases (such as spectral WCS, mappings into IRAF, etc.) to give WCS-I
a real test.  Thus far it has been used for little but sky projections.
We need to take a preliminary but serious look at these more complex,
non-sky projection WCS, see if WCS-I supports these, and then we can
standardize WCS-I.  Then we will go back and devise a new spectral WCS
standard and so on.

We are actively looking at the problem of actual FITS-based spectral
data formats and data models now and this will allow us to mentally
test the generality of WCS-I.  We will continue to look at the IRAF
issues as well (some changes to the way IRAF handles WCS are possible
here as well as possible changes to FITS WCS).

[...]

	- Doug



Date: Fri, 20 Apr 2001 16:02:00 -0400
From: Don Wells <dwells@cv.nrao.edu>
To: Doug Tody <tody@noao.edu>
Subject: RE: FW: WCS issues

Doug etal.,

Doug Tody writes:
 > > I want to emphasize the question of whether WCS-I (overall
 > > architecture of the notation) needs change.  If study of WCS-III
 > > suggests a need for revision of WCS-I, please bring this to the
 > > attention of the community ASAP.  If you think that WCS-I is OK,
 > > please annouce this so that the community can proceed with the formal
 > > review and voting process.
 > 
 > I understand...  I think it is necessary to look at the complex test
 > cases (such as spectral WCS, mappings into IRAF, etc.) to give WCS-I
 > a real test.  Thus far it has been used for little but sky projections.
 > We need to take a preliminary but serious look at these more complex,
 > non-sky projection WCS, see if WCS-I supports these, and then we can
 > standardize WCS-I.  Then we will go back and devise a new spectral WCS
 > standard and so on.

The key unsolved problem of WCS-III is the details of WCS-III.5
('Spectral coordinate variation with other coordinates'). In
particular, the major problem to be solved is the representation of
spectral coordinate variation with celestial coordinates.  

It appears to me that the latter problem can be solved by adopting a
variation on the PVi_ts/PVj_ts conventions which were adopted for TAN
(WCS-II.4.1.2).  I.e., spectral axis k would use a 40-term PVk_ts 2D
polynomial with a suitable range of t values to produce a spectral
variation with celestial (x,y) for the celestial axes defined in the
same header.  (Perhaps the PVk_ts polynomial should use the celestial
(xi,eta) rather than celestial (x,y); I leave that refinement as TBD.)
The PVk_ts values for t from 0-->m are reserved in WCS-III.2.2 for
correcting the spectral axis, so we need to offset t for the 40
additional PVk_ts values; I suggest that an offset of 20 should be
sufficient (total of 60 terms reserved).

The 2D polynomial set includes radial terms which can represent the
Fabry-Perot variation in Eq(55) in WCS-III.5. I expect that the CDi_j
matrix plus the 2D polynomial will prove to be strong enough for all
long-slit cases.

I am not aware of any need for a 3D (celestial+spectral) correction.
I.e. it appears to me that a separable correction of the sort that I
propose above (1D+2D) is sufficient for all known cases.

    +------------------------------------------------------+
    | I conclude from the above analysis that a sufficient |
    | solution for WCS-III.5 can be achieved with _no_     |
    | modification of the conventions proposed in WCS-I:   |
    | a solution using the PVk_ts terms alone will work.   |
    +------------------------------------------------------+

I ask that you (plural) consider the above conclusion.  If you agree
with it, tell me ASAP and I will announce it publicly and will ask the
community as a whole if anybody disagrees.  If the community achieves
a consensus on this point we will be free to adopt WCS-I, and probably
also WCS-II, almost immediately.

-Don



Date: Wed, 25 Apr 2001 15:20:49 -0700 (MST)
From: Frank Valdes <valdes>
To: dwells@cv.nrao.edu, tody
Cc: davis
Subject: Re: WCS issues

In this note I try to identify, from my experience, the issues associated
with spectral instruments where the dispersion coordinate has a spatial
dependence.  The ideas are based on the physical features of
instruments which use cross-dispersing elements and cameras that are found
from the far-IR through some high energy instruments.  I also include
filter scanning imaging, such as Fabry-Perot, since it is relatively
simple.

In the rest of this section I outline the components of defining a
spectroscopic WCS with some spatial component.  I am not doing this in any
FITS keyword model and leave it to later to determine how this would be
done.  In the next section I give some specific examples in case studies
for the common dispersing systems.  Finally I give some thoughts about
issues and answer Don's question about WCS-I.


First we define a coordinate reference pixel.  A natural reference point
would be the location of the zero order (undispersed) image at some spatial
point imaged on the sky.  For aperture spectra this might correspond to a
point in the undispersed image of the aperture.  This would either be the
center of the aperture or the where the center of the object falls in the
aperture.  For aperture-less 2D instruments, such as objective prisms, this
would be the zero order object positions, where multiple objects requires
multiple WCS.  However, see case 5 below for an application convention
based on a single spectral WCS.  For filter systems, such as Fabry-Perot,
this would be the position of the optical axis in the image and the
wavelength at that point of some exposure.

Next define the affine transformation so that one transformed axis is
tangent to the direction of constant dispersion coordinate at the reference
point.  The second axis is tangent to the direction in which the
undispersed reference point is dispersed.  With a single disperser the two
axes will generally be orthogonal but with cross-dispersion, such as in a
cross-dispersed echelle spectrograph, these will not be orthogonal.

Apply a general 2D distortion function to the intermediate coordinate.
This distortion function does not model or remove the shape of the zero
order entrance aperture, if there is one, or the dispersion relation.
Physically it corresponds to the distortion in the imaging system with the
dispersers and apertures removed.  It might be described in the same way as
an astrometric distortion function.  The effect of this function is to
correct for changes in the shape of the zero order spatial object (say the
entrance aperture) so that the shape is the same everywhere as it is
translated along the dispersion axis.

The effects of the dispersers and the shape of the entrance aperture are
applied with four functions.  The first is the correction to the dispersion
coordinate as a function of the cross-dispersion coordinate(s) to maintain
a line of constant wavelength.  For 2D aperture spectra this corresponds to
the cross-dispersion shape of the entrance aperture such as a curved slit
in a long-slit observation.  For a 2D filter scanning system, such as
Fabry-Perot, this corresponds to the shapes of the constant wavelength
surfaces in each filter image.

The second function is the 1D correction to the dispersion coordinate as a
function of dispersion coordinate needed to linearize the dispersion
coordinate.  This corresponds to the usual non-linear dispersion relation.

The third and fourth functions provide corrections to the cross-dispersion
coordinate as a function of cross-dispersion and dispersion coordinates.
In the absence of a cross-disperser these should be identify functions
since variations in the cross-dispersion coordinates should be due only to
the imaging distortion removed by the 2D function.  However with
cross dispersion they would have the same meaning as for the primary
disperser.

The last step is to apply the dispersion gridding function and reference
values.  This is basically what is described in WCS-III.

Multiple Transformations?

Of course one should be able to combine multiple transformations into
a single one.  However, I think it is very useful to understanding and
description to keep the terms somewhat separate.


Case Studies:

In the following case studies (x,y) refer to image pixel coordinates,
(x',y') are the transformed intermediate coordinates where x' is the
dispersion axis and y' is the spatial axis, (d,s) are the corrected
dispersion and spatial axis before applying the final coordinate reference
value and any projection or gridding transformations.  Functions pd(x',y')
are ps(x',y') are 2D polynomial functions for the dispersion and spatial
axes.  Functions fd(x'), fd(y'), fs(x'), and fs(y') are 1D functions for
the dispersion and spatial axes.  In this case the functions can include
cubic splines which have been found to be very useful for spectral
applications.

The three dimensional F-P case uses similar definitions as given in that
section.


Case 1: Long-slit spectroscopy

In the simplest case of a perfectly straight slit with undistorted camera
optics, the linear transformation corrects for rotation of the dispersion
relative to the image raster.  If the slit is not aligned with the
cross-dispersion axis this is also corrected using a skewed
transformation.  There is no imaging distortion correction (by definition)
and the only function applied is the 1D dispersion function accounting for
the non-linear dispersion of the dispersing element.

    d = fd(x')
    s = y'

The next most likely case is when there is some imaging distortion.  In
this case a 2D distortion function is is applied.  This describes all the
spatial variation, primarily the lines of constant spatial offset.
The 1D dispersion function is then added as before.

    d = fd(x') + pd(x',y')
    s = y' + ps(x',y')

The function ps(x',y') is what is sometimes refered to as the "S-distortion".
With just a single stellar observation the empirical determination of
this by tracing the profile peak reduces to fs(x').

If the slit is curved or accuracy requires this correction, then
in addition to any imaging distortion correction, the dispersion coordinate
is corrected by a 1D function of the spatial coordinate.

    d = fd(x') + fd(y') + pd(x',y')
    s = y' + ps(x',y')


Case 2: Multiobject and Integral Field Unit Spectra

The data consist of multiple spectra having a small spatial extent.  Each
spectrum has a WCS relating image (x,y) and dispersion and spatial
coordinate (d,s).  The linear transformation to x' and y' aligns the
overall dispersion and spatial axes of x' and y'.  Relative to this the
spatial coordinate origin, generally measured from the peak of the object
spatial profile, wiggles in some non-simple way.  The
mapping to d and s is given by

    d = pd(x',y') + fd(x')
    s = y' + ps(x',y') + fs(x') 

For aperture extraction of fibers and slitlets, where the extraction
aperture is described by a WCS and the extraction is over a
specified range of s, what has been used in IRAF for many years is:

    d = x' + fd(x')
    s = y' + fs(x')

>From experience I can say that the spatial offset function, fs(x'), is
generally not a simple polynomial and is best described by cubic splines.

The first generalization that would be desirable to make is that the
spatial scale, and so the aperture width, can vary due to focus
variations.  This corresponds to  the polynomial term in the s equation.  For
example a quadratic variation where the focus is best near the center of
the frame and is worse at the edges would be

    s = (a + by' + cy'^2)x' + fs(x')

The difference between this and the long-slit case is that because of
the small spatial extent one can replace the ps(x',y') with fs(x')
and use a more complex 1D function such as a cubic spline rather than
a polynomial.

Case 3: Cross-dispersed Echelle Spectra

Cross-dispersed echelle spectra, where there are many orders and a small
slit can be treated (and are treated by IRAF) in exactly the same way as
the multiobject case.  The only difference is that the transformation from
(x,y) to (x',y') includes a skew.

If one wanted to have a single WCS, as opposed to a WCS for each order,
then the form it would take is

    d(o) = (pd(x',y') + fd(x')) a/o
    s(o) = g(o) + y' + fs(x')

where 'a' is a coefficient giving the order at which the dispersion function
fd(x') is defined.  The function g(o) gives the order offset as a function of
the integer-valued order number.  The function would be based on the known
physical behavior of the cross-disperser, typically a prism or a
grating.  It would be up to an application to select the
order for some point (x,y) in the image.  This would generally be done by
finding the o which minimizes |s(o)| where s(o) = 0 at the order center
or the center of the source in the slit.

In principle, after the mapping to (x',y') the dispersion relation for
different orders is just an a/o scaling of a 1D dispersion relation at a
reference order.  However, the pd(x',y') term provides a small higher order
correction for grating and optical distortions.

Case 4: Fabry-Perot

The linear transformation is applied to (x,y,z) to produce (x',y',z').  In
this case let (x',y') be the spatial intermediate coordinates and z' be the
intermediate dispersion coordinate.  Then the celestial coordinates (l,b)
and dispersion d are obtained by

    l = pl(x',y')
    b = pb(x',y')
    d = pd(x',y',z') = fd(z') cos(arctan(R(x',y')/C))

The pl/pb would be distortion functions such as defined for the distorted
tan WCS, and the d shows an explicit function as given in WCS-III.  The
dispersion function has a 1D function for the dispersion (in wavelength) at
the origin of the (x',y') coordinates corrsponding to the optical axis on
the image.

Case 5: Aperture-less Spectroscopy

An example of this is objective prism imaging.  This would naturally
have (at least) two WCS.  One would be a celestial system giving the
zero order spatial coordinates.  The spectral WCS applies the
linear transformation converting to nearly orthogonal spatial
and dispersion coordinate relative to some the zero order
celestial coordinate.

Applications would need to interpret the WCS for some celestial
coordinate.

    (l,b) -> celestial WCS -> (x0(l,b),y0(l,b))
    (x(l,b),y(l,b)) -> spectral WCS with (x0,y0) -> (x'(l,b),y'(l,b))

    (x,y) -> spectral WCS ->  (xd',yd')

    d(l,b) = pd(xd',yd') + fd(x'(l,b))
    s(l,b) = ps(xd',yd')

The first step computes the zero order image position (x0,y0) for a desired
celestial coordinate (l,b).  This coordinate is used as the reference pixel
for the spectral WCS.  The second step creates an intermediate coordinate
relative to fixed reference point for the spectral WCS.  This is used with
the pd and ps imaging distortion functions.  Finally (d(l,b),s(l,b)) <->
(x,y) is obtained based on the reference pixel tied to a celestial
coordinate.


Issues and Thoughts

Slit spectroscopy can include a curved slit.  Therefore, one cannot just
use some standard form such as described for the distorted tan projection.

Due to the optical distortions from the additional dispersing elements
I don't believe it is feasible to define a standard 2D distortion function
such as is done for imaging systems.

Since fibers are frequently used to map the celestial field at a set of
points to a linear arrangement at the spectrograph entrance, there will
generally not be a single WCS for the frame.  Multiple WCS need to be used
to descibe the 2D format and, since instruments may include thousands of
fibers, a binary table using the Greenback conventions will be needed.

Dispersion calibrated but not resampled images, meaning a WCS with
spatial variations, are interesting for direct scientific use for
long-slit, Fabry-Perot, and objective prism data.  A WCS or
multiple WCS for multiobject, IFU, and cross-dispersed echelle data 
are useful for providing a starting point for the extraction of
spectra and conversion to other forms.  For example constructing
data cubes from IFUs and a long merged spectrum or a collection of
1D spectral orders from multi-order echelle data.

Data cubes from integral field units are "constructed" from a set of
extracted 1D spectra.  In this process the data cube WCS can be made
simple.  In particular, the spatial variation of the dispersion axis
with the spatial axes would be eliminated.

I am currently very interested in case 3.  A multi-WCS description would be
used to define the "apertures" (regions in (d,s) space) for extraction.
Initially the mask or IFU software that generates the focal plane would
provide an approximate WCS for each spectrum where the details of fd(x')
and fs('x') would be absent.  The absolute origins (the pixel reference
points) would uncertain though the relative origins would be well known.
Software would use the WCS to register to the data, measure fs('x'), and
record this improved WCS as the aperture definitions for automatic
extraction.  The initial dispersion information would be use to automate
the dispersion calibration and that could then be fed back to define the 2D
WCS for the raw exposure as well as leading to extraction to 1D spectra.
A data cube could also be built automatically for IFU instruments.

The thing limiting this from doing it entirely within a FITS WCS is
the representation of fs('x'), called the "trace" function in the
commonly used IRAF software.  As noted before, a cubic spline is best
suited for this and real data does show forms which are not simple
polynomials.


Is WCS-I a sufficient basis?

I am not completely sure.  It probably can be forced to work but maybe
there is room for a new concept.  The thing that nags at me is that
optical spectroscopic images are really degenerate in many cases.  In
other words you have an x/y 2D image that really maps to two spatial
axes on the sky and wavelength.  But then you have some axis
projection/reduction function that constrains one of the spatial
coordinates to be related to the other.  (This ignores the fact that
even with an aperture there is a multiplexing of the dispersion and
spatial coordinate that lies along the dispersion.) The concept of axes
reduction is not really in WCS-I other than through what some of us
consider is the kludge of defining dummy additional axes.  What would
be nice is to define the 3D WCS and some dimensional reduction which
would not necessarily be purely along one of the celestial axes or even
a line.


Date: Thu, 26 Apr 2001 12:04:16 -0400
From: Don Wells <dwells@cv.nrao.edu>
To: Frank Valdes <valdes@noao.edu>
Cc: dwells@NRAO.EDU, tody@tucana.tuc.noao.edu, davis@tucana.tuc.noao.edu
Subject: Re: WCS issues

Dear Frank, (and Lindsay and Doug)

I am really pleased to receive your fine summary of the technicalities
of optical spectroscopy.  IMO much of it should be incorporated into
WCS-III, because it is analogous to the detailed tutorial stuff that
is in WCS-II.

Frank Valdes writes:
 > .. some high energy instruments..

This reference surprised me, so I went looking, and immediately found
out about XMM-Newton. The relevant URL for grating spectroscopy is:

    http://xmm.vilspa.esa.es/user/uhb/xmm_uhb.html

>From the discussion there I infer that their data fits the model we
are discussing, except that they say that in their 'high time
resolution mode' they can do time-resolved spectro-photometry of
individual lines. Of course this is just adding a time axis to the WCS
header, so the model still works. I conclude that WCS-III should
include a discussion of the case of XMM-Newton, and that our
counterparts in that project should review WCS-III.

 > Of course one should be able to combine multiple transformations into
 > a single one.  However, I think it is very useful to understanding and
 > description to keep the terms somewhat separate.

Your multiple transformation description is a good thing for tutorial
purposes, and obviously is also a good way to build the software which
determines the calibrations.  However, for interchange purposes I hope
that it will be sufficient to combine the transformations.  The
community appears to think that simple polynomials should be used for
the celestial case; therefore I hope that they will be sufficient for
the spectral case too. (As always, my real goal is to get everybody to
agree on _some_ interchange standard, and I am prepared to adopt
almost anything that is technically adequate and minimizes the
negotiating effort.)

 > ..  Functions fd(x'), fd(y'), fs(x'), and fs(y') are 1D functions for
 > the dispersion and spatial axes.  In this case the functions can include
 > cubic splines which have been found to be very useful for spectral
 > applications..

Is the pixel correction function of WCS-I sufficient to support
cases where small-amplitude fine-scale corrections are needed?

 > ..The function ps(x',y') is what is sometimes refered to as the
 > "S-distortion"..

Last Friday I proposed adoption for spectral purposes of the 2-D
polynomial form we have adopted for the TAN function in WCS-II.  It
does not appear to me to be strong enough to represent the
delta-theta-as-function-of-radius form of S-distortion seen in
magnetically focussed image tubes. I did not raise a question about
this during the WCS-II negotiations because I am not aware of any
cameras of that type still in production use today (I last used them
about 25 years ago). Is this the type of distortion you are
mentioning, and are you aware of any current devices that need support
for it?

 > Case 2: Multiobject and Integral Field Unit Spectra
 > ..The data consist of multiple spectra having a small spatial extent..  

I assume that by IFU you mean data such as are shown in this URL: 

    http://www.ing.iac.es/PR/newsletter/news3/integral.html

I wonder whether we really need a FITS interchange solution for the
raw 2-D data in this case.  Isn't it more likely that we would
interchange data files like the gridded images shown in the INTEGRAL
web page?  Such data can be described well with the existing
WCS-I/-II/-III notations.

By 'multiobject' I assume that you mean conventional fibre
spectrographs such as the one used by SDSS:

http://www.apo.nmsu.edu/Telescopes/SDSS/eng.papers/19940311_FiberFeed/19940311.html

Again, I wonder whether the astronomy community really need a FITS
interchange solution for the raw 2-D data in this case.

If we could agree that there is no need to interchange single WCS
descriptions for these cases we would greatly simplify negotiations
for WCS-III.  However, if we need the capability, then probably
BINTABLE will rescue us, as we discuss below.

 > Case 3: Cross-dispersed Echelle Spectra
 > Cross-dispersed echelle spectra, where there are many orders and a small
 > slit can be treated (and are treated by IRAF) in exactly the same way as
 > the multiobject case..

Yes.

 > Case 5: Aperture-less Spectroscopy
 > 
 > ..Applications would need to interpret the WCS for some celestial
 > coordinate..

or for some spectral coordinate, e.g. imaging an emission nebula in
various lines!

Aperture-less spectroscopy is an inherently ambiguous business,
because radiation from two independent points on the sky can be imaged
to the same spectral pixels.  I think that the current WCS-III
dismisses the problem as completely insoluble, but probably it should
instead have a tutorial discussion of recommended practices for certain
particular cases.

 > Slit spectroscopy can include a curved slit.  Therefore, one cannot just
 > use some standard form such as described for the distorted tan projection.

Consider my proposal of Friday 04-20, in which there are 3 axes in the
header, 2 celestial and 1 spectral. Each of these axes would have a
40-term polynomial in the header.  I suspect (but have not yet proven)
that the celestial polynomials could be used with the WCS-II formalism
to map the spatial axis of the 2-D image onto a curved line on the
sky, and that the spectral polynomial (whose two axes are functions of
celestial coordinates) could then specify the spectral zero point of
each spectral row of the 2-D image. Do you think that this notation
would be sufficient to interchange curved-long-slit spectra?

 > Since fibers are frequently used to map the celestial field at a set of
 > points to a linear arrangement at the spectrograph entrance, there will
 > generally not be a single WCS for the frame.  Multiple WCS need to be used
 > to descibe the 2D format and, since instruments may include thousands of
 > fibers, a binary table using the Greenback conventions will be needed.

I agree. We are suggesting that an auxiliary BINTABLE with WCS keyword
values in its columns should be used for all situations in which an
n-D matrix contains a large number (n>25) of spectra and we want a WCS
notation for the spectra.  Indeed, if our community agrees that we
need to interchange WCS for _any_ cases with n>25, I would argue that
we should interchange _all_ cases of n>1 using auxiliary BINTABLEs.
This would produce a single scalable solution for the whole class of
problems: cross-dispersed echelle cases, multiple-slit cases and
multi-object fiber cases.  Indeed, even raw INTEGRAL data would work
with such BINTABLEs too.  Such a convention would support the
interchange and archiving of a wide variety of raw n-D spectroscopic
data. 

 > Dispersion calibrated but not resampled images, meaning a WCS with
 > spatial variations, are interesting for direct scientific use for
 > long-slit, Fabry-Perot, and objective prism data.

Yes.

 >  A WCS or
 > multiple WCS for multiobject, IFU, and cross-dispersed echelle data 
 > are useful for providing a starting point for the extraction of
 > spectra and conversion to other forms.  For example constructing
 > data cubes from IFUs and a long merged spectrum or a collection of
 > 1D spectral orders from multi-order echelle data.

Yes. 

 > Data cubes from integral field units are "constructed" from a set of
 > extracted 1D spectra.  In this process the data cube WCS can be made
 > simple.  In particular, the spatial variation of the dispersion axis
 > with the spatial axes would be eliminated.

This is exactly what happens in 'on-the-fly' spectral imaging with
single-dish radio telescopes: individual 1-D spectra sample the sky
irregularly and the data are gridded to a 3-D cube with a simple WCS.

 > I am currently very interested in case 3.  A multi-WCS description
 > would be used to define the "apertures" (regions in (d,s) space)
 > for extraction.

Probably you need the BINTABLE convention for this.  I suggest that you
propose the convention to the FITS community.

 > The thing limiting this from doing it entirely within a FITS WCS is
 > the representation of fs('x'), called the "trace" function in the
 > commonly used IRAF software.  As noted before, a cubic spline is best
 > suited for this and real data does show forms which are not simple
 > polynomials.

Can the pixel correction function defined in WCS-I be used for the
small-amplitude fine-scale corrections?

 > Is WCS-I a sufficient basis?
 > 
 > I am not completely sure.  It probably can be forced to work but maybe
 > there is room for a new concept..

Well, I confess that I am not completely sure either!  

My polynomial proposal of last Friday has the virtue that it closely
resembles a proposal which the community has already laboriously
negotiated and accepted. I therefore expect that it would require only
minimal further negotiating effort to win acceptance of it for the
spectroscopic case as welln.  Therefore, if it will work from a
technical point of view (aside from esthetic considerations) and will
not limit us in any data interchange cases that we are able to
anticipate, we should adopt it.

 > .. with an aperture there is a multiplexing of the dispersion
 > and spatial coordinate that lies along the dispersion..

The extreme case is objective prism imaging of galaxies. For stars you
can at least try to deconvolve the data, but for galaxies the signals
are hopelessly scrambled in most cases.

 > .. The concept of axes reduction is not really in WCS-I other than
 > through what some of us consider is the kludge of defining dummy
 > additional axes.  What would be nice is to define the 3D WCS and
 > some dimensional reduction which would not necessarily be purely
 > along one of the celestial axes or even a line..

It appears to me that the case of the straight slit at some PA is
completely covered by the dummy axis technique and CDi_j matrix. It
also appears to me that the TAN+poly representation _may_ be
sufficient to cover curved slits as well, as I suggested above.  If
technical considerations imply that we really must negotiate a
dimensional reduction description, then we should do it.  If, however,
the notational machinery we have already negotiated is _sufficient_
for interchange purposes (even if may be esthetically unpleasant for
some of us) then we should settle for the bird that is already in our
hands.

				 -=-

During the past ten days I have begun to get enquiries from members of
the FITS community regarding the state of WCS matters. I have told
these people that I am negotiating privately.  Your emails have not
been shown or forwarded to anyone.  My bet is that I can continue in
this manner for a while longer so that we can complete our discussion.

At this moment you are the only person known to me who is worrying
about FITS interchange notation for all of the spectroscopic cases.
We cannot freeze WCS-I until this analysis is done.  Thank you for
your efforts.

Regards,
Don



Date: Thu, 26 Apr 2001 17:34:05 -0700 (MST)
From: Frank Valdes <valdes>
To: dwells@cv.nrao.edu
Cc: tody, davis
Subject: Re: WCS issues

Hi Don,


Thanks for your comments and encouragement.

On the issue of dimensional reduction I agree that dummy axes are adequate
and maybe the best way to handle this.  The basic philosophical problem is
that NDIM is being used for both the WCS dimensionality and the raster
dimensionality.  If we just add WCSDIM and define the dimensions above NDIM
as being pixel 1 in the WCS that would allow a 2D raster to remain 2D while
allowing other axes for degenerate spatial, time, wavelength, etc.  WCSDIM
is used in IRAF but there is also a more general mapping as to which axes
are degenerate and what the default pixel value is for the missing axes.  I
would not suggest this extra stuff be adopt.  But just adding a single new
WCSDIM keyword  as something so simple that I don't see why anyone would
object and those who want to keep NDIM=WCSDIM are free to do so.

The issue that I am most adament about is the need for flexibility in
functional representations.  Having seen a variety of optical
spectroscopic data I have found that spline functions are required
to describe some coordinate behavior.  I don't believe it is a matter
of taste but that polynomials are simply not adequate.

Below is an example of things I have seen.  This shows the spatial offset
of the profile of a star or fiber as a function of the dispersion
coordinate.  The position is determined from centeroiding of many points at
each dispersion element and the errors can be as small as the points.  The
fraction of a pixel behavior is real and including it is necessary to
eliminate systematic effects in spectral extraction.  What you sometimes
see is fairly sharp kinks.  Using a polynomial for this requires a very
high order to follow the kink with the side-effect of residual wiggles at
other points.




         *                                  . 2
        * *                                 .
      **   *                                .
    **      *                               .
   *         ***                   *        .
  *             ***               * *       .   Deviation of object profile
 *                 *            **   *      .   in cross dispersion
*                   ****     ***     *      .
                        *****        *      .
                                      *     .
                                       *    .-2

	    dispersion axis


I have seen similar patterns in dispersion functions, the dispersion
correction as a function of pixel position on the dispersion axis.

> Can the pixel correction function defined in WCS-I be used for the
> small-amplitude fine-scale corrections?

Of course this could do it since it is totally general.  However, having
to go to another extension is something I am loathe to do.  If this
was the only option I would probably invent something to represent
the cubic spline fits as has been done already with IRAF spectral
files.  I don't want this to sound like a threat but just my reiteration
that I really want some more flexibility in functional descriptions
and find that it is needed in situations that I have encountered.

>> Of course one should be able to combine multiple transformations into
>> a single one.  However, I think it is very useful to understanding and
>> description to keep the terms somewhat separate.

> Your multiple transformation description is a good thing for tutorial
> purposes, and obviously is also a good way to build the software which
> determines the calibrations.  However, for interchange purposes I hope
> that it will be sufficient to combine the transformations.  The
> community appears to think that simple polynomials should be used for
> the celestial case; therefore I hope that they will be sufficient for
> the spectral case too. (As always, my real goal is to get everybody to
> agree on _some_ interchange standard, and I am prepared to adopt
> almost anything that is technically adequate and minimizes the
> negotiating effort.)

One problem with combining 2D and 1D functions into a single transformation
is the limited coefficient name space (ignoring the need for spline or
interpolator type of functions).  For example to get a simple 7th order
polynomial of x from eq 25 of WCS-II gives PV1_0, PV1_1, PV1_4,
PV1_7, ..., PV1_31.

> Last Friday I proposed adoption for spectral purposes of the 2-D
> polynomial form we have adopted for the TAN function in WCS-II.  It
> does not appear to me to be strong enough to represent the
> delta-theta-as-function-of-radius form of S-distortion seen in
> magnetically focussed image tubes. I did not raise a question about
> this during the WCS-II negotiations because I am not aware of any
> cameras of that type still in production use today (I last used them
> about 25 years ago). Is this the type of distortion you are
> mentioning, and are you aware of any current devices that need support
> for it?

I know that you like to think in terms of physical mechanisms.  And that is
a good way to attack a problem.  But sometimes it is not possible to
explain all the details in terms of some general physical mechanism.  The
term S-distortion is sometimes used to mean any spatial distortion and not
just image-tube distortions.  By S-distortion I mean exactly what is shown
in the figure above; the departures from a straight line of the path of a
spectrum through an image.  Some parts of the curvature are due to
"tangent" types of effects, some are due to corrector optics, and some are
due to who knows what.  As I said, the figure is not unusual for what you see
when you actually measure the position of a spectrum as it moves through an
image.

I agree that currently image tube intensifiers are not in favor.  Though
consider MAMA detectors which have the same intensification concept.  And
don't some of the HST spectrographs have some magnetic field effects?



IFS (integral field spectroscopy) and IFUs (integral field units) are a hot
area of OIR spectrographs.  It basically means sampling a nearly
continguous region of the sky with a set of apertures.  There are a variety
of types.  One is a bundle of fibers with the fiber ends being the
apertures.  A variant is to put a lenslet array (generally hexagonal or
square apertures) in front of the fibers to minimize the deadspace between
the fiber tips.  There is also just a pure lenslet array which produces a
complex image sort of like an objective prism.  Another important type is
image slicers.  By tilting thin mirrors a region of the sky can be mapped
to a set of neighboring "long-slit" images on a detector.  Having a WCS
attached to the slits or aperture spectra is part of the goal of then being
able to reconstruct a data cube solely from the information in the raw
image or auxilary tables.

Multiobject does mean independently target fibers as well as slit masks
where multiple slits are placed on different objects.


I think we are in agreement that some types of formats are not meant
to be interchanged in the sense of providing the final scientific coordinates
of position, wavelength/energy, and flux.  However what I am trying to
develop is an interchange that allows 2D multispectral images to be described
to some fairly good level such that software and archival users would
be able to say this pixel belongs to this spectrum which came from this
place in the sky and has this approximate wavelength.  I think a FITS
WCS is a good way to do this.

As far as multiple spectra, if you can define a WCS for a single spectrum
then you can have multiple WCS by using a table of CRVAL/CRPIX/etc.
** So if you can describe a single long-slit image then you've also
solved the problem for multislits, fibers, echelle orders, etc. **

I have already done a good job of it for describing things with linear
terms.  This is good enough for software to calibrate a raw image.
However, I was also hoping to use the identical format to record the result
after calibration with adjustment of the linear terms and addition of the
distortion corrections.  The only aspect for which I don't have a
reasonable answer yet is the "s-distortion" or "trace" for how a spectrum
wiggles in the cross dispersion direction.  If I had a 1D spline for the
spatial axis this would be do it.

> Probably you need the BINTABLE convention for this.  I suggest that you
> propose the convention to the FITS community.

I don't really see the need for a new convention.  It is simply the
"Greenbank" convention with columns with names such as CRPIX1, CRVAL1,
PV2_3, etc.  Each row then describes an independent WCS for the image
and provides the "aperture" description.  I also have a document that
proposes a keyword equivalent up to some number of objects.  Examples
are CRPIX1=CRP1nnnn, CRVAL1=CRV1nnnn.  The set allows up to 9999.
However, the general reaction is that if you get to n>100 this would
probably make for too long a header and a BINTABLE would be better.
I want to allow for both.  I am not trying to define a new thing but
just a keyword convention for certain types of data and the keyword
dictionary would be published.  The draft is what was advertised at
ADASS and is in the Pence's link of dictionaries.

> I agree. We are suggesting that an auxiliary BINTABLE with WCS keyword
> values in its columns should be used for all situations in which an
> n-D matrix contains a large number (n>25) of spectra and we want a WCS
> notation for the spectra.  Indeed, if our community agrees that we
> need to interchange WCS for _any_ cases with n>25, I would argue that
> we should interchange _all_ cases of n>1 using auxiliary BINTABLEs.
> This would produce a single scalable solution for the whole class of
> problems: cross-dispersed echelle cases, multiple-slit cases and
> multi-object fiber cases.  Indeed, even raw INTEGRAL data would work
> with such BINTABLEs too.  Such a convention would support the
> interchange and archiving of a wide variety of raw n-D spectroscopic
> data. 

Yes this is what I am suggesting.  The distinction that needs to be made is
that multiple WCS is intended for different views of the data such mm on
the focal plane, pixels in a larger mosaic, and traditional astronomical
coordinates such as ra/dec.  This is different than describing multiple
instances of the same view with different reference points which is
what I am suggesting.


One thing that I find missing for describing multispectral data that I
would just float for your comment is that there is no concept of a region
of validity for a WCS.  For the multispectral case, including echelle
orders, there would be coordinates giving the limits of the
cross-dispersion and dispersion axes for which the WCS, and so the
spectrum, apply.  For example if the spatial axis is arcseconds from the
center of the aperture and the spectrum is bandpass limited on the image by
filters and detector response you would say a spectrum occupies -2 to 2
arcsec and 4000 to 8000 A.  The WCS could be inverted to say the spectrum
occupies some region in the image.  I have invented such keywords in my
design (CMINi/CMAXi).

The idea that some spectral formats have limited dispersion so that
multiple spectra may appear along the same spatial coordinate (e.g.
dispersion is along lines and there are multiple non-overlapping spectra at
the same line) is used with pure lenslet IFUs, slitlets and objective
prisms.  There is slitlet data where there are three tiers of spectra.
Again, this complex pattern is what should be describably by the slit
making software as it optimizes the packing of the sky on the 2D detector
using a FITS WCS.  Internally such software already has a model of how the
objects with dispersion will fall onto the detector in order to make the
slit mask.  For objective prism data a catalog of object positions on the
sky along with knowledge of the plate scales and dispersion would allow a
WCS description for each object on the image.


>> Slit spectroscopy can include a curved slit.  Therefore, one cannot just
>> use some standard form such as described for the distorted tan projection.

> Consider my proposal of Friday 04-20, in which there are 3 axes in the
> header, 2 celestial and 1 spectral. Each of these axes would have a
> 40-term polynomial in the header.  I suspect (but have not yet proven)
> that the celestial polynomials could be used with the WCS-II formalism
> to map the spatial axis of the 2-D image onto a curved line on the
> sky, and that the spectral polynomial (whose two axes are functions of
> celestial coordinates) could then specify the spectral zero point of
> each spectral row of the 2-D image. Do you think that this notation
> would be sufficient to interchange curved-long-slit spectra?

Whether or not some standard  polynomial of fixed number of terms and form
is sufficient depends on the slit shape.  Simple curves such as a quadratic
or piece of a circle would work.  Special purpose curved slits can have a
niche such as with gravitationally lensed arcs.  But the general case is
only for the sake of argument.  ** I don't really believe a standard WCS
needs to be able to describe funny shaped spectral apertures. **  But if it
is possible then we should consider it.  Again, providing some general 1D
polynomial would address the curved slit situation along with using a
degenerate axis.

I have not actually thought in detail about whether one could  set up a
degenerate WCS so that instead of just arcseconds along the slit one could
actually get ra/dec.  This would be important for using multiple (source)
WCS in reconstructing a data cube from IFU data, especially image slicers.
I was currently only considering that each aperture would have non-WCS
keywords giving the position, dimensions, and orientation of the apertures
on the sky.

I guess this would still be necessary at least for the case of hexagonal
lenses.  After extraction from the 2D format you have a series of 1D
samples.  People who want to reconstruct a data cube then sample the 1D
fluxes into 2D hex pixels to make a spatial representation of the
observation.  But the aperture center could still be transmitted through
CRVALn information rather than APRAnnnn/APDEnnnn keywords.

				-=-

I think just having you be aware at this point and commenting on what
we send you as a knowledgable independent view point is good and I
appreciate it.


Regards,
Frank


Date: Mon, 30 Apr 2001 16:51:48 -0700 (MST)
From: Frank Valdes <valdes>
To: dwells@cv.nrao.edu, tody, davis
Subject: Re: WCS issues

Hi Don,

I have continued to think about the spectral WCS issues.  The comment I
made earlier that spectra are always inherently three dimensional led
me to realize that the path to dealing with the spatial variation of
spectral data lies is defining classes of 3D WCS.  Dimensional
reduction to 2D and 1D are then used for things like long slit and
extracted spectra.  This has the advantage of providing the full 2D
spatial information in these reduced situations.   One can  define
various functional forms for the 3D WCS as is done with the 2D
celestial WCS and the 1D dispersion WCS.

As far as functional forms I only propose a general 3D polynomial.
I believe this would be sufficient for most spectra.  I am willing
to give up on the splines and deal with deviations from what can be
produced with polynomials either with a pixel regularization for
interchange or non-WCS keywords for IRAF specific software such as
to record the trace function for extraction.  I don't think we
actually need to describe the observational formats to absolute perfection
and concerns with astrometric accuracy are not so important with
data that is primarily taken for the spectral content.

Wide-field astrometric issues would only be important for Fabry-Perot and
possibly objective prisms.  The path for Fabry-Perot would be to define a 3D
WCS (say CTYPEs of WAVE-WAV-FP, RA---TAN-FP, DEC--TAN-FP.  with the
specific FP physical behavior and built upon the celestial WCS definitions
including a distortion function.  Because it is "-TAN-FP" instead of "-TAN"
the form of the function can be defined differently if needed.

For aperture spectroscopy the placement of the aperture takes care of the
astrometry and the spatial part of the WCS would be local to the position
of the aperture.  IFS will also probably never be wide-field and, as noted
before, such data needs to be reconstructed into a collection of 1D spectra
or a regularized data cube.

The parts of my earlier thoughts based on dealing with multiple spectra
(or cross-dispersion echelle orders) as multiple WCS is independent
of the basic concept of how to deal with a single spectrum as discussed
here.

So, apart from defining new CTYPE classes for 3D WCS which are independent
of the other papers, the only issue that I want to ask be considered
is the WCS Dimensionality convention given in the appendix below.

Lindsey and Doug are seeing these proposals at the same time you are
so they may have a different view of things.  However, I know Lindsey
has also suggested and agrees with the WCS Dimensionality proposal given
here.

What do you think?

Regards,
Frank




Astronomical spectra are inherently three dimensional, consisting of two
celestial spatial axes and a dispersion axis.  One dimensional and two
dimensional spectra are formed by dimensional reduction of the spatial
axes.  Therefore the natural WCS for spectra is three dimensional.  Two
dimensional imaging can also be thought of as a 3D WCS with the dispersion
axis reduced, but for broadband observations astronomical practice is to
identify the filter rather than use a dispersion coordinate.  Of course,
one can add additional parameters by using higher dimensions but I can't
think of realistic cases where the coupling of the dimensions needs to be
addressed.  This discussion is limit to 3D spectral/spatial WCS.

When the spatial and dispersion axes are independent then the concepts and
definitions of FITS WCS described in WCS I, II, and III can be applied.
These define classes of 2D spatial systems and 1D dispersion systems.  In
this paper we define 3D WCS where the spatial and dispersion coordinates
are not independent.

We follow the general concepts for FITS WCS laid out in WCS-I.  This
provides for defining a variety of WCS identified by CTYPEn keywords.  To
start this discussion of 3D WCS let us consider a system based only on
three dimensional polynomials.  Let us call this class "P3D" and build
CTYPEn values using the prefixes from WCS II and III, namely "RA--",
"DEC-", "xLON", "xLAT", "WAVE", "FREQ", "VELO", "VRAD", "VOPT", and
"ZOPT".  We allow the spatial projections and dispersion gridding functions
that do not depend on PV coefficients.  For example, a common optical WCS
would have axes "RA---TAN-P3D", "DEC--TAN-P3D", and "WAVE-WAV-P3D" where
only the pure TAN projections without distortion terms would be used.

The standard linear transformation maps the image coordinates p_l to the
intermediate world coordinates x_l, where x_1 is the coordinate associated
with CTYPE1, x_2 with CTYPE2, and x_3 with CTYPE3.  The higher order
3D polynomial corrections are given by

    x_l' = x_l + sum_ijk {PVl_ms x_1^i x_2^j x_3^k}  i=0-I, j=0-J, and k=0-K

The values of I, J, and K are given in the parameters PV0, PV1, and PV2.
The values of m are ordered so that i varies fastest and k varies slowest
over their ranges so that m = 3 + k*(J+1)*(I+1) + j*(I+1) + i.  The polynomial
is defined so that the default values for the PV terms all default to
zero.

After the polynomial corrections the final world coordinates are
computed by applying the reference values W_l and the indicated
projection.


WCS Dimensionality

The key to the generality of the 3D WCS for spectroscopic data is the
ability to use a WCS with a higher dimensionality than the natural sampling
of data.  Currently this is accomplished by defining the image to be of the
same dimensionality as the WCS and setting dummy axes of length one for the
extra WCS axes.  The appendix describes an alternate convention.

When the dimensionality of the image is increased to match the WCS it is
highly recommended that the dummy axes be added after the natural sampling
axes and that the coordinate reference pixel values, the CRPIX, be given a
value of one.  With the convention described in the appendix the extra WCS
dimensions are required to be for pixel axes beyond the dimension of the
image.  There is no requirement on CRPIX but the same recommendation
applies.

The CD transformation rotates the three image axes to align with the three
world coordinate axes apart from the distortions handled by the non-linear
corrections and projections.  Note that the dispersion axis need not
be orthogonal to the spatial plane which may be useful in some situations.


Cases with Reduced Dimensionality

To simplify this discussion, we assume the first image and world axis is
the dispersion and that they are already aligned.  This means there are no
rotation terms between the first image axis and the spatial axes so that
the CD1_2, CD1_3, CD2_1, and CD3_1 keywords are zero or absent.  However,
the alignment of the spatial axes with the image axes is not constrained.

The spatial axes may also include a projection function.  Also for
simplicity we choose the CAR spatial projection.

1. Long Slit

Consider a rectangular slit.  The center of the slit is placed at some
position on the sky and oriented at some position angle, defined as the
direction of the long dimension of the slit measured from north and
increasing to the east.  The image of the slit is dispersed along the first
image axis of a 2D detector.

In the first example the slit is perfectly straight, the dispersion is
perfectly linear in wavelength, and there are no distortions.
The WCS at the center of the reference pixel is 5400A with right ascension
of 1h 23m 45s and declination of 1d 23m 45s.  The pixel scale is 4.3 A/pixel
and 1 arcsec/pixel.  The position angle is 30 degrees.

NDIM = 3
NAXIS1 = 1024
NAXIS2 = 1024
NAXIS3 = 1
CTYPE1 = WAVE-WAV-P3D
CTYPE2 = RA---CAR-P3D
CTYPE3 = DEC--CAR-P3D
CRPIX1 = 512
CRPIX2 = 512
CRPIX3 = 1
CRVAL1 = 5400.
CRVAL2 = 20.9375
CRVAL3 = 1.3958333333333
CD1_1  = 4.3
CD2_2  = 1.3888878249646E-4
CD2_3  = 2.4056267358801E-4
CD3_2  = -2.4056267358801E-4
CD3_3  = 1.3888878249646E-4
CUNIT1 = Angstrom
CUNIT2 = deg
CUNIT3 = deg

This is shown as a 3D image with a third axis length of 1.  Under the
WCS Dimensionality proposal NDIM would be 2 and NAXIS3 would be absent.

Since there are no polynomial coefficients the coordinates are
evaluated by the standard linear transformation and coordinate
reference values as given by

W   = CRVAL1 + CD1_1 * (p_1 - CRPIX1)
RA  = CRVAL2 + CD2_2 * (p_2 - CRPIX2)
DEC = CRVAL3 + CD3_2 * (p_2 - CRPIX2)

There are no p_3 terms because p_3 takes the value of one for the missing
image axis and CRPIX3 is set to 1.  The CD2_2 and CD3_2 terms are computed
by

CD2_2 = +/- S * sin(pa)
CD3_2 = +/- S * cos(pa)

where S is the spatial pixel scale (given in degrees/pixel) and the signs
depend on the optical projections of the celestial coordinates on
to the detector.  The other elements of the CD matrix need to be
defined to assure a valid matrix even though the dimensional reductions
eliminates use of CD2_3 and CD3_3.

In this simple example there are no dependencies between the dispersion
and spatial axes.  In this case one could also simply use a 1D WCS
for the dispersion axis and a 2D celestial WCS for the spatial axes.

A more realistic case is when the distortions in the spectrograph and
camera cause lines of constant spatial position to vary along the
dispersion axis in the image.  Also the slit may not perfectly linear, or
even curved by design, so that lines of constant wavelength will vary along
the spatial axis.  Finally, dispersing elements are typically not linear
with pixel.  Below is a representation of a 2D long slit image showing the
peak of a bright source (*) in the slit and some night sky lines (.).
X marks the reference point.


+----------------------------------------+
+           .              .             +
+       **********           .           +
+   ****    .     ****        .        **+
****       .          X***     .   ****  + Position
+        .                *********      + along the slit
+      .                       .         +
+   .                         .          +
+ .                         .            +
+----------------------------------------+
                Dispersion


The WCS might be given by

NDIM   = 3
NAXIS1 = 1024
NAXIS2 = 1024
NAXIS2 = 1
CTYPE1 = WAVE-WAV-P3D
CTYPE2 = RA---CAR-P3D
CTYPE3 = DEC--CAR-P3D
CRPIX1 = 512
CRPIX2 = 512
CRPIX3 = 1
CRVAL1 = 5400.
CRVAL2 = 20.9375
CRVAL3 = 1.3958333333333
CD1_1  = 4.3
CD2_2  = 1.3888878249646E-4
CD2_3  = 2.4056267358801E-4
CD3_2  = -2.4056267358801E-4
CD3_3  = 1.3888878249646E-4
CUNIT1 = Angstrom
CUNIT2 = deg
CUNIT3 = deg
PV1_0  = 7
PV1_1  = 2
PV2_0  = 3
PV1_5  = ...
PV1_6  = ...
PV1_7  = ...
PV1_8  = ...
PV1_9  = ...
PV1_10 = ...
PV1_19 = ...
PV2_5  = ...
PV2_6  = ...

The polynomial terms show only the non-zero values.  In this example
the dispersion is described by variation in wavelength with position
along the dispersion of 7th order and the zero point of the dispersion
varies as a quadratic with position along the slit.  The spatial
distortion is a cubic function of the position along the first image
dimension with no other dependence.  The functional forms are

x_1' = x_1 + PV1_5 * (x_1)^2 + ... PV1_10 * (x_1)^7 + PV1_19 * (x_2)^2

x_2' = x_2 + PV2_5 * (x_1)^2 + PV2_6 * (x_1)^3


2. One Dimensional Spectra

One dimensional spectra are typically extracted from 2D or 3D images.
While there are then no spatial axes, use of a 3D WCS will define the
celestial coordinate for the 1D spectrum.  A spectrum is generally extracted
along a line of constant spatial coordinate.  In this case there will
be no dependence of the spatial coordinate on position.  So one has the
option of using a 3D WCS or a 1D coordinate axis for the dispersion and
a 2D celestial WCS for the reduced spatial axes.

The example below is for an extraction of the 2D long slit example given
earlier.  The extraction followed the object peak line of the figure (line
of constant spatial position) and the non-linear wavelength dispersion
function was removed by resampling to a uniform wavelength sampling.
Summation across the dispersion was done along lines of constant wavelength
as schematically shown by the figuer.

NDIM = 3
NAXIS1 = 1024
NAXIS2 = 1
NAXIS3 = 1
CTYPE1 = WAVE-WAV-P3D
CTYPE2 = RA---CAR-P3D
CTYPE3 = DEC--CAR-P3D
CRPIX1 = 512
CRPIX2 = 1
CRPIX3 = 1
CRVAL1 = 5400.
CRVAL2 = 20.9375
CRVAL3 = 1.3958333333333
CD1_1  = 4.3
CD2_2  = 1.3888878249646E-4
CD2_3  = 2.4056267358801E-4
CD3_2  = -2.4056267358801E-4
CD3_3  = 1.3888878249646E-4
CUNIT1 = Angstrom
CUNIT2 = deg
CUNIT3 = deg

Note that under the WCS Dimensionality proposal NDIM could be 1 and
the NAXIS2 and NAXIS3 keywords would be absent.


3. Fibers

Optical fibers scramble the spatial information so that the spatial
axes cannot be related to celestial coordinates.  However, the center
of the aperture still has a celestial coordinate which the WCS can
provide.  The orientation of the celestial axes on the detector is
then irrelevant



APPENDIX: Proposal on WCS Dimensionality

Currently the dimensionality of the WCS is assumed to match the
dimensionality of the image raster as specified by the NDIM keyword.  If
the WCS is of a higher dimensionality than the image the current
recommendation is to add dummy image dimensions.  This is done by setting
NDIM to the higher dimensionality of the WCS and adding NAXISn keywords
with a value of 1 for the extra dimensions.

This proposal defines the WCS dimensionality to be the maximum of the image
dimensionality, as specified by NDIM, and the highest index of the CTYPEn
keywords.  Thus, if NDIM is 2 but there is a CTYPE3 keyword then the WCS
dimensionality would be 3.  To evaluate the WCS, pixel coordinates of one
would be used for the higher dimensions.

Software that deals with FITS files where the WCS dimensionality indicated
by the CTYPE keywords is higher than NDIM may chose to internally represent
the image as a raster with the higher dimensionality and axes lengths of 1
or treat the image raster with the dimensionality specified by NDIM but
evaluate coordinates with dummy pixel coordinates of 1 for the extra
dimensions of the WCS.

This proposal follows the basic FITS principle that existing FITS files
using the the extra raster dimensions of length 1 are still be valid.  New
files created using this proposal will require reading software to adopt
one of the easy to implement options described here.

In earlier data there was the possibility that no CTYPE keywords would be
present.  This proposal continues to support this possiblity.  However, in
order to define a WCS dimensionality higher than the image dimensionality
given by NDIM requires that CTYPE keywords be used (technically only the
one for the highest dimensionality).  Not including CTYPE keywords is a
practice which is discouraged.


Date: Mon, 30 Apr 2001 17:38:24 -0700 (MST)
From: Frank Valdes <valdes>
To: dwells@cv.nrao.edu, tody, davis
Subject: Re: WCS issues

Don,

Doug pointed out that it sounds like I am saying there is no need for
a WCS which is not polynomial.  To clarify,  what I am describing is
one representation, based on polynomials, to describe the
interdependence of spatial and dispersion axes to a reasonable level
of detail in raw 2D or 3D observational data.

When extracting a 1D spectrum in such a way that you avoid resampling,
such as is done by summing purely along rows or columns over a limited
spatial extent where you accept that any spatial variation might
slightly degrade the spectral resolution, you have a 1D spectrum which
is not dispersion resampled.  You do this to avoid correlations in the
noise and interpolation smoothing.  For this type of common situation
(this is how IRAF extracts spectra), astronomers do want a highly
accurate 1D WCS for dispersion where the data itself is not resampled.
For this situation, now the provence of the WCS III paper, the 1D
polynomial proposed for the non-linear corrections (section 2.2) is not
sufficient.  In this case a sampled system, such as splines or fully
sampled pixels is desirable.  I would still want to provide such a high
accuracy 1D representation for extracted 1D spectra.

Frank


From valdes Fri May 18 15:43:22 2001
Date: Fri, 18 May 2001 15:43:20 -0700 (MST)
From: Frank Valdes <valdes>
To: dwells@cv.nrao.edu, tody, davis
Subject: Re: WCS issue

Hi Don, Doug, and Lindsey,

I have been working hard on expanding and experimenting with some of the
ideas I mentioned earlier.  [Ideas concerning spectral WCS and coupling of
spectral and spatial WCS.  FV]  I now have some definite proposals.  I
believe these address the questions of how to deal with real spectral data;
in particular, how to deal with coupled spatial and spectral axes.  In the
process I am proposing a more general mechanism for distortion corrections
to all the defined coordinate system types.

In order to make sure the scheme works I have a prototype implementation
through the IRAF MWCS routines.  The examples included are all from real
FITS images (though not real data) using the keywords shown.

See

http://iraf.noao.edu/projects/spectroscopy/formats/pc.html
http://iraf.noao.edu/projects/spectroscopy/formats/spec3d.html

Frank


>From dwells@fits.cv.nrao.edu Wed Jun  6 08:09:14 2001
Date: Wed, 6 Jun 2001 11:09:11 -0400
From: Don Wells <dwells@cv.nrao.edu>
To: Frank Valdes <valdes>
Cc: dwells@NRAO.EDU, tody, davis
Subject: Re: WCS issue

Dear Frank, (and Doug and Lindsey)

Frank Valdes writes:
 > ..  I now have some definite proposals. [regarding] .. how to deal
 > with coupled spatial and spectral axes.. I am proposing a more
 > general mechanism for distortion corrections to all the defined
 > coordinate system types..

I have been wishing that somehow we could avoid making changes to the
WCS-I document, because any proposed changes will lead to more rounds
of discussion and negotiation.  However, I have hesitated to put WCS-I
to a vote precisely because I have been suspicious that the needs of
WCS-III will require us to make some changes to WCS-I.  Indeed, you
have confirmed my suspicions. Your general-purpose distortion
polynomial and implied-axis proposals would definitely require changes
to WCS-I, and probably to WCS-II also (because you propose
incorporating the functionality of the TAN polynomials into the place
currently occupied by the pixel regularization matrix).  So, I
conclude that we must have those rounds of discussion and negotiation.

Your ideas appear to me to offer a valid and general solution to the
whole class of instrumental distortion problems, and therefore they
should be considered seriously by our committees.

 > .. have .. prototype implementation .. IRAF MWCS routines..

Excellent!

 > http://iraf.noao.edu/projects/spectroscopy/formats/pc.html
 > http://iraf.noao.edu/projects/spectroscopy/formats/spec3d.html
 > .. your comments would be useful..

Nice work. Elegant web pages. My notes, questions and suggestions:

[In the current draft the suggested DCF name is PLYs.  Previously it
was PCn (where n is a number).  This is what is refered to in the
following comments.  Use of PCn leads to confusion about what is
the function name and what are the coefficient keywords. FV]

* The discussion in Section 2 of 'pc.html' regarding the proposed
  '-PC1' suffix to CTYPEis values confuses me.  It seems to say that
  the suffix can be of the form '-PCi', but all examples shown in
  'pc.html' use only '-PC1', and the examples leave me wondering what
  would be the purpose of a '-PC2' suffix.  Either some additional
  clarifying text/examples should be added, or the suffix should be
  changed, perhaps to just '-PC'.  Actually, one could ask what
  purpose the suffix serves--as far as I can see, the only information
  it conveys is that 'there _might_ be pixel correction terms
  associated with this axis'.  You have designed the PCi_ks terms to
  default to zero, and it is also true that it will not be a disaster
  if a reader algorithm ignores them, so maybe the suffix is not
  really needed?  The FITS community has not always been consistent in
  this matter---sometimes we have hard indications of what is present,
  and other times we leave it to reader algorithms to discover the
  features.  Like a lawyer, I could argue the case for either side in
  that debate.  If you favor continuing with a positive indication in
  your proposal, I recommend that you also consider an alternative
  notation in which a separate logical-valued keyword would signal
  that PCi_ks keywords may be present.

* The implicit NAXISi=1 scheme which you propose in 'spec3d.html' is a
  clever and elegant use of FITS notation; in my opinion it is quite
  consistent with the spirit of FITS.  I don't think it violates any
  existing rules.

* Your 'spec3d.html' text refers to the 'NDIM' keyword; I presume that
  'NDIM' is the IRAF name.  You should change the text to use 'NAXIS'.

* Figure 2 in your 'pc.html' text shows values for CTYPEi; I recommend
  that lowercase 's' be appended to 'CTYPEi' to remind the reader of
  the multiple WCS capability.  Actually, more generally, I recommend
  that 's' be appended to WCS keywords everywhere in your two
  documents.

* I recommend that all of your header examples include the NAXIS and
  NAXISi keywords and values.

* The last paragraph of Section 1 of 'pc.html' asserts that the pixel
  regularization matrix must be 'a fully populated image', which I
  take to mean that you are asserting 1-to-1 pixel association (full
  size).  This is not true: the scheme described in WCS-I.A includes
  an interpolation capability, so regularization matricies would be
  much smaller than images in practice.

Regards,
Don
-- 
  Donald C. Wells      Scientist - GBT Project        dwells@nrao.edu
                    http://www.cv.nrao.edu/~dwells
  National Radio Astronomy Observatory                +1-434-296-0277
  520 Edgemont Road,   Charlottesville, Virginia       22903-2475 USA



>From valdes Wed Jun  6 09:28:28 2001
Date: Wed, 6 Jun 2001 09:28:26 -0700 (MST)
From: Frank Valdes <valdes>
To: dwells@cv.nrao.edu
Subject: Re: WCS issue
Cc: tody, davis

Hi Don,

Thanks for looking at the proposals and the positive comments.

 > * The discussion in Section 2 of 'pc.html' regarding the proposed
 >   '-PC1' suffix to CTYPEis values confuses me.  It seems to say that
 >   the suffix can be of the form '-PCi', but all examples shown in
 >   'pc.html' use only '-PC1', and the examples leave me wondering what
 >   would be the purpose of a '-PC2' suffix.  Either some additional
 >   clarifying text/examples should be added, or the suffix should be
 >   changed, perhaps to just '-PC'.  Actually, one could ask what
 >   purpose the suffix serves--as far as I can see, the only information
 >   it conveys is that 'there _might_ be pixel correction terms
 >   associated with this axis'.  You have designed the PCi_ks terms to
 >   default to zero, and it is also true that it will not be a disaster
 >   if a reader algorithm ignores them, so maybe the suffix is not
 >   really needed?  The FITS community has not always been consistent in
 >   this matter---sometimes we have hard indications of what is present,
 >   and other times we leave it to reader algorithms to discover the
 >   features.  Like a lawyer, I could argue the case for either side in
 >   that debate.  If you favor continuing with a positive indication in
 >   your proposal, I recommend that you also consider an alternative
 >   notation in which a separate logical-valued keyword would signal
 >   that PCi_ks keywords may be present.

I will clarify this and give examples in the web doc.  The suffix is
important for the following reasons.  The suffix not only implies there are
likely to be coefficients (though as you say this could be done
differently) but also identifies which coordinates are coupled.  Since
there may be more than one group of coupled axes you need to have different
suffixes.  As a contrived example, suppose you have six world axes of RA,
DEC, WAVELENGTH, TIME, TEMPERATURE, and POLARIZATION.  Now suppose the
RA/DEC/WAVELENGTH are coupled as in my other examples but are independent
of the other axes.  Now also suppose TEMPERATURE and POLARIZATION are
coupled and that TIME is non-linear.  Then in my proposal you would have

RA---TAN-PC1
DEC--TAN-PC1
WAVE-WAV-PC1
TEMPERATURE-PC2
POLARIZATION-PC2
TIME-PC3

PC1 --> independent variables are intermediate coordinates from RA,
	DEC, and WAVE and cartisian radius of these
PC2 --> independent variables are temperature and polarization and radius
PC3 --> independent variable is time

The last case of TIME-PC3 just indicates a TIME coordinate type with a 1D
non-linear polynomial.  As you noted this last case would not really be
necessary if you said all axes may have PC coefficients but I think it is
clearer to have the indication.

Another reason for the suffix that I had in mind is to allow new functional
forms, possibly also using the PC keywords.  At one point I thought there
would be (and might still be) an FP suffix for Fabry-Perot.  It was only an
epiphany that I realized that in any realistic case a short polynomial
expansion gives the same result so FP need not require another functional
form.  Another example as a random thought is WAVE-WAV-TR1 might be a
trignometric polynomial.

In my early thinking I was going to suggest using PV for the functional
coefficients.  But then I realised that PVs are used not only in the
non-linear intermediate coordinate transformation (such as proposed in
WCS-III and the distorted tan in WCS-II) but that some of the ideal
projection models also used them.  So in order to have, for instance,
RA---AIR-PC1 you need the PC for the coupling/non-linear transformation and
the PV for the theta_b parameter.

I should note that Lindsey made a strong statement to me that she prefers
all the information about the transformation(s) be part of the CTYPE
keyword and not some combination such as another keyword that says "now
modify the CTYPE by applying this other transformation".  I would be
comfortable either way but I do see the logic of having all the information
in the CTYPE keyword as long as it does not become overloaded.  Now we are
suggesting three pieces of information encoded in CTYPE:
"coordinate_type-ideal_projection-intermediate_coupling_transformation".

 > * Your 'spec3d.html' text refers to the 'NDIM' keyword; I presume that
 >   'NDIM' is the IRAF name.  You should change the text to use 'NAXIS'.

My mistake, I did mean NAXIS.

 > * Figure 2 in your 'pc.html' text shows values for CTYPEi; I recommend
 >   that lowercase 's' be appended to 'CTYPEi' to remind the reader of
 >   the multiple WCS capability.  Actually, more generally, I recommend
 >   that 's' be appended to WCS keywords everywhere in your two
 >   documents.

There is always the conflict between showing the most common form and
always showing the general form.  I can take your recommentation.

 > * I recommend that all of your header examples include the NAXIS and
 >   NAXISi keywords and values.

OK

 > * The last paragraph of Section 1 of 'pc.html' asserts that the pixel
 >   regularization matrix must be 'a fully populated image', which I
 >   take to mean that you are asserting 1-to-1 pixel association (full
 >   size).  This is not true: the scheme described in WCS-I.A includes
 >   an interpolation capability, so regularization matricies would be
 >   much smaller than images in practice.

Sorry, there is so much to read that sometimes I skim and miss something.
There papers have lots of details.  You are right so my criticism
is simply the need to have an extra extension.

Regards,
Frank


>From dwells@fits.cv.nrao.edu Wed Jun  6 11:25:41 2001
Date: Wed, 6 Jun 2001 14:25:38 -0400
From: Don Wells <dwells@cv.nrao.edu>
To: Frank Valdes <valdes@noao.edu>
Cc: dwells@NRAO.EDU, tody@tucana.tuc.noao.edu, davis@tucana.tuc.noao.edu
Subject: Re: WCS issue

Frank, (Doug, Lindsey)

Frank Valdes writes:
 > 
 >  > * The discussion in Section 2 of 'pc.html' regarding the proposed
 >  >   '-PC1' suffix to CTYPEis values confuses me..
 > 
 > I will clarify.. The suffix.. not only implies there are likely to
 > be coefficients..but also identifies which coordinates are
 > coupled. Since there may be more than one group of coupled axes you
 > need to have different suffixes..

OK. I understand and agree.

 > .. Lindsey.. prefers all the information about the
 > transformation(s) be part of the CTYPE keyword and not some
 > combination..  we are suggesting three pieces of information
 > encoded in CTYPE:
 > "coordinate_type-ideal_projection-intermediate_coupling_transformation".

OK. Sounds reasonable to me.

Regards,
Don



>From valdes Wed Jun  6 11:50:47 2001
Date: Wed, 6 Jun 2001 11:50:45 -0700 (MST)
From: Frank Valdes <valdes>
To: dwells@cv.nrao.edu
Subject: Re: WCS issue
Cc: davis, tody

Don,

I have clarified the web pages based on your comments.  I now understand
the reason for your suggestion about including NAXIS keywords in the
examples.  I had left them out to imply that either the dummy axes
approach or the WCS dimensionality idea could be used.  But if we are
to propose that NAXIS may disagree with CTYPE the examples should
illustrate how that would work.

Frank

From dwells@fits.cv.nrao.edu Wed Jun 20 08:20:03 2001
Date: Wed, 20 Jun 2001 11:20:01 -0400
From: Don Wells <dwells@cv.nrao.edu>
To: Frank Valdes <valdes@noao.edu>
Cc: dwells@NRAO.EDU, tody@tucana.tuc.noao.edu, davis@tucana.tuc.noao.edu
Subject: Re: WCS issue

Dear Frank, (and Doug and Lindsey)

Frank Valdes writes:
 > I will clarify this and give examples in the web doc.  The suffix is
 > important for the following reasons.  The suffix not only implies there are
 > likely to be coefficients (though as you say this could be done
 > differently) but also identifies which coordinates are coupled.  Since
 > there may be more than one group of coupled axes you need to have different
 > suffixes..

Yes, I see that. 

.  As a contrived example, suppose you have six world axes of RA,
 > DEC, WAVELENGTH, TIME, TEMPERATURE, and POLARIZATION.  Now suppose the
 > RA/DEC/WAVELENGTH are coupled as in my other examples but are independent
 > of the other axes.  Now also suppose TEMPERATURE and POLARIZATION are
 > coupled and that TIME is non-linear.  Then in my proposal you would have
 > 
 > RA---TAN-PC1
 > DEC--TAN-PC1
 > WAVE-WAV-PC1
 > TEMPERATURE-PC2
 > POLARIZATION-PC2
 > TIME-PC3
 > 
 > PC1 --> independent variables are intermediate coordinates from RA,
 > 	DEC, and WAVE and cartisian radius of these
 > PC2 --> independent variables are temperature and polarization and radius
 > PC3 --> independent variable is time

Section 3 of your "Functional corrections and axis coupling.."
document says ".. keywords PC1_ns, where 1 is the axis number..".
Does this mean that the above example should be:

 > RA---TAN-PC1-PC2-PC3
 > DEC--TAN-PC1-PC2-PC3
 > WAVE-WAV-PC3
 > TEMPERATURE-PC4-PC5
 > POLARIZATION-PC4-PC5
 > TIME-PC6

with six polynomials

PC1_ns    2D poly corrects RA using RA, DEC
PC2_ns    2D poly corrects DEC using RA, DEC
PC3_ns    3D poly corrects WAVE using RA, DEC, WAVE
PC4_ns    2D poly corrects TEMPERATURE using TEMPERATURE, POLARIZATION
PC5_ns    2D poly corrects POLARIZATION using TEMPERATURE, POLARIZATION
PC6_ns    1D poly corrects TIME using TIME

 > I should note that Lindsey made a strong statement to me that she prefers
 > all the information about the transformation(s) be part of the CTYPE
 > keyword..

I think that the above CTYPE notation is comparatively easy to read,
and that that is a good justification for it.

Regards,
Don

-- 
  Donald C. Wells      Scientist - GBT Project        dwells@nrao.edu
                    http://www.cv.nrao.edu/~dwells
  National Radio Astronomy Observatory                +1-434-296-0277
  520 Edgemont Road,   Charlottesville, Virginia       22903-2475 USA



>From valdes Wed Jun 20 11:56:21 2001
Date: Wed, 20 Jun 2001 11:56:19 -0700 (MST)
From: Frank Valdes <valdes>
To: dwells@cv.nrao.edu
Subject: Re: WCS issue
Cc: tody, davis

Hi Don,

> .  As a contrived example, suppose you have six world axes of RA,
>  > DEC, WAVELENGTH, TIME, TEMPERATURE, and POLARIZATION.  Now suppose the
>  > RA/DEC/WAVELENGTH are coupled as in my other examples but are independent
>  > of the other axes.  Now also suppose TEMPERATURE and POLARIZATION are
>  > coupled and that TIME is non-linear.  Then in my proposal you would have
>  > 
>  > RA---TAN-PC1
>  > DEC--TAN-PC1
>  > WAVE-WAV-PC1
>  > TEMPERATURE-PC2
>  > POLARIZATION-PC2
>  > TIME-PC3
>  > 
>  > PC1 --> independent variables are intermediate coordinates from RA,
>  > 	DEC, and WAVE and cartisian radius of these
>  > PC2 --> independent variables are temperature and polarization and radius
>  > PC3 --> independent variable is time
> 
> Section 3 of your "Functional corrections and axis coupling.."
> document says ".. keywords PC1_ns, where 1 is the axis number..".
> Does this mean that the above example should be:
> 
>  > RA---TAN-PC1-PC2-PC3
>  > DEC--TAN-PC1-PC2-PC3
>  > WAVE-WAV-PC3
>  > TEMPERATURE-PC4-PC5
>  > POLARIZATION-PC4-PC5
>  > TIME-PC6
> 
> with six polynomials
> 
> PC1_ns    2D poly corrects RA using RA, DEC
> PC2_ns    2D poly corrects DEC using RA, DEC
> PC3_ns    3D poly corrects WAVE using RA, DEC, WAVE
> PC4_ns    2D poly corrects TEMPERATURE using TEMPERATURE, POLARIZATION
> PC5_ns    2D poly corrects POLARIZATION using TEMPERATURE, POLARIZATION
> PC6_ns    1D poly corrects TIME using TIME
                     http://www.cv.nrao.edu/~dwells

You are confused between the notation for the intermediate transformation
function (the suffix) and the keywords for the coefficients.  


CTYPE1 = RA---TAN-PC1
CTYPE2 = DEC--TAN-PC1
CTYPE3 = WAVE-WAV-PC1
CTYPE4 = TEMPERATURE-PC2
CTYPE5 = POLARIZATION-PC2
CTYPE6 = TIME-PC3

The number in this case is the index for three different
transformations.  PC1 is the first transformation of 3 variables (not
counting the radius variable) involving axis 1, 2, and 3 (RA, DEC, and
WAVE).  PC2 is the second transformation of 2 variables involving
axes 4 and 5, and PC3 is the third transformation involving axis 6
alone.  For this transformation suffix the meaning is PC=polynomial
correction and three different polynomials.

For the PC1 transformation the **keywords** for axes 1, 2, and 3 would be
PC1_0s, PC1_1s, ..., PC2_0s, PC2_1s, ..., and PC3_0s, PC3_1s,....
So each axes has possible keywords PCi_ns (i=axis, n=coeff#, s=multiple WCS).

Rather than use PC for both it could be

CTYPE1 = RA---TAN-FEE
CTYPE2 = DEC--TAN-FEE
CTYPE3 = WAVE-WAV-FEE
CTYPE4 = TEMPERATURE-FIE
CTYPE5 = POLARIZATION-FIE
CTYPE6 = TIME-FOE

where the coefficient keywords for transformation FEE are PC1_ns,
PC2_ns, PC3_ns and for FIE are PC4_ns, PC5_ns, and for FOE are PC6_ns.

Yours,
Frank


>From dwells@fits.cv.nrao.edu Wed Jun 20 13:17:25 2001
Date: Wed, 20 Jun 2001 16:17:22 -0400
From: Don Wells <dwells@cv.nrao.edu>
To: Frank Valdes <valdes@noao.edu>
Cc: dwells@NRAO.EDU, tody@tucana.tuc.noao.edu, davis@tucana.tuc.noao.edu
Subject: Re: WCS issue

Frank,

Frank Valdes writes:
 > You are confused between the notation for the intermediate transformation
 > function (the suffix) and the keywords for the coefficients.  
 > 
 > CTYPE1 = RA---TAN-PC1
 > CTYPE2 = DEC--TAN-PC1
 > CTYPE3 = WAVE-WAV-PC1
 > CTYPE4 = TEMPERATURE-PC2
 > CTYPE5 = POLARIZATION-PC2
 > CTYPE6 = TIME-PC3
 > 
 > The number in this case is the index for three different
 > transformations.  PC1 is the first transformation of 3 variables (not
 > counting the radius variable) involving axis 1, 2, and 3 (RA, DEC, and
 > WAVE).  PC2 is the second transformation of 2 variables involving
 > axes 4 and 5, and PC3 is the third transformation involving axis 6
 > alone.  For this transformation suffix the meaning is PC=polynomial
 > correction and three different polynomials.

 > For the PC1 transformation the **keywords** for axes 1, 2, and 3 would be
 > PC1_0s, PC1_1s, ..., PC2_0s, PC2_1s, ..., and PC3_0s, PC3_1s,....
 > So each axes has possible keywords PCi_ns (i=axis, n=coeff#, s=multiple WCS).

OK, I see, you have an implied convention.  It seems to me that the
fact that the strings 'PC1', 'PC2', .., are used for two different
purposes is going to be an endless source of confusion and potential
error.  Doesn't your implied convention also carry with it a
limitation that an axis cannot be a part of more than one
transformation?  I.e., isn't the notation that I suggested not only
more readable and less likely to lead to error, but also more general?

				 -=-

Another issue is the Cartesian radius terms. In a 2D polynomial that
is correcting a Schmidt camera image for distortion it is clear that
you can make a good numerical map by using enough terms. In a 3D
polynomial that is correcting a scanning Fabry-Perot cube, I am less
sure that the 3D radius from the reference point is the right form for
the distortion.  Have I misunderstood something here too?

regards,
Don



>From valdes Wed Jun 20 13:51:28 2001
Date: Wed, 20 Jun 2001 13:51:24 -0700 (MST)
From: Frank Valdes <valdes>
To: dwells@cv.nrao.edu
Subject: Re: WCS issue
Cc: tody, davis

Don,

>  > You are confused between the notation for the intermediate transformation
>  > function (the suffix) and the keywords for the coefficients.  
>  > 
>  > CTYPE1 = RA---TAN-PC1
>  > CTYPE2 = DEC--TAN-PC1
>  > CTYPE3 = WAVE-WAV-PC1
>  > CTYPE4 = TEMPERATURE-PC2
>  > CTYPE5 = POLARIZATION-PC2
>  > CTYPE6 = TIME-PC3
>  > 
>  > The number in this case is the index for three different
>  > transformations.  PC1 is the first transformation of 3 variables (not
>  > counting the radius variable) involving axis 1, 2, and 3 (RA, DEC, and
>  > WAVE).  PC2 is the second transformation of 2 variables involving
>  > axes 4 and 5, and PC3 is the third transformation involving axis 6
>  > alone.  For this transformation suffix the meaning is PC=polynomial
>  > correction and three different polynomials.
> 
>  > For the PC1 transformation the **keywords** for axes 1, 2, and 3 would be
>  > PC1_0s, PC1_1s, ..., PC2_0s, PC2_1s, ..., and PC3_0s, PC3_1s,....
>  > So each axes has possible keywords PCi_ns (i=axis, n=coeff#, s=multiple WCS).
> 
> OK, I see, you have an implied convention.  It seems to me that the
> fact that the strings 'PC1', 'PC2', .., are used for two different
> purposes is going to be an endless source of confusion and potential
> error.  Doesn't your implied convention also carry with it a
> limitation that an axis cannot be a part of more than one
> transformation?  I.e., isn't the notation that I suggested not only
> more readable and less likely to lead to error, but also more general?

Yes I see the possiblity for confusion.  First one needs to settle on a
short (two character) root for the function parameters.  As I said
earlier, I was going to suggest PV but then I realized that for those
basic, ideal transformations in the WCS papers where PV is used, that
would cause a conflict.  So one general rule is that at each level of
transformation they have to use different parameter roots.  PC seemed
nice since P is parameter or polynomial and C is coefficient,
correction, or coupling.  So let us use PC for the keyword root.

Then the distortion transformation would use PC parameter/coefficient
keywords.  The transformation name, which is specified as a suffix to
the CTYPE string, can be something else.  Instead of PC maybe POLYi?
[Now suggested to be PLYs. FV]

It is correct that an axis can only be associated with one distortion/coupling
transformation.  But that is not really a problem because if an axis
is required to be part of two other axes then that implies the two other
axes are also coupled and so you would use a 3D coupling.  By not using
coefficients for certain axes (defaulting to zero) it is possible to
get A=f(A,C), B=f(B,C), and C=f(C) as well as A=f(A,B,C), etc.

To do as you suggested, TRANS1-TRANS2-..., is certainly possible but then
you would need additional layers of keywords for each level of transformation.
It is more general but I don't see that the added complexity is warrented.
By allowing future additions of transformation types beyond the polynomial
one I defined, if you need POLY->TRIG->X this set of transformations would
become a new defined transformation, say PTX and would appear as -PTX.
I think describing the process in paper I should just be a single
linear transformation, followed by a single multvariable distortion/coupling
transformation, followed by the final ideal transformation to world coordinates.
I don't think we would regret not allowing an arbitrary sequence of
intermediate transformation.  However that could certainly be considered.


> Another issue is the Cartesian radius terms. In a 2D polynomial that
> is correcting a Schmidt camera image for distortion it is clear that
> you can make a good numerical map by using enough terms. In a 3D
> polynomial that is correcting a scanning Fabry-Perot cube, I am less
> sure that the 3D radius from the reference point is the right form for
> the distortion.  Have I misunderstood something here too?
> 

This is related to the question that not all variables have to appear in
every transformation.  While it is possible to use a 3D radius, and
maybe some distortion that is not the basic optical one might need it,
but in practice
one selects only certain variables to be part of the radius.  If you look
closely at the example for FP the radius is defined to be computed from

r = sqrt (c1*X^2 + c2*Y^2 + c3*Z^2)

but by setting the c3 coefficient to 0 or absent would mean that the
radius variable is actually only a 2D radius of the spatial axes.  So
the by selecting which coefficients are present one has a 2D radius
appearing only where it is needed.

Yours,
Frank



>From dwells@fits.cv.nrao.edu Wed Jun 20 15:04:07 2001
Date: Wed, 20 Jun 2001 18:04:04 -0400
From: Don Wells <dwells@cv.nrao.edu>
To: Frank Valdes <valdes@noao.edu>
Cc: dwells@NRAO.EDU, tody@tucana.tuc.noao.edu, davis@tucana.tuc.noao.edu
Subject: Re: WCS issue

Dear Frank,

Frank Valdes writes:
 > ..  First one needs to settle on a short (two character) root for
 > the function parameters..  PC seemed nice since P is parameter or
 > polynomial and C is coefficient, correction, or coupling.  So let
 > us use PC for the keyword root..

Right. I would have made the exact same choice.  The notion that you
might use some other root string in the future is exactly the sort of
thing that WCS-I is for.  I infer that you are proposing the
convention that a CTYPE value string should have only one of these
appended transformation codes?

 > It is correct that an axis can only be associated with one
 > distortion/coupling transformation.  But that is not really a
 > problem because if an axis is required to be part of two other axes
 > then that implies the two other axes are also coupled and so you
 > would use a 3D coupling.  By not using coefficients for certain
 > axes (defaulting to zero) it is possible to get A=f(A,C), B=f(B,C),
 > and C=f(C) as well as A=f(A,B,C), etc.

OK, I think that I understand it now.  I still say that using -PC1 and
-PC2 to identify two transformations and PC1_is, PC2_is, PC3_is
keywords for the first transformation and PC4_is and PC5_is keywords
for the second transformation is sure to lead to confusion (especially
for newcomers), but I agree that _it_will_work_. Several complicated
examples with careful, redundant explanations may be sufficient.  I
dare not propose using alphabetic characters for the transformation
suffix code because of possible confusion with the 's' characters.
Maybe insertion of an underscore would help, e.g. '-PC_1', but that is
inconsistent with other FITS usage. Again, the bottom line is that I
agree that your design will work, and that is what our community
desperately needs at this moment to get the WCS process moving again..

 > I don't think we would regret not allowing an arbitrary sequence of
 > intermediate transformation.  However that could certainly be
 > considered.

OK. 

 > > Another issue is the Cartesian radius terms.. I am [unsure] that
 > > the 3D radius from the reference point is the right form for the
 > > distortion.  Have I misunderstood something here too?
 > 
 > .. for FP the radius is defined to be computed from
 > 
 > r = sqrt (c1*X^2 + c2*Y^2 + c3*Z^2)
 > 
 > but by setting the c3 coefficient to 0 or absent would mean that the
 > radius variable is actually only a 2D radius of the spatial axes..

My error -- I had not understood this implications of equation (2)
even though you stated it explicitly in the text. 

My compliments -- your design is clever.

				 -=-

I am, of course, still tempted to propose orthogonal polynomials (Chebyshev
on axes, radial part of Zernike for radius) rather than ordinary
polynomials, but the fact is that when I floated that trial balloon at the
Sonthofen BoF three years ago I was shot down.  Sigh... it is a fact that
the transformations can be derived by any means we please and then
transformed to ordinary polynomials for export, and the FITS value fields
are able to convey full double precision which should generally avoid
truncation errors.  Perhaps some remarks along these lines should be added
to the text?

				 -=-

You have told me in previous messages that slit trajectories can be
'jagged' functions, and that you want splines to represent them. Are
you now content to export polynomials of order high enough to match
the capability of the splines you are using?  

Regards,
Don
-- 
  Donald C. Wells      Scientist - GBT Project        dwells@nrao.edu
                    http://www.cv.nrao.edu/~dwells
  National Radio Astronomy Observatory                +1-434-296-0277
  520 Edgemont Road,   Charlottesville, Virginia       22903-2475 USA



>From valdes Wed Jun 20 16:20:30 2001
Date: Wed, 20 Jun 2001 16:20:28 -0700 (MST)
From: Frank Valdes <valdes>
To: dwells@cv.nrao.edu
Subject: Re: WCS issue
Cc: davis, tody

Dear Don,

> You have told me in previous messages that slit trajectories can be
> 'jagged' functions, and that you want splines to represent them. Are
> you now content to export polynomials of order high enough to match
> the capability of the splines you are using?  

For the problem of coupled axes I don't know any way other than
polynomials to have a general coupling.  With polynomials you have
coeff * x^i * y^j * z^k ...  This has the advantages that with all
variables fixed but one, you get a 1D polynomial and the full function
is the product of 1D functions in each axis.  I have no problem with
other forms of polynomials.  However, with the form I proposed that
provides for normalization of the variables I don't see a real need for
other forms.

For interpolation or sampling functions it is harder to understand how
to have coupled relations.  Basically a 2D or 3D interpolation/sampling
function will be like the pixel regularization image.

I note that in IRAF when we get to 2D and higher fitting functions
we currently only use polynomials.

My attitude is that provided the door is open to define additional
intermediate transformation forms, the multi-variate polynomial is
sufficient for me.  I have no idea how to propose anything else for
coupling axes.

It is only the case of 1D functions or non-coupled axes where splines
or sampled functions might be desirable.  I don't feel very strongly
about this any more but I believe Doug still wants to think about this.

Frank   



>From dwells@fits.cv.nrao.edu Wed Jun 20 16:40:35 2001
Date: Wed, 20 Jun 2001 19:40:32 -0400
From: Don Wells <dwells@cv.nrao.edu>
To: Frank Valdes <valdes@noao.edu>
Cc: dwells@NRAO.EDU, davis@tucana.tuc.noao.edu, tody@tucana.tuc.noao.edu
Subject: Re: WCS issue

Dear Frank,

Frank Valdes writes:
 > For the problem of coupled axes I don't know any way other than
 > polynomials to have a general coupling..

Agreed.

 > ..  I have no problem with other forms of polynomials.  However,
 > with the form I proposed that provides for normalization of the
 > variables I don't see a real need for other forms..

With axis normalization agreed to, the only issue is orthogonality
(correlation) of the terms, and I think that that is really more an
esthetic issue than a technical one, for the reasons I gave in my
previous message.

 > It is only the case of 1D functions or non-coupled axes where splines
 > or sampled functions might be desirable.  I don't feel very strongly
 > about this any more but I believe Doug still wants to think about this.

The 'n' in your PCi_ns keywords can go to 9999 for the case of no
alternative WCS descriptions, and that should provide a very high
order 1D polynomial indeed. Even the 999 limit is quite generous for
1D cases (probably equivalent to a spline with knots every 4th pixel
for big CCDs). I hope that trucation is not a problem for such high
orders. 

-Don 



>From valdes Wed Jun 20 17:05:07 2001
Date: Wed, 20 Jun 2001 17:05:06 -0700 (MST)
From: Frank Valdes <valdes>
To: dwells@cv.nrao.edu
Subject: Re: WCS issue
Cc: davis, tody

Dear Don,

The problem with representing things with sharp features, such
as abrupt kinks, is not that you can't make a polynomial of a
high enough order but that the polynomial then has high frequency
wiggles.  For illustration suppose you have a sawtooth
    .
   . .
  .   .
 .     ,
.       .


This can be smoothly represented with two linear functions.  But
with a polynomial you would get wiggles.  Even if you used an
order that makes the function go through every "pixel", evaluation
of subpixel points will not match the true linear behavior very
well.  Polynomials and splines are different beasts and you can't
make one be the other by having more terms.  So it is not a question
of having a high enough order (which the 999 limit is indeed quite
high) but of having a local low order behavior (which is what splines
are meant to do).

So I don't think you can just say to use a high enough order
of polynomial and make the concern go away.  It is just hard for me
to say if this is really a problem in real data.

Frank



>From tody@noao.edu Wed Jun 20 17:57:11 2001
Date: Wed, 20 Jun 2001 17:57:03 -0700 (MST)
From: Doug Tody <tody@noao.edu>
To: Frank Valdes <valdes@noao.edu>
cc: dwells@cv.nrao.edu, Lindsey Davis <davis@noao.edu>
Subject: Re: WCS issue


I strongly agree with this (and yes we have seen such data).  In general
polynomials of more than about 6th order are likely to have serious
nonphysical errors even though they may fit the actual data points well.

There really is no question, to me at least, that one needs to support
sampled WCS (of which the spline is one case; a fully sampled curve is
another) to fit general 1D data.  Polynomials of reasonable order plus
the usual family of sampled WCS will fit 99.9% of the data we are likely
to see.


>From tody@noao.edu Wed Jun 20 20:43:30 2001
Date: Wed, 20 Jun 2001 20:43:21 -0700 (MST)
From: Doug Tody <tody@noao.edu>
To: Don Wells <dwells@cv.nrao.edu>
cc: Frank Valdes <valdes@noao.edu>, Lindsey Davis <davis@noao.edu>
Subject: Re: WCS issue

> Doug Tody writes:
>  > There really is no question, to me at least, that one needs to
>  > support sampled WCS (of which the spline is one case; a fully
>  > sampled curve is another) to fit general 1D data.  Polynomials of
>  > reasonable order plus the usual family of sampled WCS will fit
>  > 99.9% of the data we are likely to see.
> 
> So does the WCS community need to do anything new for interchange of
> such functions?  WCS-I.A ('Pixel regularization image') is defining a
> piecewise linear approximation.  Is it sufficient for these purposes?

Don - I don't want to give you a definite answer now - but we will try to
do so when we get back to you next week.  The quick answer is that none
of us here like the separate pixel regularization matrix, and we would
like to find a better way to solve this problem, or at least get things
to the point where we can do everything we want to do without having to
support a PRM (I suppose it could remain in the standard but we would like
to avoid having to use it except in the most extreme cases - certainly
one should not have to use it for common cases).

One wants the WCS to be neat and self-contained in one extension if at
all possible.  This does not necessarily mean avoiding some form of sampled
vector or matrix, just that in 99% of cases where one needs such a thing,
we would prefer to store everything in the same FITS extension.  As an
example, thus far in IRAF MWCS we have gotten by pretty well just storing
the coefficients or samples in header keywords.  For spectra, one can also
imagine schemes where a binary data vector is used - similar to a PRM,
but stored in the same extension as the other spectral vectors and possibly
part of the spectral WCS itself rather than the core.

I could say more but as you know Frank has been looking into this in-depth
recently for the spectral case and we need to chat about this locally
before I say more.  One way or the other, we will get back to you before
the end of the month as promised (completing this analysis is probably
the best way to deal with the meta-discussions currently underway on
fitswcs).  We will try to get back to you sooner if we can so that we have
a little more time to iterate.

Doug

[The latest draft includes more discussion about our thinking on the
pixel regularization image and the possiblity of sampled DCFs to
provide samplied type of distortion corrections.]