Fix and extend `wcs1d-fits` loader for multi-D WCS #1009

dhomeier · 2023-02-14T02:05:58Z

To address #592 (comment); quick fix for enabling writing of spectral cubes with celestial components in WCS.
Reading in with the full WCS is now supported as well, plus adding mask and uncertainty to the attributes read from file.
The latter at this point are still limited to the format storing them all in different ImageHDUs; an alternative common format storing the fields as sub-arrays of an N+1 image yet to be implemented.
Fixes #592

codecov · 2023-02-14T02:13:35Z

Codecov Report

Merging #1009 (1d5a59f) into main (8179313) will increase coverage by 0.24%.
The diff coverage is 84.88%.

@@            Coverage Diff             @@
##             main    #1009      +/-   ##
==========================================
+ Coverage   70.01%   70.26%   +0.24%     
==========================================
  Files          64       64              
  Lines        4346     4419      +73     
==========================================
+ Hits         3043     3105      +62     
- Misses       1303     1314      +11

Impacted Files	Coverage Δ
specutils/io/default_loaders/wcs_fits.py	`45.71% <84.88%> (+10.75%)`	⬆️

... and 3 files with indirect coverage changes

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

pllim · 2023-02-14T12:43:04Z

So this fixes #592 ?

dhomeier · 2023-02-14T15:00:45Z

This should work for your example.
Reading the full WCS in rather than only the spectral dimension may of course have consequences if spectrum.wcs is used at a later point. Another immediate change is that this will now only accept the official FITS spectral axis keywords recognised my WCSLIB, FREQ, AFRQ, ENER, WAVN,VRAD, WAVE, VOPT, ZOPT, AWAV, VELO, BETA.
The tests are adapted, but if downstream applications have used names like wavelength or wavenumber, this may need additional pre-processing of the header.

dhomeier · 2023-02-14T16:18:51Z

I was somewhat surprised when updating the writer, that there are apparently WCS that have a spectral component, yet
wcs.spectral.all_pix2world(pix) returns a different (wrong) result than wcs.all_pix2world(pix) (the former fails on all other tests that only input a 1D, pure spectral WCS). But the check should probably be there anyway, if it ever has to deal with non-FITS WCS.

nmearl · 2023-02-15T19:48:25Z

specutils/tests/test_loaders.py

@@ -707,11 +707,10 @@ def test_tabular_fits_compressed(compress, tmpdir):


 @pytest.mark.parametrize("spectral_axis",
-                         ['wavelength', 'frequency', 'energy', 'wavenumber'])
+                         ['WAVE', 'FREQ', 'ENER', 'WAVN'])


Just some minor context: if I recall, the spectral axis names were meant to correspond with the class properties and not the FITS naming standard. Either way, it's not an issue, but figured the history of the usage might add some context.

dhomeier · 2023-02-15T20:49:17Z

I have added read and write for spectrum.mask and given the HDUs names, handling it I think more similar to how the cubeviz loader treats its files. Things like detection strategy for the spectrum HDU (and mask? Should it be found under other names as well?) may need some more discussion, also perhaps more test coverage for the naming schemes.
In principle this could probably extended in a similar way to also handle uncertainties.

dhomeier · 2023-02-15T20:52:08Z

Apparently something in oldestdeps is pulling numpy 1.24.2 in, overriding numpy==1.19 and breaking the old scipy version. I guess that means 1.19 is not really supported anymore, but don't know what it is that's requiring a newer version.

pllim · 2023-02-15T21:08:00Z

Apparently something in oldestdeps is pulling numpy 1.24.2

Hmm easy way out is to bump min version of numpy though even that might not work. Another thing you can do is to hunt down what is overwriting the numpy pin and pin that package to an older version as well.

extended in a similar way to also handle uncertainties

Definitely needed if you want it to be generic enough since Spectrum1D does have an uncertainty component. I just didn't do it downstream because apparently model fitting in Cubeviz does not generate that.

dhomeier · 2023-03-03T15:24:49Z

Marking this as ready for review now, as the multi-HDU format for masks and uncertainties is functional (unless there are more labels the loader needs to identify).
The format reading different attributes from data slices in a single HDU, which would work similarly to 6dfgs

specutils/specutils/io/default_loaders/sixdfgs_reader.py

Lines 188 to 189 in 6eb7f96

    
           flux = hdu.data[0] * Unit("count") / w.wcs.cunit[0] 
        
           uncertainty = VarianceUncertainty(hdu.data[1])

could still be added ~~here or~~ in a follow-up PR.

specutils/io/default_loaders/wcs_fits.py

pllim · 2023-03-08T13:57:06Z

specutils/io/default_loaders/wcs_fits.py

