Dealing with inconsistent WCS data

[tiagopereira]

Perhaps I'm just unlucky, but the first time I try to use the ndcube spectrograph class I get this error:

>>> s = read_iris_spectrograph_level2_fits("iris_l2_20180210_033712_3690115104_raster_t000_r00000.fits")
(...)
SingularMatrixError: ERROR 3 in wcsset() at line 2218 of file cextern/wcslib/C/wcs.c:
Linear transformation matrix is singular.
ERROR 3 in linset() at line 679 of file cextern/wcslib/C/lin.c:
PCi_ja matrix is singular.

I traced the error and it seems like WCS from ndcube.utils.wcs does not like the WCS matrix from that file, which is:

PC1_1   =        1.00000000000 /
PC1_2   =        0.00000000000 /
PC2_1   =        0.00000000000 /
PC2_2   =       0.999938666821 /
PC3_1   =        0.00000000000 /
PC3_2   =      0.0112707177177 /
PC3_3   =       0.999938666821 /
PC2_3   =     -0.0112707177177 /

This does NOT seem singular to me. In any case, read_iris_spectrograph_level2_fits should not bail out. It should deal gracefully with the error and continue (even if without WCS).

[DanRyanIrish]

Hi @tiagopereira. Could supply the full header for this file?

[tiagopereira]

The FITS file is here. Header is:

XTENSION= 'IMAGE   '           / IMAGE extension
BITPIX  =                   16 / Number of bits per data pixel
NAXIS   =                    3 / Number of data axes
NAXIS1  =                 1027 /
NAXIS2  =                  548 /
NAXIS3  =                   15 /
PCOUNT  =                    0 / No Group Parameters
GCOUNT  =                    1 / One Data Group
BSCALE  =                 0.25 /
BZERO   =                 7992 /
CDELT1  =      0.0259600002319 /
CDELT2  =             0.332700 /
CDELT3  =        0.00000000000 /
CRPIX1  =              1.00000 /
CRPIX2  =              274.500 /
CRPIX3  =              8.00000 /
CRVAL1  =        1331.68977015 /
CRVAL2  =             -6.80353 /
CRVAL3  =             -178.529 /
CTYPE1  = 'WAVE    '           /
CTYPE2  = 'HPLT-TAN'           /
CTYPE3  = 'HPLN-TAN'           /
CUNIT1  = 'Angstrom'           /
CUNIT2  = 'arcsec  '           /
CUNIT3  = 'arcsec  '           /
PC1_1   =        1.00000000000 /
PC1_2   =        0.00000000000 /
PC2_1   =        0.00000000000 /
PC2_2   =       0.999938666821 /
PC3_1   =        0.00000000000 /
PC3_2   =      0.0112707177177 /
PC3_3   =       0.999938666821 /
PC2_3   =     -0.0112707177177 /

[DanRyanIrish]

It looks like a corrupted header of some kind... Is the problem here:

PC3_3   =       0.999938666821 /
PC2_3   =     -0.0112707177177 /

I don't see equivalent entries in the header of another IRIS SJI FITS file I have locally.

[tiagopereira]

Why is that a problem? This is not an SJI file, it is a spectrograph file. That header is in the first extension. I looked at other spectrograph files and the headers looked very similar, no missing entries.

[DanRyanIrish]

Ah, I thought it was SJI. Never mind then.

[tiagopereira]

I found the problem here. It affects all sit and stare raster files, and happens because CDELT3=0. It seems that astropy.wcs.WCS does not like this and throws a singular matrix error. I think in these cases CTYPE3 should be Time instead of HPLN-TAN, but we can't change IRIS's format, so must find a workaround instead. I don't know what is best. A simple fix would be to set CDELT=1e-10 or something. We could also set CTYPE3=TIME.

[DanRyanIrish]

Thanks @tiagopereira for figuring this out. That does sound like you had to go digging in the weeds!

I agree with what say about the difference between the ideal and practical solutions. I’m not sure what the best way forward is either. I’ll give it sone thought. But perhaps your idea of setting CDELT3 to 1e-10 will end up being the best though...

[steinhh]

Sorry for getting into this discussion extremely late, but here's a rant from the FITS fanatic: The ideal solution is to fix the FITS files. They're not set in stone. It's perfectly doable - Martin sits next-door to me :). It will just take some time before the change is effectuated throughout the data set through a reprocessing. Those FITS files are invalid according to the FITS standard - which is a bad, bad thing. PCi_j matrices should be invertible. Doesn't hurt to have a temporary work-around, but having invalid WCS info may throw a wrench into third-party applications that don't care which instrument this file comes from.

I suppose it's desirable to make it possible for third-party applications to show IRIS data alongside data from other instruments, and that may require valid WCS info.

CTYPE3 should be HPLN-TAN, by the way. How else are you going to specify the HPLN-TAN coordinate for those observations? Values for TIME can (and should, really) be specified in one or both of these ways: Either as a "phantom dimension", i.e. CTYPE4='TIME' and with the appropriate PCi_j cross-terms, or in an alternate coordinate system. Regarding the "should, really" parenthesis: If knowing relative time values for each exposure of an SJI sequence image is important for scientific analysis, then this is also the case (believe it or not) for these files.

[DanRyanIrish]

Hi @steinhh. Thanks for your contribution. I completely agree that the "correct" thing to do is to fix the FITS headers, although that is beyond the competency of IRISpy. However, if the IRIS team were willing to do this that would solve the problem not only for IRISpy, but as you say, for any instrument-non-specific routines that rely on valid FITS headers. I would very much advocate for this option with a temporary fix in IRISpy until it is propagated through the IRIS dataset.

[byrdie]

Is there currently a workaround for this problem? I see that @tiagopereira suggests "setting CDELT3 to 1e-10" but I do not know how this needs to be implemented. I am currently trying to do statistical analyses on large portions of the IRIS dataset, and this issue is slightly annoying...

[DanRyanIrish]

Hi @byrdie. Are you using IRISpy in your work or does this issue have consequences beyond IRISpy? There is currently not a fix but I'll have a quick look to see how this fix could be implemented.