Singular matrix error when setting WCS

[rhiannonlynne]

Description

I have a series of images, all of which produce an error when trying to set a WCS:

SingularMatrixError: ERROR 3 in wcsset() at line 2395 of file cextern/wcslib/C/wcs.c:
Linear transformation matrix is singular.
ERROR 3 in linset() at line 677 of file cextern/wcslib/C/lin.c:
PCi_ja matrix is singular.

Expected behavior

I expected to get an astropy.wcs.WCS object to use to convert sky and pixel coordinates.

Actual behavior

When reading the image either with CCDData (which creates a WCS automatically) or with astropy.io.fits.open(), and then passing the header to astropy.wcs.WCS myself, I get the SingularMatrixError message above.

I also get a related error if I just run wcslint on the fits file directly:

  WCS key ' ':
    - RADECSYS= 'FK5' / Coordinate system, per TCC ObjSys
      the RADECSYS keyword is deprecated, use RADESYSa.
    - 'datfix' made the change 'Invalid parameter values: MJD-OBS and
      DATE-OBS are inconsistent'.
    - 'celfix' made the change 'Linear transformation matrix is
      singular'.

HDU 1 (LAST_FRAME):
  WCS key ' ':
    - 'datfix' made the change 'Invalid parameter values: MJD-OBS and
      DATE-OBS are inconsistent'.
    - 'celfix' made the change 'Linear transformation matrix is
      singular'.

Steps to Reproduce

Example image is available at https://epyc.astro.washington.edu/~lynnej/astropy_issue2_image.fits

Either just opening the file and passing the header to WCS directly, or first editing the header for the issues that wcslint reported result in the SingularMatrixError.

import os
from astropy.io import fits 
from astropy import wcs

d = fits.open(os.path.join(ddir, imname))
header = d[0].header
del header['MJD-OBS']
header['RADESYSa'] = header['RADECSYS']
del header['RADECSYS']
header
wcs.WCS(header, naxis=(1,2))

System Details

Darwin-19.2.0-x86_64-i386-64bit
Python 3.7.6 (default, Jan 8 2020, 13:42:34)
[Clang 4.0.1 (tags/RELEASE_401/final)]
Numpy 1.18.1
astropy 4.0
Scipy 1.3.1

[bsipocz]

Running the example snippet on master I get the same traceback (with wcslib line numbers slightly changed):

~/munka/devel/astropy/astropy/wcs/wcs.py in __init__(self, header, fobj, key, minerr, relax, naxis, keysel, colsel, fix, translate_units, _do_set)
    497 
    498             if naxis is not None:
--> 499                 wcsprm = wcsprm.sub(naxis)
    500             self.naxis = wcsprm.naxis
    501 

SingularMatrixError: ERROR 3 in wcsset() at line 2507 of file cextern/wcslib/C/wcs.c:
Linear transformation matrix is singular.
ERROR 3 in linset() at line 677 of file cextern/wcslib/C/lin.c:
PCi_ja matrix is singular.

[dhomeier]

specutils had a similar problem showing up on astropy/specutils#160.

[mcara]

I wonder if this ever worked or if it is supposed to even work. A quick look at the header and WCSLIB code shows that the presence of CNPIX in the header triggers special processing in WCSLIB for DSS. I am not familiar with this format but it is likely that the header simply does not comply with expected description of the WCS.

Specifically, it looks like there is not enough information to infer plate scale and WCSLIB defaults to 0 => singular PC matrix.

The "fix" below is a nonsense and I have no idea what it means but it makes your header not to crash when creating a WCS object:

In [11]: header['XPIXELSZ'] = 1
    ...: header['YPIXELSZ'] = 1
    ...: for k in range(5):
    ...:     header['AMDY{:d}'.format(k)] = np.random.rand()
    ...:     header['AMDX{:d}'.format(k)] = np.random.rand()
    ...:     header['AMDREX{:d}'.format(k)] = np.random.rand()
    ...:     header['AMDREY{:d}'.format(k)] = np.random.rand()
    ...: wcs.WCS(header)