-    verbose : bool
+    hdu : int, str or None, optional
+        The index or name of the HDU to load into this spectrum
+        (default: find 1st applicable `ImageHDU`).


This confuses me. Is this loader only for image? A spectrum can also come as table HDU, right?

wcs1d handles only spectra with flux represented as image and spectral axis defined by the WCS; tables are covered by tabular-fits. That already supports masks and uncertainties, more or less, which are simply additional table columns in that representation.

specutils/io/default_loaders/wcs_fits.py

pllim · 2023-03-08T14:19:02Z

specutils/io/default_loaders/wcs_fits.py

+    uncertainty_name : str or `None`, optional
+        HDU name to store uncertainty under (default set from type; `None`: do not save)
+    unit : str or :class:`~astropy.units.Unit`, optional
+        Unit for the flux (and associated uncertainty; defaults to `spectrum.flux.unit`)


Unless spectrum.flux.unit actually generates API doc, use double-backticks. I cannot check because RTD build failed on your PR. Also see #1028.

spectrum.flux does, but not specifically unit, is that what you mean? I don't know if the loaders themselves generate API docs, or rather where they might be linked from ;-)
Thanks for pointing out the RTD failure; throws in all current PRs and does not seem related – from a quick glance some parsing error related to reproject, but needs further investigation – oh, it already does!

Actually it seems none of the functions in default_loaders have any docs built at all, though it would certainly be useful to have them and link them in https://specutils.readthedocs.io/en/latest/spectrum1d.html#list-of-loaders.

pllim · 2023-03-15T17:42:02Z

Sorry, I think my fixes to the tests caused a conflict here, so you need to rebase and resolve the conflict.

dhomeier · 2023-03-15T18:07:27Z

Doh, that seems to have screwed up more commits than I even had in this branch!

specutils/io/default_loaders/wcs_fits.py

specutils/tests/test_loaders.py

pllim · 2023-03-15T19:02:03Z

specutils/tests/test_loaders.py

@@ -875,7 +1026,7 @@ def test_wcs1d_fits_compressed(compress, tmp_path):
    # Try again without compression suffix:
    with warnings.catch_warnings():
        warnings.simplefilter('ignore', FITSFixedWarning)
-        os.system(f'mv {tmpfile}{ext[compress]} {tmpfile}')
+        os.system(f'mv {tmpfile}.{ext[compress]} {tmpfile}')


Side note: Why call os.system instead of shutil? Anyways, I guess fixing this is out of scope.

Perhaps because I (or whoever wrote the tests originally) did not know any better at the time, or perhaps because os.system is needed anyway to do the compression – or is there also some kind of exec functionality with shutil?

I don't see how "mv" cares about compression? https://docs.python.org/3/library/shutil.html#shutil.move

Not at all, I guess 😄; just meant that it's kind of the direct counterpart to
os.system(f'{compress} {tmpfile}') just above, and that cannot be replaced with shutil.

Anyway there are a couple of those compress-rename duos, so refactoring them would be better left to a separate PR.

pllim

I left what general review comments I could. Now it is up to specutils maintainers to review. Thanks!

rosteen · 2023-03-16T19:37:03Z

specutils/io/default_loaders/wcs_fits.py

@@ -157,12 +246,16 @@ def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, **kwargs)
        raise ValueError(f'Only Spectrum1D objects with valid WCS can be written as wcs1d: {err}')


I'm not that familiar with the workings of WCS (and especially wcs.to_fits()) - is there any realistic scenario for setting hdu to a value other than 0, given that the initial hdulist is always created from the spectrum's WCS anyway? I guess what I'm asking is, is there really utility in having empty HDUs in front of the specified hdu index? There are no options to specify the HDU numbers for uncertainty and mask, so it's not like the user could write out in order [uncertainty, mask, flux], for example. Maybe we just get rid of the hdu argument altogether?

As long as this is not considered "public" API and ok to break just like that, I agree.

I'm not that familiar with the workings of WCS (and especially wcs.to_fits()) - is there any realistic scenario for setting hdu to a value other than 0, given that the initial hdulist is always created from the spectrum's WCS anyway?

It's not as much associated to using a WCS, but rather to the flux being stored in an ImageHDU, which can be the PrimaryHDU or an extension, while tables can only be stored in extensions. However for symmetry reasons or whatever sometimes writing the image to an extension, too, and using just a basic PrimaryHDU, is desired (there was one comment somewhere sketching just this scenario, but I cannot find it either in the Slack discussion or any related issue now).
For the record: that example is right in #592 (comment) that set this whole thing off.
I think originally accepting any value for hdu was just introduced as a poor substitute for not being able to write multiple spectra to one file/hdulist, but there is probably no real use case for this. For hdu=1 I'd say there is, and this should still remain available.

