Update 'wcs_autolinking' code to handle N-D cases

[kakirastern]

Description

To rewrite the __init__ of the WCSLink class to handle the case of N-D data cubes where N >= 4. While retaining the linking of the wcs spatial coordinates, support has been added for other quantities as dimensions as well including but not limited to time and wavelength.

[codecov]

Codecov Report

Merging #2161 into master will increase coverage by 0.04%.
The diff coverage is 98.07%.

@@            Coverage Diff             @@
##           master    #2161      +/-   ##
==========================================
+ Coverage   87.85%   87.89%   +0.04%     
==========================================
  Files         246      246              
  Lines       22730    22764      +34     
==========================================
+ Hits        19969    20009      +40     
+ Misses       2761     2755       -6     

Impacted Files	Coverage Δ
glue/plugins/wcs_autolinking/wcs_autolinking.py	97.54% <98.07%> (+6.63%)	⬆️
glue/plugins/dendro_viewer/data_factory.py	100.00% <0.00%> (+2.27%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b2ca7aa...222d31e. Read the comment docs.

[astrofrog]

Thanks for working on this! Before I review this further, my main comment at the moment is that it would be good to make it so that the celestial axes are optional - that is if you had e.g. a time dataset and a time + spectral dataset, the time axes would match. Or if you had one dataset with celestial and spectral axes, and a 1D spectrum, the 1D spectrum would match.

Another big picture comment - a spectral axis can have different physical types, so we need to think about how to match those (e.g. if you have a frequency and a wavelength dataset). Maybe the matching could be done based on high level classes rather than or in addition to physical types?

Finally, it would be great to add tests for this, including for the use cases I mentioned above. To be clear, this is already an improvement on what we have, so if there is no time to improve things along the lines mentioned above, we could always merge this and then improve later - but I think it would be good to already add tests for a number of cases even if we have to xfail them for now.

[kakirastern]

Hi @astrofrog I think I would prefer to have this PR merged as soon as possible if that's okay with you, and then for me to open another PR to extend it to include the most generalized use case where we have a data cube having unequal number of axes and no spatial dimensions, but with say temporal, spectral, and maybe even radial velocity/redshift/cz as the third component plus etc. (I.e., making the spatial axes optional without loss of generality.)

[astrofrog]

Thanks for the updates and tests! See below for a few more tests to add (but you can xfail the ones that don't work)

I also wonder if we can improve the testing to avoid repetition and just be able to assert in a single line which dimensions are linked - but that can wait to a separate PR 😄

[kakirastern]

Hi @astrofrog

Actually I took the time to modify the code to make the celestial axes optional. But does that now mean the following tests is no longer needed?

def test_wcs_autolink_dimensional_mismatch():
    # No links should be found because the WCS don't actually have well defined
    # physical types.
    wcs1 = WCS(naxis=1)
    wcs1.wcs.ctype = ['FREQ']
    wcs1.wcs.set()
    data1 = Data()
    data1.coords = wcs1
    data1['x'] = [1, 2, 3]
    wcs2 = WCS(naxis=3)
    wcs2.wcs.ctype = 'DEC--TAN', 'FREQ', 'RA---TAN'
    wcs2.wcs.set()
    data2 = Data()
    data2.coords = wcs2
    data2['x'] = np.ones((2, 3, 4))
    dc = DataCollection([data1, data2])
    links = wcs_autolink(dc)
    assert len(links) == 0

I am asking because this test is the only one failing after the modifications have been made.

[kakirastern]

Another big picture comment - a spectral axis can have different physical types, so we need to think about how to match those (e.g. if you have a frequency and a wavelength dataset). Maybe the matching could be done based on high level classes rather than or in addition to physical types?

This will need a lot more work and is potentially outside the scope of GSoC... But we could come back and work on this later if that is possible.

[kakirastern]

PR is ready for another review

[astrofrog]

This is almost ready! Just a few comments below - thanks!

[kakirastern]

Am seeing the following in the error log entry for the pyqt513 test:

ERROR: THESE PACKAGES DO NOT MATCH THE HASHES FROM THE REQUIREMENTS FILE. If you have updated the package versions, please update the hashes. Otherwise, examine the package contents carefully; someone may have tampered with them.
    PyQt5==5.13.* ...

[astrofrog]

Just a couple of comments below - thanks!

[kakirastern]

Could you take a quick look again and let me know what errors you get if you do just use the default forwards and backwards?

So I have tried reverting and the error is gone now... But I bumped into another issue, which is some missing edges in the connection graph:

Think some of the changes we have introduced after a previous code review is causing this. Will need to back track and debug soon.

[kakirastern]

And one of the CI tests are failing now...

================================================================== short test summary info ==================================================================
FAILED glue/plugins/wcs_autolinking/tests/test_wcs_autolinking.py::test_2d_and_1d_data_cubes_with_no_celestial_axes - assert 0 == 1
============================================== 1 failed, 11 passed, 1 xfailed, 1 xpassed, 13 warnings in 1.05s ==============================================

Let me sort these issues out and get back to you shortly.

[kakirastern]

So I am finally able to fetch the issue output (as it is not technically an error) which is causing the missing edges between some of the nodes:

em.wl matches with em.wl
custom:pos.helioprojective.lat matches with custom:pos.helioprojective.lat
custom:pos.helioprojective.lon matches with custom:pos.helioprojective.lon
time matches with time
slicing_axes1: [3, 2, 1, 0]
slicing_axes2: [3, 2, 1, 0]
slices1: [slice(None, None, None), slice(None, None, None), slice(None, None, None), slice(None, None, None)]
slices2: [slice(None, None, None), slice(None, None, None), slice(None, None, None), slice(None, None, None)]
cid1_slices: [Pixel Axis 3, Pixel Axis 2, Pixel Axis 1, Pixel Axis 0]
cid2_slices: [Pixel Axis 3, Pixel Axis 2, Pixel Axis 1, Pixel Axis 0]
None None None None
self._physical_types_1: ['custom:pos.helioprojective.lat', 'custom:pos.helioprojective.lon', 'em.wl', 'time']
self._physical_types_2: ['custom:pos.helioprojective.lat', 'custom:pos.helioprojective.lon', 'em.wl', 'time']

The None None None None line of code above as output of get_cids_and_functions is where the problem lies.

[kakirastern]

This only happens if I have 4 matching axes, but not when 3 axes are matching:

custom:pos.helioprojective.lon matches with custom:pos.helioprojective.lon
custom:pos.helioprojective.lat matches with custom:pos.helioprojective.lat
time matches with time
slicing_axes1: [2, 1, 0]
slicing_axes2: [2, 1, 0]
slices1: [slice(None, None, None), slice(None, None, None), slice(None, None, None)]
slices2: [slice(None, None, None), slice(None, None, None), slice(None, None, None), 0]
cid1_slices: [Pixel Axis 2 [x], Pixel Axis 1 [y], Pixel Axis 0 [z]]
cid2_slices: [Pixel Axis 2, Pixel Axis 1, Pixel Axis 0]
pixel_cids1: [Pixel Axis 2 [x], Pixel Axis 1 [y], Pixel Axis 0 [z]]
WARNING: ErfaWarning: ERFA function "utctai" yielded 1 of "dubious year (Note 3)" [astropy._erfa.core]
WARNING:sunpy:ErfaWarning: ERFA function "utctai" yielded 1 of "dubious year (Note 3)"
WARNING: ErfaWarning: ERFA function "taiutc" yielded 1 of "dubious year (Note 4)" [astropy._erfa.core]
WARNING:sunpy:ErfaWarning: ERFA function "taiutc" yielded 1 of "dubious year (Note 4)"
[Pixel Axis 2 [x], Pixel Axis 1 [y], Pixel Axis 0 [z]] [Pixel Axis 2, Pixel Axis 1, Pixel Axis 0] <function get_cids_and_functions.<locals>.forwards at 0x7f8ffea61940> <function get_cids_and_functions.<locals>.backwards at 0x7f8f3f071310>
self._physical_types_1: ['custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat', 'time']
self._physical_types_2: ['custom:pos.helioprojective.lat', 'custom:pos.helioprojective.lon', 'time']

[kakirastern]

I was able to discover the underlying error regarding the get_cids_and_functions I mentioned previously:

Traceback (most recent call last):
  File "/Users/krisstern/glue/glue/utils/misc.py", line 55, in result
    return func(*args, **kwargs)
  File "/Users/krisstern/glue/glue/app/qt/application.py", line 257, in _choose_load_data
    app.add_datasets(data)
  File "/Users/krisstern/glue/glue/app/qt/application.py", line 1412, in add_datasets
    run_autolinker(self.data_collection)
  File "/Users/krisstern/glue/glue/dialogs/autolinker/qt/autolinker.py", line 89, in run_autolinker
    suggestions = find_possible_links(data_collection)
  File "/Users/krisstern/glue/glue/core/autolinking.py", line 26, in find_possible_links
    links = function(data_collection)
  File "/Users/krisstern/glue/glue/plugins/wcs_autolinking/wcs_autolinking.py", line 223, in wcs_autolink
    link = WCSLink(data1, data2)
  File "/Users/krisstern/glue/glue/plugins/wcs_autolinking/wcs_autolinking.py", line 157, in __init__
    pixel_cids1, pixel_cids2, forwards, backwards = get_cids_and_functions(
  File "/Users/krisstern/glue/glue/plugins/wcs_autolinking/wcs_autolinking.py", line 32, in get_cids_and_functions
    forwards(*pixel_input)
  File "/Users/krisstern/glue/glue/plugins/wcs_autolinking/wcs_autolinking.py", line 16, in forwards
    return pixel_to_pixel(wcs1, wcs2, *pixel_input)
  File "/opt/miniconda3/envs/astronomy-2020/lib/python3.8/site-packages/astropy/wcs/utils.py", line 798, in pixel_to_pixel
    return wcs_out.world_to_pixel(*world_outputs)
  File "/opt/miniconda3/envs/astronomy-2020/lib/python3.8/site-packages/astropy/wcs/wcsapi/high_level_api.py", line 148, in world_to_pixel
    raise ValueError("Number of world inputs ({}) does not match "
ValueError: Number of world inputs (2) does not match expected (3)

[kakirastern]

Code patched...

But I am xfail-ing one of the new tests for now which should be passing. From experience it could be due to the way the wcs objects used in this particular test case is initiated. The test concerned is test_2d_and_1d_data_cubes_with_no_celestial_axes.

[kakirastern]

Also a caveat, the code currently does not work with 4.1rc1 but with 4.2dev...

[kakirastern]

@astrofrog PR is ready for another review. I have tried simplifying the code structure for the wcs_autolinking.py module to tidy things up.

[astrofrog]

This is good to go, aside from a couple of cleanup items below.

[astrofrog]

Looks great, thank you for addressing all my comments! This is a great contribution 🏆