WARNING: FITSFixedWarning: RADECSYS= 'FK5' / Coordinate system, per TCC ObjSys
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Invalid parameter value: inconsistent date '2020-01-13T01:08:52.636000''. [astropy.wcs.wcs]
Out[11]:
WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN'  'DEC--TAN'
CRVAL : 1.0432582341946395e-308  2.14739352155e-313
CRPIX : 1290.6021065659497  671.0189656991273
PC1_1 PC1_2  : 0.0008094668730622751  -0.0012388367513387422
PC2_1 PC2_2  : -0.0013736800038028686  0.0008669474893213419
CDELT : -0.00014110467379405977  0.00014110467379405977
NAXIS : 1024  1024

---

Alternatively, remove CNPIX1 and CNPIX2 header keywords to fix the crash:

In [12]: d = fits.open(os.path.join(ddir, imname))
    ...: header = d[0].header
    ...: del header['CNPIX*']
    ...: del header['MJD-OBS']
    ...: header['RADESYSa'] = header['RADECSYS']
    ...: del header['RADECSYS']
    ...: header
    ...: wcs.WCS(header, naxis=(1,2))
WARNING: FITSFixedWarning: 'datfix' made the change 'Set DATE-REF to '1858-11-17' from MJD-REF.
Set MJD-OBS to 58861.047831 from DATE-OBS'. [astropy.wcs.wcs]
Out[12]:
WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN'  'DEC--TAN'
CRVAL : 326.07266745682  -1.3576216518211
CRPIX : 475.2301142  475.2737288
CD1_1 CD1_2  : 7.555095534183e-05  0.0
CD2_1 CD2_2  : -0.0  7.5640676530211e-05
NAXIS : 1024  1024

[mcara]

@dhomeier in that specutils example CDELT is set to 0. I wonder if that is the reason for the crash.

[pllim]

