WCS: exception if MJD-OBS missing when some other headers present

[bmerry]

Description

If a FITS header has a spectral axis and OBSGEO-* keywords but no MJD-DATE, then _get_components_and_classes crashes because it tries to construct a Time from a NaN MJD. I'm not aware of any requirement in the FITS specification to provide this header, but I didn't comb the entire spec.

Expected behavior

I'm not familiar with the details of SpectralCoord at all so I don't know what the appropriate fix is, but I would expect the resulting WCS to be basically usable.

Actual behavior

When running the test code below, received an exception (see first comment - I forgot to paste it when filing the bug).

Steps to Reproduce

#!/usr/bin/env python3

import astropy.wcs
import astropy.io.fits

header = astropy.io.fits.Header()
header['OBSGEO-X'] = 5109360.133321234
header['OBSGEO-Y'] = 2006852.586042906
header['OBSGEO-Z'] = -3238948.127478876
header['DATE'] = '2020-05-14T08:03:43.930'
header['SPECSYS'] = 'TOPOCENT'
header['TIMESYS'] = 'UTC'
header['EQUINOX'] = 2000.0
header['RADESYS'] = 'FK5'
header['NAXIS'] = 3
header['NAXIS1'] = 100
header['NAXIS2'] = 100
header['NAXIS3'] = 1
header['BUNIT'] = 'Jy/beam'
header['CDELT1'] = 1.0
header['CDELT2'] = 1.0
header['CDELT4'] = 1.0
header['CUNIT1'] = 'deg'
header['CUNIT2'] = 'deg'
header['CUNIT3'] = 'Hz'
header['CTYPE1'] = 'RA---SIN'
header['CTYPE2'] = 'DEC--SIN'
header['CTYPE3'] = 'FREQ'
header['CRVAL1'] = 0.0
header['CRVAL2'] = 0.0
header['CRVAL3'] = 1e9
header['CRPIX1'] = 50.0
header['CRPIX2'] = 50.0
header['CRPIX3'] = 1.0
wcs = astropy.wcs.WCS(header)
wcs.world_axis_object_classes

System Details

Linux-5.4.0-60-generic-x86_64-with-glibc2.29
Python 3.8.5 (default, Jul 28 2020, 12:59:40)
[GCC 9.3.0]
Numpy 1.19.1
astropy 4.3.dev450+gaa71c2f45
Scipy 1.5.2
Matplotlib 3.3.0

[pllim]

I don't see any traceback provided above, but this is what I see with astropy 4.3.dev:

>>> wcs = astropy.wcs.WCS(header)
WARNING: FITSFixedWarning: The WCS transformation has more axes (4) than the image it is associated with (3) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'obsfix' made the change 'Set OBSGEO-L to    21.443889 from OBSGEO-[XYZ].
Set OBSGEO-B to   -30.711055 from OBSGEO-[XYZ].
Set OBSGEO-H to     1083.594 from OBSGEO-[XYZ]'. [astropy.wcs.wcs]

wcs.world_axis_object_classes gives this traceback:

.../astropy/time/core.py in _get_time_fmt(self, val, val2, format, scale, precision, in_subfmt, out_subfmt)
    429             try:
--> 430                 return cls(val, val2, scale, precision, in_subfmt, out_subfmt)
    431             except UnitConversionError:

.../astropy/time/formats.py in __init__(self, val1, val2, scale, precision, in_subfmt, out_subfmt, from_jd)
    148         else:
--> 149             val1, val2 = self._check_val_type(val1, val2)
    150             self.set_jds(val1, val2)

.../astropy/time/formats.py in _check_val_type(self, val1, val2)
    424         if val1.dtype.kind == 'f':
