Incorrect transform from FITS with N image axes and N+1 World axes

[jehturner]

Hi @dsberry & co.,

Our pipeline at Gemini is writing out a WCS for long-slit spectra that has axes for wavelength, RA & Dec -- ie. one more axis than the 2D image array. Obviously both RA & Dec may vary along the slit (or not), depending on its position angle, so we have both CD2_2 & CD3_2 terms (in addition to CD1_1 for wavelength). When trying to display such an image with DS9 v8.0.1 or later (which was refactored to use AST), @chris-simpson and I have discovered that (in addition to a problem in DS9 itself) AST seems to implement the transforms described by the FITS header incorrectly, such that the 3rd World axis gets assigned a constant* value.

Here is a simplified example that illustrates the same behaviour using 1 pixel axis mapping to 2 World axes, rather than 2 to 3:

#include "fitsio.h"
#include "ast.h"

int main (int argc, char *argv[]) {

  fitsfile *fptr;
  AstFitsChan *fitschan;
  AstFrameSet *wcsinfo;
  char *header, *ctype="LINEAR";
  int n, nkeys, wcsaxes, crpix1, status=0;
  long naxes[1]={10};
  double crval1, crval2, cd1_1, cd2_1, pix[1], ipix[1], world[2]={};

  fits_create_file(&fptr, "!tmpast.fits", &status);

  fits_create_img(fptr, 16, 1, naxes, &status);

  wcsaxes = 2;
  crpix1 = 1;
  crval1 = 1.0;
  crval2 = 1.0;
  cd1_1 = 1.0;
  cd2_1 = 1.0;

  fits_update_key(fptr, TINT, "WCSAXES", &wcsaxes, NULL, &status);
  fits_update_key(fptr, TINT, "CRPIX1", &crpix1, NULL, &status);
  fits_update_key(fptr, TSTRING, "CTYPE1", ctype, NULL, &status);
  fits_update_key(fptr, TSTRING, "CTYPE2", ctype, NULL, &status);
  fits_update_key(fptr, TDOUBLE, "CRVAL1", &crval1, NULL, &status);
  fits_update_key(fptr, TDOUBLE, "CRVAL2", &crval2, NULL, &status);
  fits_update_key(fptr, TDOUBLE, "CD1_1", &cd1_1, NULL, &status);
  fits_update_key(fptr, TDOUBLE, "CD2_1", &cd2_1, NULL, &status);

  fits_hdr2str(fptr, 0, NULL, 0, &header, &nkeys, &status);

  fits_close_file(fptr, &status);

  astBegin;

  fitschan = astFitsChan( NULL, NULL, "" );
  astPutCards(fitschan, header);
  free(header);
  wcsinfo = astRead(fitschan);

  astShow(wcsinfo);
  printf("\n");

  for(n=1; n < 4; n++) {

    pix[0] = (double) n;

    /* Fwd: */
    astTranN(wcsinfo, 1, 1, 1, pix, 1, 2, 1, world);

    /* Inv: */
    astTranN(wcsinfo, 1, 2, 1, world, 0, 1, 1, ipix);

    printf("Pix %.3f -> World (%.3f, %.3f) -> pix %.3f\n",
           pix[0], world[0], world[1], ipix[0]);
    printf("       - should be (%.3f, %.3f)\n\n",
        (pix[0]-crpix1)*cd1_1 + crval1, (pix[0]-crpix1)*cd2_1 + crval2);

  }

  astEnd;

}

The output is as follows:

 Begin FrameSet 	# Set of inter-related coordinate systems
#   Title = "2-d coordinate system" 	# Title of coordinate system
#   Naxes = 2 	# Number of coordinate axes
#   Domain = "LINEAR-LINEAR" 	# Coordinate system domain
#   Lbl1 = "LINEAR" 	# Label for axis 1
#   Lbl2 = "LINEAR" 	# Label for axis 2
 IsA Frame 	# Coordinate system description
    Nframe = 2 	# Number of Frames in FrameSet
    Base = 1 	# Index of base Frame
    Currnt = 2 	# Index of current Frame
    Lnk2 = 1 	# Node 2 is derived from node 1
    Frm1 = 	# Frame number 1
       Begin Frame 	# Coordinate system description
          Title = "Pixel Coordinates" 	# Title of coordinate system
          Naxes = 1 	# Number of coordinate axes
          Domain = "GRID" 	# Coordinate system domain
#         Lbl1 = "Pixel axis 1" 	# Label for axis 1
          Ax1 = 	# Axis number 1
             Begin Axis 	# Coordinate axis
                Label = "Pixel axis 1" 	# Axis Label
             End Axis
       End Frame
    Frm2 = 	# Frame number 2
       Begin Frame 	# Coordinate system description
          Ident = " " 	# Permanent Object identification string
       IsA Object 	# AST Object
#         Title = "2-d coordinate system" 	# Title of coordinate system
          Naxes = 2 	# Number of coordinate axes
          Domain = "LINEAR-LINEAR" 	# Coordinate system domain