Thanks for the explanation...I guess it makes sense to leave it in then for now.

specutils/io/default_loaders/wcs_fits.py

rosteen · 2023-03-16T20:34:51Z

Thanks @dhomeier, this is a great improvement. I had one question and left one suggestion, and there's one suggestion from @pllim remaining, but other than that I think this looks good. I appreciate the test coverage.

I was going to say, maybe we change the name from wcs1d_fits to just wcs_fits if it handles multi-D WCS now, but perhaps we go with the standard answer of "1D means one spectral dimension" like in Spectrum1D.

dhomeier · 2023-03-17T12:13:42Z

I was going to say, maybe we change the name from wcs1d_fits to just wcs_fits if it handles multi-D WCS now, but perhaps we go with the standard answer of "1D means one spectral dimension" like in Spectrum1D.

I had been thinking about that some times; to me the distinctive property was actually usually having the flux as Image data rather than a table column, so having something like image-fits as counterpart to tabular-fits seemed more sensible.
But I also think the "1D spectral" designation is sensible enough, and changing it is probably not worth deprecating the old label.

rosteen · 2023-03-17T18:06:41Z

@dhomeier Can I commit the last suggestions for you and merge this?

rosteen

Approved, need to commit last two suggestions before merging.

Co-authored-by: P. L. Lim <2090236+pllim@users.noreply.github.com>

pllim · 2023-03-17T18:24:48Z

Hmm does this mean we should ditch Jdaviz custom writer and use this?

dhomeier · 2023-03-17T18:28:23Z

@dhomeier Can I commit the last suggestions for you and merge this?

I've included them in some more local updates I had made; should be good to merge then.

rosteen · 2023-03-17T18:33:03Z

@dhomeier Can I commit the last suggestions for you and merge this?

I've included them in some more local updates I had made; should be good to merge then.

Blast, I went ahead and committed and merged, thinking you might be done and logged off for the day. Are your other updates worth opening a follow-up PR for? Sorry about that.

dhomeier · 2023-03-17T18:37:01Z

Ah, remote tests are taking too long here (spPlate really bogs everything down).
Will think about the follow-up, the extra test might be good to have, but now I am going to be offline for a while.

pllim · 2023-03-17T21:35:49Z

Thanks, @dhomeier ! I am trying out this improved format over at spacetelescope/jdaviz#2094 but I encountered error. Perhaps we need an extra step to convert bool mask to uint16 (I guess it could be uint8 but uint16 is common in bit-plane format)?

>       fitted_spectrum.write(out_fn, format="wcs1d-fits", overwrite=True)

jdaviz/configs/default/plugins/model_fitting/tests/test_fitting.py:229: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
astropy/nddata/mixins/ndio.py:103: in __call__
    self.registry.write(self._instance, *args, **kwargs)
astropy/io/registry/core.py:383: in write
    return writer(data, *args, **kwargs)
specutils/io/default_loaders/wcs_fits.py:278: in wcs1d_fits_writer
    hdulist[-1].data = spectrum.mask
astropy/utils/decorators.py:848: in __set__
    ret = self.fset(obj, val)
...
E           KeyError: 'bool'

https://github.com/spacetelescope/jdaviz/actions/runs/4451633462/jobs/7818514722?pr=2094

dhomeier added bug io spectral-cube labels Feb 14, 2023

dhomeier mentioned this pull request Feb 14, 2023

Writer needed for cube data #592

Closed

pllim mentioned this pull request Feb 14, 2023

Add cube writer to properly write out cube generated by model fitting spacetelescope/jdaviz#2012

Merged

12 tasks

nmearl reviewed Feb 15, 2023

View reviewed changes

dhomeier force-pushed the wcs1d-3d branch 3 times, most recently from de33ac4 to 20f9de7 Compare February 18, 2023 23:08

pllim mentioned this pull request Feb 22, 2023

Is it possible to write out 3D Spectrum1D and have it roundtrip? #951

Closed

dhomeier mentioned this pull request Feb 23, 2023

Ensure masks round-trip from Spectrum1D (or other nddata-derived objects) glue-viz/glue-astronomy#23

Open

This was referenced Feb 23, 2023

Revamp Spectrum1D tabular-fits writer #905

Open

jdaviz-cube writer should also write out uncertainty spacetelescope/jdaviz#2044

Open

dhomeier changed the title ~~Fix wcs1d-fits loader for multi-D WCS~~ Fix and extendwcs1d-fits loader for multi-D WCS Mar 3, 2023

dhomeier marked this pull request as ready for review March 3, 2023 15:19

