Details
-
Type:
Bug
-
Status: Done
-
Resolution: Done
-
Fix Version/s: None
-
Component/s: afw
-
Labels:None
-
Story Points:4
-
Sprint:Arch 2018-11-05, Arch 2018-11-12
-
Team:Architecture
Description
While reviewing DM-16253 I re-processed DECam data through processCCD to test the metadata fix, but found that the reported wcs of the images was occasionally catastrophically wrong. The bug appears seemingly randomly when the fits header is read to create a wcs:
for _i in range(40):
|
wcs_test = butler_test.get('calexp_wcs', {'visit': 411735, 'ccdnum': 10})
|
print(wcs_test.pixelToSky(-.5,-.5))
|
(137.919516, -1.319873)
|
(nan, nan)
|
(137.919516, -1.319873)
|
(138.379127, +87.701332)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(228.379127, -0.000000)
|
(137.919516, -1.319873)
|
(138.379127, +87.701332)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(138.379127, +87.701332)
|
(137.919516, -1.319873)
|
(228.379127, -0.000000)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(nan, nan)
|
(138.379127, +87.701332)
|
(318.379127, -87.701332)
|
(nan, nan)
|
(228.379127, -0.000000)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(137.919516, -1.319873)
|
(228.379127, -0.000000)
|
(137.919516, -1.319873)
|
(228.379127, -0.000000)
|
(228.379127, -0.000000)
|
(nan, nan)
|
(137.919516, -1.319873)
|
The above values are random each time, and the order is not repeatable.
If the full exposure is read in and the wcs is obtained from exposure.getWcs(), then it is correct each time. Also, the above command when run on data ingested in late January 2018 also produces correct results.
Attachments
Issue Links
- relates to
-
DM-16429 Ensure WCS (and other Exposure components) are retrieved properly when loaded individually
- Done
Activity
I have attached a YAML dump of the file header. To make it easier to demonstrate the problem here is some code:
#!/usr/bin/python3
|
import yaml |
import lsst.afw.geom as afwGeom |
|
|
file = "fitsheader-decam-calexp-0412037_10.yaml" |
with open(file) as fd: |
header = yaml.load(fd) |
|
|
for i in range(40): |
wcs = afwGeom.makeSkyWcs(header) |
print(wcs.pixelToSky(-.5, -.5)) |
Running this I get lots of variation. It is the makeSkyWcs constructor that is causing the variation. If you move the constructor out the loop you get repeatable answers. It also seems to me that the first number calculated is always repeatable "(137.917775, -1.316585)" but after that things go random, suggesting some internal state not being reset.
If I use astshim to parse the header and transform the same pixel I get an output value of "137d55m03.728s -1d18m59.4449s". The value I get first from the pixelToSky output is "137d55m03.9907s -1d18m59.7077s" and this is actually what I get from astshim if I request pixel (0.5, 0.5) so the two are consistent modulo a one pixel shift.
Ian Sullivan what value do you get from the exposure WCS for that pixel?
My conclusion is that the instability is coming from somewhere above astshim.
Thanks for the simple example Tim Jenness; that's super useful.
I converted your YAML to a string using fits::makeLimitedFitsHeader, then tried the following:
$ cat astshim_test.py
|
import astshim
|
|
|
f = open("fitsheader-decam-calexp-0412037_10.str", "r")
|
fitshead = f.read()
|
|
|
for _ in range(400):
|
ss = astshim.StringStream(fitshead)
|
fc = astshim.FitsChan(ss, "Encoding=FITS-WCS, IWC=1, SipReplace=0")
|
print(fc.read().show())
|
|
|
$ python astshim_test.py | grep PV1_3 | sort | uniq
|
PV1_3 = 0.00020875199697911739 # Projection parameter 3 for axis 1
|
PV1_3 = 0.0017097517848014832 # Projection parameter 3 for axis 1
|
PV1_3 = 0.0018924465402960777 # Projection parameter 3 for axis 1
|
PV1_3 = -0.0024933274835348129 # Projection parameter 3 for axis 1
|
PV1_3 = 0.007853679358959198 # Projection parameter 3 for axis 1
|
PV1_3 = -0.011195987462997437 # Projection parameter 3 for axis 1
|
PV1_3 = -0.016240030527114868 # Projection parameter 3 for axis 1
|
PV1_3 = -0.019780918955802917 # Projection parameter 3 for axis 1
|
PV1_3 = 0.029940336942672729 # Projection parameter 3 for axis 1
|
PV1_3 = -0.040052175521850586 # Projection parameter 3 for axis 1
|
PV1_3 = -0.7469334602355957 # Projection parameter 3 for axis 1
|
PV1_3 = -0.9999995231628418 # Projection parameter 3 for axis 1
|
PV1_3 = 0.9999995231628418 # Projection parameter 3 for axis 1
|
PV1_3 = 0 # Projection parameter 3 for axis 1
|
PV1_3 = 1.0147895812988281 # Projection parameter 3 for axis 1
|
PV1_3 = 1.0185579796633307e-312 # Projection parameter 3 for axis 1
|
PV1_3 = 1.0397779375729834e-312 # Projection parameter 3 for axis 1
|
PV1_3 = 1.0609978954826362e-312 # Projection parameter 3 for axis 1
|
PV1_3 = 1.0822178533922889e-312 # Projection parameter 3 for axis 1
|
PV1_3 = 1.1458777271212471e-312 # Projection parameter 3 for axis 1
|
PV1_3 = 1.2851807146773715e-70 # Projection parameter 3 for axis 1
|
PV1_3 = 1.2862567834555961e+248 # Projection parameter 3 for axis 1
|
PV1_3 = 1 # Projection parameter 3 for axis 1
|
PV1_3 = -2.2954349517822266 # Projection parameter 3 for axis 1
|
PV1_3 = 2.3364124312642001e-307 # Projection parameter 3 for axis 1
|
PV1_3 = 2.3364158264574656e-307 # Projection parameter 3 for axis 1
|
PV1_3 = 2.3364616615665505e-307 # Projection parameter 3 for axis 1
|
PV1_3 = 2.5691505243038807e+161 # Projection parameter 3 for axis 1
|
PV1_3 = 2.6204526022630148e-310 # Projection parameter 3 for axis 1
|
PV1_3 = 4.6548099670614214e-310 # Projection parameter 3 for axis 1
|
PV1_3 = 4.8909910913784047e+199 # Projection parameter 3 for axis 1
|
PV1_3 = 57073.03125 # Projection parameter 3 for axis 1
|
PV1_3 = 6.0134693029269245e-154 # Projection parameter 3 for axis 1
|
PV1_3 = 6306.33203125 # Projection parameter 3 for axis 1
|
PV1_3 = 6.7786571791757129e+78 # Projection parameter 3 for axis 1
|
PV1_3 = 6.9359554423490891e-310 # Projection parameter 3 for axis 1
|
PV1_3 = 7.4057653104688004e-312 # Projection parameter 3 for axis 1
|
PV1_3 = 8.4879831638610893e-314 # Projection parameter 3 for axis 1
|
PV1_3 = -nan # Projection parameter 3 for axis 1
|
(Apologies for my naivety in extracting values from astshim; I've never attempted to use this API before, and I didn't want to spend long understanding it.)
The above is effectively replicating what's happening in the guts of afw, and I'm pretty sure that's where the problem is arising: every time the WCS is incorrect, its because FitsChan.read() is giving back something with PV1_3 or PV2_3 set to garbage.
Based on the above, I think I disagree with Tim Jenness's conclusions — this looks like it is happening in astshim or maybe starlink_ask.
I note that Valgrind says a bunch of things like:
==153626== Use of uninitialised value of size 8
|
==153626== at 0x5A947EB: __cos_avx (s_sin.c:510)
|
==153626== by 0x5A53F9D: sincos (s_sincos.c:43)
|
==153626== by 0x9C4E9975: astEraS2c (s2c.c:21)
|
==153626== by 0x9C18A2B3: Transform (sphmap.c:1282)
|
==153626== by 0x9BFFEB9B: astTransform_ (mapping.c:24060)
|
==153626== by 0x9BD74724: Transform (cmpmap.c:3790)
|
==153626== by 0x9BFFEB9B: astTransform_ (mapping.c:24060)
|
==153626== by 0x9BD7473D: Transform (cmpmap.c:3791)
|
==153626== by 0x9BFFEB9B: astTransform_ (mapping.c:24060)
|
==153626== by 0x9BD7473D: Transform (cmpmap.c:3791)
|
==153626== by 0x9BFFEB9B: astTransform_ (mapping.c:24060)
|
==153626== by 0x9BD7473D: Transform (cmpmap.c:3791)
|
which look suspicious, but I'm not sufficiently familiar with this part of the codebase to dive in further without significant effort.
Russell Owen or Tim Jenness — do you have time to look at this further? I realise you're both busy, but this looks like a major problem and you have the most expertise here.
Thanks John Swinbank, that's interesting. Here's my code for completeness:
#!/usr/bin/python3
|
import yaml |
import numpy as np |
from astropy.coordinates import Angle |
import astshim as ast |
import lsst.afw.geom as afwGeom |
|
|
file = "tests/data/fitsheader-decam-calexp-0412037_10.yaml" |
with open(file) as fd: |
header = yaml.load(fd) |
|
|
point = [-0.5, -0.5] |
for i in range(40): |
wcs = afwGeom.makeSkyWcs(header) |
ra, dec = wcs.pixelToSky(*point) |
print(Angle(ra.asDegrees(), unit="deg").to_string("deg"), |
Angle(dec.asDegrees(), unit="deg").to_string("deg")) |
|
|
# Now do the conversion natively using AST
|
header = header.toOrderedDict() |
|
|
|
|
def makeWcs(header): |
fc = ast.FitsChan(ast.StringStream()) |
for k in header: |
v = header[k] |
comment = "" |
if isinstance(v, int): |
fc.setFitsI(k, v, comment)
|
elif isinstance(v, float): |
fc.setFitsF(k, v, comment)
|
elif isinstance(v, str): |
fc.setFitsS(k, v, comment)
|
else: |
# print(f"Skipping key {k}") |
pass |
fc.setCard(0) |
wcs = fc.read() |
fc.emptyFits()
|
return wcs |
|
|
point = [0.5, 0.5] |
input = np.array([ |
*point |
])
|
for i in range(40): |
wcs2 = makeWcs(header) |
result = wcs2.applyForward(input) |
print(Angle(result[0], unit="rad").to_string("deg"), |
Angle(result[1], unit="rad").to_string("deg")) |
(it has both versions in it)
If I use makeLimitedFitsHeader I do get the instability. I see that the FITS header coming back from that call does not handle booleans properly but fixing that doesn't make a difference. If I build the FitsChan myself one line at a time everything works but if we build a FITS header and pass that in directly things go wonky. I'm going to have to try this in native AST to see if it's astshim that's the problem.
I can confirm that I see the problem in pyast and it does work with astSetFitsX but not when reading from the header as a single text string. I have sent my test code to David Berry.
Thank you for chasing this Tim Jenness!
The problem is that makeLimitedFitsHeader writes some of the QV parameters as integer 0 rather than floating point 0.0. AST then gets a bit confused and undefined behavior occurs. I have reported this to David Berry but I will try to fix afw so we do it correctly at our end. The reason it worked for me in using astSetFitsX was that I was handling the data types properly from the property list.
I hadn't realized calexp_wcs was using the header, and that itself is a problem we need to fix.
Russell Owen would you be able to review this? It is mostly additions to the tests and a couple of minor changes to the format strings to ensure a decimal point is written out.
Pull Request: https://github.com/lsst/afw/pull/414/files
David Berry has also made a change to AST upstream for next time we sync up our version.
It would be nice in the future if we could have a helper function somewhere that would directly convert a PropertyList to a FitsChan without going through the string phase. This would be relatively easy to write (I have the bare bones of it in Python on this ticket) although I'm not sure where it would go.
I hadn't realized calexp_wcs was using the header, and that itself is a problem we need to fix.
Agreed (do we have a ticket for this?). But I'm also worried about the complexity of metadata flow through the system: it's not obvious, and I don't know where it's documented, that looking at the header is not an appropriate way to get the WCS. It seems like we need a “metadata model” to go with the science data model.
Regardless, thank you very much Tim Jenness for taking care of this ticket!
Thanks Tim Jenness for fixing this.
Adding another vote to the "calexp_wcs should return the true wcs, always". I too was unaware that that was the behavior. I don't think it's affected jointcal's output, but I use the magic metadata getters exclusively there, so I don't have to read in the full exposure object for a big speedup. Someone who can explain it in more detail, please file a Critical ticket?
Per the discussion thread on Slack here, it sounded to me like calexp_wcs was designed to use the header information in order to bypass reading in an entire Exposure just to get at the wcs from a single FITS extension. I'm not sure if that's a design we want to uphold or not, but we should at least be aware of that before we go changing how calexp_wcs works.
Also wanted to chime in and say thank you everyone (and especially Tim!) for getting to the bottom of this. It explains why I had a few CCDs showing up in very wrong, nonsensical locations... but only sometimes.
Just created a ticket for the deeper problem: DM-16429, and added some notes there about how we should go about implementing it; it should definitely be possible to avoid the potential inefficiency Meredith Rawls just noted.
The demo is currently failing on this branch because some of the declinations in the catalog have shifted by a tiny amount (suggesting that we are putting coordinates in catalogs derived at some point from WCS in headers and not the real WCS).
This looks great. Two small suggestions on github.
Russell Owen this evening I have had a go at using the PropertyList->FitsChan conversion directly without going through the string. It seems to work fine but it will need another review. If people prefer I merge the reviewed code and then put this improvement on a separate ticket please let me know.
(but I must remember to fix the demo – I assume I just copy in the new versions of the catalogs – it still surprises me that they are changing at all).
Merged. Jenkins passed.
My understanding is that Exposure.getWcs() will use the full serialized WCS retrieved from a table extension, and that butler calexp_wcs will try to calculate an approximate WCS from the FITS headers. One difference between that January files and the modern files is that the new files have distortion terms included as QV headers whereas no distortion was present previously. I'm a bit surprised that QV is being used rather than PV but maybe AST is using that as a default somewhere. My guess would be that Starlink AST is hitting some edge cases with the distortion handling.