Update WCSLIB to 7.4 with fix for wcshdo export of complete TPD coefficients

[dhomeier]

Description

#11216 (comment) described an inconsistency in WCSLIB's iparm parameter lists in dis.c and wcshdr.c, introduced somewhere between 6.4 (c754932) and 7.1 (56fbc3d), when iparm[I_DOCORR] was removed from dis.c, but not in wcshdr.c. As a result, that parameter is not correctly checked to switch on writing of the AUX.* and TPD.FWD.* coefficients by wcs.to_header(), which consequently creates an incomplete distortion model (applying no distortion at all) in the header.
The differences to the correct WCS with distortion are found when directly creating a new WCS from that header or from a Cutout2D produced from the original, either directly or when written to file and read back (see added test_distortion_header for examples).

Fixes part of #11216 (old files written with missing 'DOCORR: 0' keyword will still not be read in correctly with current WCSLIB, which has 'DOCORR: 1' as default since 6.3).
Since we want to create correct headers when linking a system WCSLIB, this requires an upstream fix, but as I don't know the best way to send patches to https://www.atnf.csiro.au/people/mcalabre/WCS/, I have added a proposed fix in the bundled wcslib here as a suggestion.

I'd also like to propose keeping the info comment à la Translated from DSS to TPD by WCSLIB 7.3.0 in wcs.to_header(), as the issues in #11216 show that this information can be rather useful. This would require updating a number of tests checking for (no) COMMENT cards in the header though, so maybe it could be alternatively be written to its own card like
WCSVER = 7.3.0 / Translated from DSS to TPD by WCSLIB 7.3.0

[mcara]

Even if we are going to handle this internally in astropy, should we patch WCSLIB here in cextern ? I think this would make future maintenance more difficult unless this is going to be fixed upstream in the next release and this PR is a temporary fix.

[dhomeier]

Yes, this is only meant as a kind of template how to fix it in WCSLIB. Ideally upstream could apply something equivalent not too long from now, then the new version might even be added to our cextern as part of this PR.
Otherwise I think we could not fully fix this internally as I don't believe we can really control if people build with ASTROPY_USE_SYSTEM_WCSLIB=1, and anyway linking the system is the official recommendation now, right?
If the upstream fix should take too long, all that could be done internally is probably checking if wcs.to_header() has output the complete coefficients, and issue a warning otherwise.

In addition, if 'DOCORR: 1' is confirmed to be taken as default, and there really still are that many files around created with an earlier version, we should discuss again if we keep it at some instructions in the docs how to fix those manually, or offer a fixing option for WCS(hdr).

Setting up an email to mcalabre to discuss the first part.

[mcara]

@dhomeier Thank you!

[pllim]

I am tentatively marking this for 4.3. Thanks!

[dhomeier]

Big thanks to Mark Calabretta for quickly merging the patch and getting it into WCSLIB 7.4, released yesterday!
I have updated our bundled copy now (tested with both the bundled and external 7.4 under macOS Mojave).

[dhomeier]

gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -fPIC -DECHO -DWCSTRIG_MACRO -DASTROPY_WCS_BUILD -D_GNU_SOURCE -DNDEBUG -DHAVE_SINCOS -UDEBUG -I/tmp/pip-build-env-47a6b2im/overlay/lib/python3.7/site-packages/numpy/core/include -Icextern/wcslib/C -Iastropy/wcs/include -I/__w/astropy/astropy/.tox/py37-test/include -I/opt/python/cp37-cp37m/include/python3.7m -c cextern/wcslib/C/wcshdr.c -o build/temp.linux-x86_64-3.7/cextern/wcslib/C/wcshdr.o -Wno-strict-prototypes -Wno-unused-function -Wno-unused-value -Wno-uninitialized
1467
cextern/wcslib/C/wcshdr.c: In function ‘wcshdo_format’:
1468
cextern/wcslib/C/wcshdr.c:1949:3: error: ‘for’ loop initial declarations are only allowed in C99 mode
1469
for (int i = 0; i < nval; i++) {
1470
^
1471
cextern/wcslib/C/wcshdr.c:1949:3: note: use option -std=c99 or -std=gnu99 to compile your code
1472
error: command 'gcc' failed with exit status 1

Anyone familiar with the compiler on CI / 32-bit?

[dhomeier]

Was adding FITSFixedWarning to the base-level ignores the right way to deal with it? Seems a bit blunt.

In addition to the library update, I made one change to wcs.to_header saving the version info from the wcshdo comment to a new keyword WCSLIBV (keeping it as a comment created massive trouble with various tests, and in principle the version is easier to retrieve this way). Up to discussion of course, but the original issue shows that having the WCSLIB version stored in the file could be quite useful. One point for keeping it as a comment would be consistency with other applications using WCSLIB, which will likely have kept the header in its original state.

And the original issue in #11216 is still open, since an old-style WCS distortion header is still read in assuming DOCORR=1 as default. That leaves some documentation work to point out the different versions and describe how to fix older files, either through a convenience method, or by documenting an example.
Somewhat related, I found that some 3rd-party software may not be able to understand the TPD distortion model (yet) – e.g. the current SAOImageDS9 apparently can correctly process the AMD[XY]* coefficients in original DSS files, but not files written back with the TDP coefficients, simply ignoring distortion for the latter. This is not nearly as bad as mixing up corrections and final coordinates, but can produce offsets of a few arcsec over a plate field. So it should be documented as well (for plain image copies and some cutouts one could write back the original header to avoid this, but Cutout2D will in general require a new header from wcs.to_header()).

[mcara]

Anyone familiar with the compiler on CI / 32-bit?

Not familiar, but I suggest we declare int i at the beginning of the function. This loop seems to be the only loop in the entire WCSLIB C code that declares loop variable inside the loop. Most likely this is a "bug/typo". Fixing this in the WCSLIB code should not pose problems here, I think.

[mcara]

Was adding FITSFixedWarning to the base-level ignores the right way to deal with it? Seems a bit blunt.

@nden What do you think?

In addition to the library update, I made one change to wcs.to_header saving the version info from the wcshdo comment to a new keyword WCSLIBV (keeping it as a comment created massive trouble with various tests, and in principle the version is easier to retrieve this way). Up to discussion of course, but the original issue shows that having the WCSLIB version stored in the file could be quite useful. One point for keeping it as a comment would be consistency with other applications using WCSLIB, which will likely have kept the header in its original state.

First, I did not know that we removed comments obtained from to_header().

Second, about WCSLIBV... I wish FITS keywords could be longer. Then, while I think this is, in principle, a great addition, there may be issues the header returned from to_header is going to be merged with a FITS header that contains several "alternate" WCSes each created with a different version or with an older version for which we did not write WCSLIBV before. Definitely, we will need to take care of this in stwcs.altwcs. I am not sure yet if there will be issues in astropy itself. CC: @nden

In any case, I would prefer to have this addition be done in a different PR separate from WCSLIB update.

[dhomeier]

First, I did not know that we removed comments obtained from to_header().

Second, about WCSLIBV... I wish FITS keywords could be longer. Then, while I think this is, in principle, a great addition, there may be issues the header returned from to_header is going to be merged with a FITS header that contains several "alternate" WCSes each created with a different version or with an older version for which we did not write WCSLIBV before. Definitely, we will need to take care of this in stwcs.altwcs. I am not sure yet if there will be issues in astropy itself. CC: @nden

As I read the discussion back in #3905 the primary incentive was to remove the empty lines, and beyond that, that COMMENT keywords are not part of the WCS standard, which I don't entirely follow... But it sounds indeed as if the comment may still cause problems in stwcs.

Anyway keeping the version info as a COMMENT instead requires two more wcs tests to be adapted, easily enough dealt with, plus test_wcs_keyword_removal_for_wcs_test_files in nddata. So either COMMENT would have to be added to the list of ccddata._KEEP_THESE_KEYWORDS_IN_HEADER or the test be modified somehow.
Not sure what's the better solution, and this could probably touch on your concerns about having multiple "competing" WCSes. Though I think in general any number of comment keys in a header should not break any reader.
Update: the error arose in fact from (other) COMMENT cards in the original header from 2wcses.hdr; so I tend to ignore comments in the union set(new_header) & set(new_wcs_header) for the new PR.

[mcara]

I do not think comments should cause technical problems. However, imagine you will get a file with multiple (let's say 4) WCSes in the header and when you inspect that file you will find:

COMMENT  WCS header keyrecords produced by WCSLIB 5.8
COMMENT  WCS header keyrecords produced by WCSLIB 7.1
COMMENT  WCS header keyrecords produced by WCSLIB 7.3

How would you know which WCS was created with what WCSLIB. I think it just creates confusion and does not really help unless you always have only one WCS. By the way, in the example above, assume 4th WCS was created "by hand" by an user directly editing the header (or by some function in a third-party package like stwcs or astropy's Cutout2D).

I think making 'WCSVERINFO' (or a shorter version) a part of the standard would help. Then you could assign some values like WCSLIB 5.8 or STWCS 1.0.0 or astropy Cutout2D version ... That would be helpful.

Maybe we should add another parameter to to_header(..., return_wcslibver={None, 'keyword', 'comment', 'string'}) (or something like that) and if user says 'comment' - add a card to the returned header. If 'string' return two values: header (without comment) and a separate string value with version information. Then users can do whatever they want with it.

However, I do not like my own suggestion (I do not think it add much practical value). Also, for HST images, we convert PC matrix to CD matrix. Should this kind of modification still be considered WCSLIB? If not, how do I delete the correct comment from the header? Moreover, SIP and table distortion corrections are implemented independently in astropy. So, for best results one should also record the astropy version used to create the header.

[mcara]

Maybe comment could be made more informative, like this:

COMMENT  WCSNAME='MYWCS1' header keyrecords produced by: WCSLIB=7.3, astropy=4.2, stwcs=1.0, ...

The more I think about it, the less I like the idea of adding a new header keyword that, IMO, will contribute to more confusion.

[dhomeier]

I do not think comments should cause technical problems. However, imagine you will get a file with multiple (let's say 4) WCSes in the header and when you inspect that file you will find:

COMMENT  WCS header keyrecords produced by WCSLIB 5.8
COMMENT  WCS header keyrecords produced by WCSLIB 7.1
COMMENT  WCS header keyrecords produced by WCSLIB 7.3

I think several different WCSs created with all different versions or libraries in one file or even HDU are bound to create confusion in the first place. The comments are appended to the end of each header section created by a WCS, so ideally they would stay in that place to be somewhat associated with the corresponding WCS. But I admit that we probably cannot guarantee that order to be preserved though all the potential header updates. Still, I think having the different WCS version with no direct trace of their origin, but say, one implementing the pre-6.3 behaviour of the TPD coefficients, and another the current one, would be even more confusing. For the COMMENT solution I am thinking more of additional information that one can look up in case of trouble, than something to be automatically parsed by readers.

For the alternative of defining a new keyword, this could make the information less ambiguous, but then storing them for each WCS would probably require another class of structured keywords; in short I do not like the effort in defining this in a FITS- and WCS-standard-conforming way either.

In conclusion I probably like your last suggestion to add WCSNAME and info on other creators to the comment best.

[mcara]

@dhomeier Thank you! Could you please rebase and reorder/squash/drop some of the commits in such a way that you have maybe three commits left: adding WCSLIB 7.4, fixing C code, and updating tests?

[dhomeier]

Added a warning now - tested locally with external WCSLIB 7.3.1; I think all CI jobs are set up to build against the bundled version, so they will miss this.
WCSLIB < 7.1 will still produce functional TPD headers, but incompatible with later versions. In principle this might be also fixed in to_header by inserting DOCORR=0 keywords, but that would result in those headers in turn to not be readable with the same WCSLIB version. Deserves probably better coverage in the docs.

[nden]

Was there a patch to wcslib in this PR? If so, I'd like to see it as a separate commit and also document it so next time there is a wcslib release we can follow up. (It should really be fixed in wcslib.)
The question is where to document it. One option is #9018.

[dhomeier]

Was there a patch to wcslib in this PR? If so, I'd like to see it as a separate commit and also document it so next time there is a wcslib release we can follow up. (It should really be fixed in wcslib.)

Patch to wcshdr.c ended up in d0a629b after my squashing attempts.
There is no word yet if it is acknowledged as a bug, since it only affected build with some older compilers.
Nothing I'd consider to constitute a new minor version in any case.

The question is where to document it. One option is #9018.

You mean just record it in an issue? I was thinking of something equivalent to note_sip.rst. Will require a bit more writing work, but I think the relevant parts of 'Historical idiosyncrasies' are documented now in https://www.atnf.csiro.au/people/mcalabre/WCS/wcslib/dis_8h.html and could just be linked.

In the end it perhaps boils down to where to place a recipe like

if 'DQ1.DOCORR' not in header:
    header.insert('CQDIS1', ('DQ1', 'DOCORR:  0'), after=True)

for fixing deprecated headers.

[nden]

wcslib 7.4 was patched before merging.

[dhomeier]

@nden could you hold this for another day? Mark has sent me a mail about the C99 issues, and I'd like to look into fixing the build of cextern/wcslib instead of patching the source.

[nden]

Oops, too late. Can you do it on a separate PR or do you want me to revert it?

[dhomeier]

Indeed, as WCSLIB's configure probably would have worked out on it's own, it just takes the -std=c99 or -std=gnu99 option to let older gcc handle the inside declaration. Tested this with gcc 4.6.2 on an older Linux box, so adding the above to cfg['extra_compile_args'] instead of patching the source code should be the cleaner (and more future-proof) solution.

[nden]

So may be create a PR for this and remove the patch there?

[dhomeier]

Oh, well. Whatever is easier. It's just another one-line change, but the cextern part of d0a629b would be reverted as well. Probably a new PR then.

[pllim]

With this, I now see a new warning 'datfix' made the change: Success pop up in an existing test over at Ginga. I wonder if any of you understand why. Just for academic curiosity. ejeschke/ginga#938

[dhomeier]

I found on possible change to return value of datfix

+       if (status == FIXERR_NO_CHANGE) status = FIXERR_SUCCESS;

Maybe you can bug @mcara about it, who seems to have been mentioned in this context 😄

• In datfix(), ensure that 0 is returned if an informational message
  is set in wcsprm::err. Consequent on feedback independently from
  Mihai Cara and Bruce Merry.

– it looks like in your example there were no changes actually made? This would previously have returned -1, but now it has apparently become indistinguishable from 0 (successful changes); not sure if this was fully intended.

[pllim]

Where do you even report these things upstream to wcslib? Does it have a GitHub?

And thanks for the very useful info, @dhomeier !

[dhomeier]

Ha, I suspect it has been developed and archived safely in a set of tarfiles on Mark Calabretta's disk for a quarter-century!
I think email to Mark is the usual feedback method; I was planning to send him a reply about #11260 (comment) once we have made a decision on #11376, so if you don't want to swamp him I can mention it then.

[mcara]

I found on possible change to return value of datfix

+       if (status == FIXERR_NO_CHANGE) status = FIXERR_SUCCESS;

I am sorry, I probably missed most of the discussion on datfix. What is this change above about? Is it in WCSLIB 7.4, or is this a proposed fix in WCSLIB 7.4? The code above looks weird: replacing one status code with another one. You might as well simply search/replace FIXERR_NO_CHANGE with FIXERR_SUCCESS in WCSLIB but I doubt that was the intention in WCSLIB.

So what is going on @pllim with the datfix message? When and how does it happen?

[pllim]

Thanks, @dhomeier ! (I'll let you bring it up since you already plan to contact him.)

@mcara , the change is already in astropy dev. This new warning caused failures in existing Ginga and SunPy tests.

[dhomeier]

@mcara these are the changes in question, which are part of the entire 7.4 update

astropy/cextern/wcslib/C/wcsfix.c

Line 363 in d971b1c

if (status == FIXERR_NO_CHANGE) status = FIXERR_SUCCESS;

astropy/cextern/wcslib/C/wcsfix.c

Line 586 in d971b1c

if (status == FIXERR_NO_CHANGE) status = FIXERR_SUCCESS;

I am wondering of the motivation was to output the diagnostic message here

astropy/cextern/wcslib/C/wcsfix.c

Lines 181 to 191 in d971b1c

	if (stat[ifix] == FIXERR_NO_CHANGE) {
	// No change => no message.
	wcserr_copy(0x0, info+ifix);
	} else if (stat[ifix] == 0) {
	// Successful translation, but there may be an informative message.
	if (wcs->err && wcs->err->status < 0) {
	wcserr_copy(wcs->err, info+ifix);
	} else {
	wcserr_copy(0x0, info+ifix);
	}

but then, was only the old return value of -1 treated as success, and 0 as warning level?
Seems like one might rather want to check for wcs->err->status < 0 in the stat[ifix] == FIXERR_NO_CHANGE case as well...