--> 425             val1, val2 = super()._check_val_type(val1, val2)
    426         elif (not orig_val2_is_none

.../astropy/time/formats.py in _check_val_type(self, val1, val2)
    262         if not (ok1 and ok2):
--> 263             raise TypeError('Input values for {} class must be finite doubles'
    264                             .format(self.name))

TypeError: Input values for mjd class must be finite doubles

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
<ipython-input-3-2d4ce764bcb7> in <module>
----> 1 wcs.world_axis_object_classes

.../astropy/wcs/wcsapi/fitswcs.py in world_axis_object_classes(self)
    332     @property
    333     def world_axis_object_classes(self):
--> 334         return self._get_components_and_classes()[1]
    335 
    336     @property

.../astropy/wcs/wcsapi/fitswcs.py in _get_components_and_classes(self)
    420 
    421                 earth_location = EarthLocation(*self.wcs.obsgeo[:3], unit=u.m)
--> 422                 obstime = Time(self.wcs.mjdobs, format='mjd', scale='utc',
    423                                location=earth_location)
    424                 observer_location = SkyCoord(earth_location.get_itrs(obstime=obstime))

.../astropy/time/core.py in __init__(self, val, val2, format, scale, precision, in_subfmt, out_subfmt, location, copy)
   1522                 self._set_scale(scale)
   1523         else:
-> 1524             self._init_from_vals(val, val2, format, scale, copy,
   1525                                  precision, in_subfmt, out_subfmt)
   1526             self.SCALES = TIME_TYPES[self.scale]

.../astropy/time/core.py in _init_from_vals(self, val, val2, format, scale, copy, precision, in_subfmt, out_subfmt)
    370 
    371         # Parse / convert input values into internal jd1, jd2 based on format
--> 372         self._time = self._get_time_fmt(val, val2, format, scale,
    373                                         precision, in_subfmt, out_subfmt)
    374         self._format = self._time.name

.../astropy/time/core.py in _get_time_fmt(self, val, val2, format, scale, precision, in_subfmt, out_subfmt)
    436                 # easier for user to see what is wrong.
    437                 if len(formats) == 1:
--> 438                     raise ValueError(
    439                         f'Input values did not match the format class {format}:'
    440                         + os.linesep

ValueError: Input values did not match the format class mjd:
TypeError: Input values for mjd class must be finite doubles

[bmerry]

Whoops, I pasted the sample code twice instead of pasting the traceback. I get the same exception though. I'll just edit the report to remove the duplicated sample code.

[dhomeier]

Are you getting the same error with DATE-OBS instead of DATE? The standard specifies in 9.2.2 that a time reference point, if present,

must be expressed in one of ISO-8601, JD or MJD

where Table 22 could be understood as requiring an ISO-8601 time to be put in DATE-OBS (or possibly DATE-AVG or another of the observation-related keywords). However the table also only includes MJD-OBS etc., so by the same logic one could expect MJD-DATE to be ignored for this purpose as well.

[bmerry]

Are you getting the same error with DATE-OBS instead of DATE

No, because wcslib's fixers compute MJD-OBS from DATE-OBS. I actually ran into this issue because I was using an observation which has a DATE-OBS but I'd disabled the fixers due to #10365. For now I've just reenabled the fixers.

Section 9.2.2 deals with DATEREF/MJDREF, so it doesn't seem relevant here. Just to be sure that we're talking about the same section numbers, I'm referring to https://fits.gsfc.nasa.gov/standard40/fits_standard40aa-le.pdf.

I see 9.8 says:

One or more of the informational keywords DATE-xxxx and/or MJD-xxxx should be present in all HDUs whenever a meaningful value can be determined.

If I add a DATE-AVG then my header is compliant with that, and the crash still occurs.

[dhomeier]

Yes, I thought the "Time reference" referred to the time point for converting between the different coordinate systems, but that was a mistake. Primarily I was looking for some rule that would accept MJD, but not ISO-8601 formats, but that is not the case, and indeed, as you pointed out, wcslib is already converting all DATE-* keys to the corresponding MJD-* ones in wcs.wcs.mjd*.
But when setting header['MJD-DATE'], I am getting the same exception, so that part at least looks consistent.
If I am not mistaken

astropy/astropy/wcs/wcsapi/fitswcs.py

Lines 421 to 424 in aef31be

	earth_location = EarthLocation(*self.wcs.obsgeo[:3], unit=u.m)
	obstime = Time(self.wcs.mjdobs, format='mjd', scale='utc',
	location=earth_location)
	observer_location = SkyCoord(earth_location.get_itrs(obstime=obstime))

requires a reasonably accurate observation date to compute the observer location relative to a celestial coordinate system, so it seems justified to not silently drop in a generic DATE – that is defined as "Date of HDU creation", so in principle might be delayed by an arbitrary time past the actual observation.
However DATE-AVG/MJD-AVG, and with lower priority -BEG/-END, would certainly acceptable replacements, so _get_components_and_classes should perhaps be modified to search for those in order, if self.wcs.mjdobs has no valid value.

[dhomeier]

Beyond that, I see your use case has a DATE-OBS, so enabling its conversion to MJD is a WCSLIB internal issue, right?
That would leave Astropy with the question if it should really raise an exception if no observation time can be determined –since observer and observer_location were only added in 93076aa, I think one could just proceed with observer = None as already in the case for missing wcs.wcs.obsgeo.

[bmerry]

I think one could just proceed with observer = None as already in the case for missing wcs.wcs.obsgeo.

That is one solution I had in mind. An alternative might be to create the observer without an obstime, but I don't know if that's just going to cause an exception down the road when you try to change frames (possibly it makes sense for some values of SPECSYS?)

However DATE-AVG/MJD-AVG, and with lower priority -BEG/-END, would certainly acceptable replacements

In fact the spec says "In order to compute the velocities it is necessary to have the date and time of the observation; the keywords DATE-AVG and MJD-AVG are reserved for this purpose." That sounds like they should be considered before DATE-OBS/MJD-OBS.

[dhomeier]

Beyond that, I see your use case has a DATE-OBS, so enabling its conversion to MJD is a WCSLIB internal issue, right?

Although _get_components_and_classes could also fall back to wcs.wcs.dateobs or wcs.wcs.dateavg...

[dhomeier]

I think one could just proceed with observer = None as already in the case for missing wcs.wcs.obsgeo.

That is one solution I had in mind. An alternative might be to create the observer without an obstime, but I don't know if that's just going to cause an exception down the road when you try to change frames (possibly it makes sense for some values of SPECSYS?)

observer[_location] is created as SkyCoord, which will assume obstime=J2000.000 as default, so that would almost certainly create some non-sensical results (e.g. it would not trigger the observer.transform_to(ICRS) check somewhere further down, but of course return completely different sky coordinates). And observer is set to None afterwards if it cannot be transformed to ICRS, so we can better leave it at that right away.