#         Lbl1 = "LINEAR" 	# Label for axis 1
#         Lbl2 = "LINEAR" 	# Label for axis 2
          Ax1 = 	# Axis number 1
             Begin Axis 	# Coordinate axis
                Label = "LINEAR" 	# Axis Label
                Symbol = "LINEAR" 	# Axis symbol
             End Axis
          Ax2 = 	# Axis number 2
             Begin Axis 	# Coordinate axis
                Label = "LINEAR" 	# Axis Label
                Symbol = "LINEAR" 	# Axis symbol
             End Axis
       End Frame
    Map2 = 	# Mapping between nodes 1 and 2
       Begin CmpMap 	# Compound Mapping
          Nin = 1 	# Number of input coordinates
          Nout = 2 	# Number of output coordinates
       IsA Mapping 	# Mapping between coordinate systems
          MapA = 	# First component Mapping
             Begin PermMap 	# Coordinate permutation
                Nin = 1 	# Number of input coordinates
                Nout = 2 	# Number of output coordinates
             IsA Mapping 	# Mapping between coordinate systems
                Out1 = 1 	# Output coordinate 1 = input coordinate 1
                Out2 = -1 	# Output coordinate 2 = constant no. 1
                In1 = 1 	# Input coordinate 1 = output coordinate 1
                Nconst = 1 	# Number of constants
                Con1 = 1 	# Constant number 1
             End PermMap
          MapB = 	# Second component Mapping
             Begin WinMap 	# Map one window on to another
                Nin = 2 	# Number of input coordinates
                IsSimp = 1 	# Mapping has been simplified
             IsA Mapping 	# Mapping between coordinate systems
                Sft2 = 1 	# Shift for axis 2
             End WinMap
       End CmpMap
 End FrameSet

Pix 1.000 -> World (1.000, 2.000) -> pix 1.000
       - should be (1.000, 1.000)

Pix 2.000 -> World (2.000, 2.000) -> pix 2.000
       - should be (2.000, 2.000)

Pix 3.000 -> World (3.000, 2.000) -> pix 3.000
       - should be (3.000, 3.000)

It appears that the PermMap is assigning Out2 to a constant rather than input coordinate 1 for some reason, though I'm not very familiar with all the AST structures, which seem quite complicated.

This happens with both AST 8.6.3 (used in DS9) and the latest 9.1.1. I haven't tested the latest master branch because it seems to require some Starlink versions of autotools stuff, which I'm no expert with generally, so it's not trivial to build.

*In our real-world example, I do actually see a small amount of variation in the 3rd World co-ordinate for very large pixel offsets from CRPIX2, but it remains almost constant for valid offsets, with the PermMap assigning that axis to a constant in the same way as above.

Am I right in thinking this is a bug / limitation, or are we doing something wrong here? Our FITS header seems valid to me.

Thanks!

James.

[dsberry]

Hi @jehturner. Thanks for pointing this out. I think the problem has two parts - one in AST and one in your FITS headers.

1. AST: If a FITS-WCS header defines more WCS axes than pixel axes (i.e. the header does NOT use the usual trick of defining a degenerate pixel axis that spans only one pixel) then, internally, the FitsChan class invents extra degenerate pixel axes itself and feeds them a value of 1.0 when doing the forward (pixel->sky) transformation. However, the default value for CRPIX specified in the FITS-WCS papers is 0.0, not 1.0. So it is more appropriate for FitsChan to feed missing pixel axes a value of 0.0 rather than 1.0. I have made this change to AST (available now in v9.1.2 on github or the AST web page).

2. Your headers: again it's a defaults issue. FITS-WCS paper I says that the default value for all missing CDi_j keywords is zero. So in your case a value of zero is used for CD1_2 and CD2_2. This means the CD matrix cannot be inverted, which breaks FITS-WCS requirements. The answer is to assign a value to CD2_2 in your header. Since the missing second pixel axis always has a value zero (if you use the fix to AST described above), the actual value of CD2_2 doesn't matter since it is always multiplied by zero. But a non-zero CD2_2 value is required to make the CD matrix invertable. This is all discussed further in FITS-WCS paper II section 7.4.3 (particularly the paragraphs preceding equation 195).

So I think the answer is a) use AST V9.1.2 and b) set CD2_2 = 1 in your header.

David

[jehturner]

Thanks, @dsberry; those comments have helped us understand the problem a bit better (sorry for overlooking some key details from the WCS papers). In addition to your change of default for CRPIX, I see that you have changed how the number of axes is determined, based on CRVAL rather than CRPIX. Without that change, the CD matrix was not getting included in the above AST transformation, until I added CRPIX explicitly, thus the second World co-ordinate was not just off by 1 at (the implicit) CRPIX2 -- it also was not varying with the pixel co-ordinate as it should have.

Now my above example works with 9.1.2, but I think we should have defined CRPIX2 as well as CD2_2 in our header, to go with the implicit degenerate axis and define a reversible transformation fully. Since the degenerate axis has a single pixel co-ordinate of 1 and CD2_2 is now 1.0, we'll need CRPIX2=1 to avoid affecting the second World co-ordinate. This is all fine and can be made to work, but with CRPIX2 present, astTranN now expects 2 input co-ordinates (although NAXIS=1), so I have to feed it a second dummy input of 1. Shouldn't it accept a single input when there's only one real pixel axis and fix the degenerate second co-ordinate at 1 internally, similar to the above? Might the current behaviour not break assumptions based on NAXIS in application software? Anyway, this is not causing us any immediate problem -- in fact, it seems to be what DS9 expects and I think adding CD3_3 and CRPIX3 to our real data resolves our ds9 problem for now. Thanks again for providing a new release to try so quickly!

James.

[dsberry]

The change to the way the number of axes is determined only affects the code that normalises the input headers into the form expected by the rest of the FitsChan class - specifically the bit that converts a CD matrix into an equivalent PC matrix. The previous version was failing to find and convert all elements of the CD matrix. The primary determination of the number of pixel axes is still based on the presence of CRPIX headers (done at the end of function FitsToStore).

If you choose to include CRPIX2 in the header then you get a 2-input Mapping. So it's then up to you to assign a suitable value to both inputs when calling astTranN. If you omit CRPIX2, then the FitsChan class effectively does what you suggest - it automatically adds an extra pixel axis internally and assigns it the constant value zero.

I'm not sure this answers your question though. The point about breaking assumptions based on NAXIS was why the WCSAXES keyword was introduced. I'm trying to think of an explicit situation where the current behaviour would be problematic.