If it is indeed the case of invalid WCS, then the expectation is that the user fixes the data. astropy.wcs should not be making assumptions for invalid cases (let's learn that lesson from FITS!).

[Cadair]

in that specutils example CDELT is set to 0. I wonder if that is the reason for the crash.

Yes, this causes this error. I have seen it in some more frustrating solar data.

[dhomeier]

in that specutils example CDELT is set to 0. I wonder if that is the reason for the crash.

Yes, this causes this error. I have seen it in some more frustrating solar data.

OK, that does sound bad. We had been discussing if explicitly setting CDELT to 1 would help in the present case, but it seemed counterintuitive, as the respective wcs->cdelt[i] appear to be either explicitly set to 1 or tested for 0 already within wcsset.

So is there a good reference on what constitutes a valid set of WCS keywords? Here the best take seems to be @mcara's suggestion to remove CNPIX1 and CNPIX2, but has anyone an idea why? The errors from wcslint and wcslib are not terribly illuminating...

[mcara]

Here the best take seems to be @mcara's suggestion to remove CNPIX1 and CNPIX2, but has anyone an idea why?

Because removing these keywords does not any longer trigger special DSS WCS processing in the WCSLIB. Then your WCS is assumed to be a regular FITS WCS and so stuff from the WCS papers (I, II, III, etc.) applies to your (modified) WCS. Specifically, in this case, when CNPIX1,2 are removed, WCSLIB is no longer trying to read AMD* coefficients (I have no idea what they are but it looks like they are assumed to be 0 when missing) and instead it is loading CD matrix specified in the @rhiannonlynne's header (which is no longer singular).

[nden]

@dhomeier A valid set of WCS keywords is defined in the WCS papers which is what wcslib imlements. wcslib also is capable of reading some non-standard conventions. CNPIX1/2 are definitely part of a non-standard convention but I am not sure which one. Perhaps @rhiannonlynne can point us to this WCS convention? All I found through google is https://iraf.net/forum/viewtopic.php?showtopic=104279

[rhiannonlynne]

Thanks @mcara .. deleting the "CNPIX1" and "CNPIX2" keywords did indeed let the header be parsed by WCS.

d = fits.open(os.path.join(ddir, imname))
header = d[0].header
del header['MJD-OBS']
header['RADESYSa'] = header['RADECSYS']
del header['RADECSYS']
del header['CNPIX1']
del header['CNPIX2']
header
wcs.WCS(header, naxis=(1,2))

resulted in

WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN'  'DEC--TAN'  
CRVAL : 326.07266745682  -1.3576216518211  
CRPIX : 475.2301142  475.2737288  
CD1_1 CD1_2  : 7.555095534183e-05  0.0  
CD2_1 CD2_2  : -0.0  7.5640676530211e-05  
NAXIS : 1024  1024

Why doesn't wcslint or WCS itself complain about these header keys though?

As far as what convention it was based on ... I don't know.
The images are from NIC-FPS (NIR camera) at Apache Pointe Observatory.

[mcara]

Why doesn't wcslint or WCS itself complain about these header keys though?

Probably because CNPIX1 and CNPIX2 keywords are supported for DSS images. When I look at your header I see these keywords and I immediately wonder whether it is a DSS image. Then I look for AMD* and ?PIXELSZ coefficients (because if, i.e., CNPIX1 is present then I know I am supposed to look for them) but I do not find any. It seems that the default action in that case is to assume all missing coefficients are zero => singular matrix => error that you see. WCSLIB cannot know that you added CNPIX1 and CNPIX2 "unintentionally". Also, if CNPIX? is present, CD (or PC) matrix seems to be computed from AMD* and ?PIXELSZ coefficients instead of being loaded from the header. Why? I have no idea. Maybe that is how DSS headers are supposed to work.

[pllim]

Can someone ping someone who works for APO?

[rhiannonlynne]

I think my question should be rephrased as : "given that CNPIX* header keywords imply that you should have additional specific keys, could an error be given that indicates such?"
something like
"CNPIXi header keywords imply that a DSS WCS is in use and requires AMD* or ?PIXELSZ coefficients, but these are missing."

I can send a message to our local APO people and see if I can find out more.

[mcara]

Ah, I see. That is a tough question. We had these kind of issues with SIP and it is a mess (a sample of the kind of discussion is here #4761 but there are many more issues). If "standard" (DSS, APO???) says it is OK to have missing coefficients then there is little that can be done. But this (CNPIX) is totally out of my league.

[mcara]

"CNPIXi header keywords imply that a DSS WCS is in use and requires AMD* or ?PIXELSZ coefficients, but these are missing."

That's the thing: is it actually required to be present?

[rhiannonlynne]

I don't know, as I can't find a standard anywhere. However, if not having it crashes wcslib, then I think it's required for wcslib at least.
And if it's not actually required but somehow ending up in the SingularMatrixError because the CDi_j values aren't being evaluated, then that seems like a bug.

[mcara]

if not having it crashes wcslib, then I think it's required for wcslib at least

In that case the fix for your particular issue should be in not deleting CNPIX1,2 but rather in adding the missing ?PIXELSZ and AMD* coefficients.

[mcara]

if it's not actually required but somehow ending up in the SingularMatrixError because the CDi_j values aren't being evaluated, then that seems like a bug.

I do not think it works like that. I think the CD matrix is not being read because it has detected singular matrix before even getting to the CD matrix or maybe regardless of CD matrix. For example, if the matrix "A" computed from AMD* is null and total coordinate transformation matrix is A*CD, then the resultant matrix will be singular regardless of CD.

[mcara]

For reference: https://archive.stsci.edu/dss/booklet_n.pdf

Take a look at how xi and eta are defined: if first coefficients are 0 then the linear transformation matrix ought to be null.

I, in particular, like:

The above information is provided for those who wish to compute object coordinates once a subimage has been uncompressed. We suggest examining the functions readhr (found in the file header.c of the GetImage software) and amdpos (found in the file astrmcal.c of the GetImage software) to see how the keywords are used in an actual algorithm.

[mcara]

To me, at this point, this issue looks more like a bug in whatever software was used to create the WCS in @rhiannonlynne's image: either CNPIX? was accidentally left over or AMD* are not properly set.

[mcara]

Here is the software that illustrates this kind of WCS: https://gsss.stsci.edu/Software/GetImage/GetImage.htm

[pllim]

Looks like the data creator made invalid WCS? Should we close?

[github-actions]

Hi humans 👋 - this issue was labeled as Close? approximately 17 hours ago. If you think this issue should not be closed, a maintainer should remove the Close? label - otherwise, I will close this issue in 7 days.

If you believe I commented on this issue incorrectly, please report this here

[github-actions]

I'm going to close this issue as per my previous message, but if you feel that this issue should stay open, then feel free to re-open and remove the Close? label.

If this is the first time I am commenting on this issue, or if you believe I closed this issue incorrectly, please report this here