dhomeier force-pushed the wcs1d-3d branch from 90daabd to 5f79132 Compare March 3, 2023 16:02

pllim changed the title ~~Fix and extendwcs1d-fits loader for multi-D WCS~~ Fix and extend wcs1d-fits loader for multi-D WCS Mar 3, 2023

pllim reviewed Mar 8, 2023

View reviewed changes

specutils/io/default_loaders/wcs_fits.py Outdated Show resolved Hide resolved

pllim reviewed Mar 8, 2023

View reviewed changes

specutils/io/default_loaders/wcs_fits.py Show resolved Hide resolved

pllim reviewed Mar 8, 2023

View reviewed changes

specutils/io/default_loaders/wcs_fits.py Show resolved Hide resolved

pllim reviewed Mar 8, 2023

View reviewed changes

dhomeier added 3 commits March 15, 2023 18:47

WIP: fix wcs1d-fits loader for celestial WCS

8d788d0

DNM: read in full WCS from wcs1d-fits

d98d48a

Support HDU naming, add mask HDU r/w for wcs1d-fits

6f4908d

dhomeier force-pushed the wcs1d-3d branch 2 times, most recently from d5705fa to 696a0b3 Compare March 15, 2023 18:06

Handle custom uncertainty names and test them

abb8c1e

dhomeier force-pushed the wcs1d-3d branch from 696a0b3 to abb8c1e Compare March 15, 2023 18:24

pllim reviewed Mar 15, 2023

View reviewed changes

specutils/io/default_loaders/wcs_fits.py Outdated Show resolved Hide resolved

pllim reviewed Mar 15, 2023

View reviewed changes

specutils/tests/test_loaders.py Outdated Show resolved Hide resolved

pllim reviewed Mar 15, 2023

View reviewed changes

dhomeier mentioned this pull request Mar 15, 2023

For tabular_fits output, force uncertainty to StdDev #1027

Merged

rosteen reviewed Mar 16, 2023

View reviewed changes

specutils/io/default_loaders/wcs_fits.py Outdated Show resolved Hide resolved

rosteen approved these changes Mar 17, 2023

View reviewed changes

Apply suggestions from code review

1d5a59f

Co-authored-by: P. L. Lim <2090236+pllim@users.noreply.github.com>

rosteen merged commit f5db745 into astropy:main Mar 17, 2023

pllim mentioned this pull request Mar 17, 2023

Cubeviz to use wcs1d-fits from specutils to write out fitted cube spacetelescope/jdaviz#2094

Draft

13 tasks

dhomeier mentioned this pull request Mar 29, 2023

Support write and read of boolean masks in wcs1d-fits #1051

Merged

dhomeier mentioned this pull request Jun 5, 2023

Cycle 3 Funding: Homeier | General Maintainer astropy/astropy-project#317

Open

		@@ -157,12 +246,16 @@ def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, **kwargs)
		raise ValueError(f'Only Spectrum1D objects with valid WCS can be written as wcs1d: {err}')

Fix and extend wcs1d-fits loader for multi-D WCS #1009

Fix and extend wcs1d-fits loader for multi-D WCS #1009

Conversation

dhomeier commented Feb 14, 2023 • edited

codecov bot commented Feb 14, 2023 • edited

Codecov Report

pllim commented Feb 14, 2023

dhomeier commented Feb 14, 2023

dhomeier commented Feb 14, 2023

Choose a reason for hiding this comment

dhomeier commented Feb 15, 2023

dhomeier commented Feb 15, 2023

pllim commented Feb 15, 2023

dhomeier commented Mar 3, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dhomeier Mar 8, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pllim commented Mar 15, 2023

dhomeier commented Mar 15, 2023

Choose a reason for hiding this comment

dhomeier Mar 17, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dhomeier Mar 17, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pllim left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dhomeier Mar 17, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rosteen commented Mar 16, 2023

dhomeier commented Mar 17, 2023

rosteen commented Mar 17, 2023

rosteen left a comment

Choose a reason for hiding this comment

pllim commented Mar 17, 2023

dhomeier commented Mar 17, 2023

rosteen commented Mar 17, 2023

dhomeier commented Mar 17, 2023

pllim commented Mar 17, 2023

Fix and extend `wcs1d-fits` loader for multi-D WCS #1009

Fix and extend `wcs1d-fits` loader for multi-D WCS #1009

dhomeier commented Feb 14, 2023 •

edited

codecov bot commented Feb 14, 2023 •

edited

dhomeier commented Mar 3, 2023 •

edited

dhomeier Mar 8, 2023 •

edited

dhomeier Mar 17, 2023 •

edited

dhomeier Mar 17, 2023 •

edited

dhomeier Mar 17, 2023 •

edited