Add fast method to lookup/estimate WCS limits in world coordinates

[dhomeier]

Adding a method ~ndcube.NDCube.axis_world_coords_limits to return the extension (minimum/maximum in all axes) in world coordinates as laid out in #413 (comment). This returns the appropriate limits as SkyCoord, SpectralCoord objects in the same format as axis_world_coords, to be built into the __str__ representation:

>>> cube.axis_world_coords_limits(max_size=1000)
(<SkyCoord (ICRS): (ra, dec) in deg
     [( 89.44873826, 22.18535743), (169.21775477, 84.43524424)]>,
 <SpectralCoord 
    (target: <ICRS Coordinate: (ra, dec, distance) in (deg, deg, kpc)
                 (90.05, 30.1, 1000.)
              (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                 (0., 0., 0.)>)
   [1.02e-09, 3.00e-09] m>)

The function is largely a clone of axis_world_coords, but adding a max_size parameter above which the call to the full _generate_world_coords evaluation is replaced by estimating the limits from only a subset of pixel coordinates.
My timings indicate that below WCS sizes of some 10000 the former has negligible impact, from some 100000 it approaches 0.1 seconds, so that's probably the limit that __str__ or __repr__ should set when calling it.

TODO etc.:

• I've included basic tests for different dimensionality (just noting that 2d is still missing); not sure if/which of the split/spliced test cases need to be covered in addition. Something for extra_coords and combined_wcs should probably be added, although that is not the core purpose of this function.
• I'd like to return a status flag indicating if the exact/brute force method or an approximation has been used; wondering if this could be stored in a new NDCubeBase property. While at that, it's probably worth discussing if the entire (axes_coords, ) for the limits should simply be cached in a property as well (only to be updated whenever axis_world_coords_limits with higher precision/max_size), so __str__ would only need to read that back.
• As the test for the WCS extending over the pole highlights, simply evaluating min() and max() will just return the pixel closest to the pole for SkyCoord.dec.max() (which lies somewhere along a crpix[i] axis in this example), and whatever pixels fall closest to the meridian on both sides for SkyCoord.ra minmax. Both are of course not representative of the WCS extension on the sky in any way, but I don't have a good idea how to find a better representation in this case. Maybe falling back to the pixel corners would be an option; shifting RA to a range [-180 deg, 180 deg] would also be helpful.

[pep8speaks]

Hello @dhomeier! Thanks for updating this PR.

• In the file ndcube/conftest.py:

Line 544:1: E303 too many blank lines (3)

Comment last updated at 2021-10-27 16:15:21 UTC

[dhomeier]

Doh, taking the extrema on pixel_to_world_values() does not work as expected!
Or at least, in the ndcube_2d_ln_lt case lon wraps around 360˚ in values, but is mapped to [-16", ..., 20"] in high level objects, having the wrong min() and max() returned.
Overcoming this would require identifying the extrema in all Tx, Ty, ra or whatever dimensions separately and then reconstruct the SkyCoord from those, very nasty. But even that would not really resolve cases where the high level coordinate like ra itself wraps over the meridian as in the test_axis_world_coords_rapol case. Basically, how does one define a truly general method to determine a "maximum" in a spherical, periodic coordinate?

[dhomeier]

The different output in the py38-devdeps test probably highlights further issues with the extra_coords treatment.
I don't really understand how the existing axis_world_coords handles those as wcs input or what it is supposed to return. E.g. adding the simple_extra_coords_3d from the tests to a NAXIS (4, 3, 5) ndcube returns

>>> from ndcube.tests import test_ndcollection
>>> cube = test_ndcollection.cube1
>>> cube._extra_coords = ndcube.ExtraCoords.from_lookup_tables(('time', 'hello', 'bye'), (0, 1, 2), (list(range(2)) * u.pix, list(range(3)) * u.pix, list(range(4)) * u.pix))
>>> cube.axis_world_coords()[1]
<SpectralCoord [1.02e-09, 1.04e-09, 1.06e-09, 1.08e-09, 1.10e-09] m>
>>> cube.axis_world_coords(wcs=cube.extra_coords)
(<Quantity [0., 1., 2., 3.] pix>,)
>>> cube.axis_world_coords(wcs=cube.extra_coords.wcs)
(<Quantity [ 0.,  1.,  2.,  3., nan] pix>,)

whereas the same with different units produces

>>> cube._extra_coords = ndcube.ExtraCoords.from_lookup_tables(('time', 'hello', 'bye'), (0, 1, 2), (list(range(2)) * u.second, list(range(3)) * u.pix, list(range(4)) * u.micron))
>>> cube.axis_world_coords(wcs=cube.extra_coords)
(<Quantity [ 0.,  1., nan, nan, nan] s>, <Quantity [0., 1., 2.] pix>, <Quantity [0., 1., 2., 3.] micron>)
>>> cube.axis_world_coords(wcs=cube.extra_coords.wcs)
(<Quantity [ 0.,  1., nan, nan] s>, <Quantity [0., 1., 2.] pix>, <Quantity [ 0.,  1.,  2.,  3., nan] micron>)

Maybe the option to use a different wcs is better just removed from axis_world_coords_limits until this is sorted out.

[dhomeier]

The different output in the py38-devdeps test probably highlights further issues with the extra_coords treatment.
I don't really understand how the existing axis_world_coords handles those as wcs input or what it is supposed to return.

Still holds even for consistent output from that function. Testing the length of the returned output on devdeps has

        coords = ndcube_3d_l_ra_dec.axis_world_coords(wcs=ndcube_3d_l_ra_dec.extra_coords)
>       assert len(coords) == 1
E       assert 3 == 1

which passes on all other platforms and with my local astropy 5.0dev1016 installation with WCSLIB 7.7; no idea what might cause the different behaviour above.

[dhomeier]

which passes on all other platforms and with my local astropy 5.0dev1016 installation with WCSLIB 7.7; no idea what might cause the different behaviour above.

As @Cadair pointed out, the difference is the dev version of gwcs being installed on devdeps; the different behaviour can be isolated to the fact that for gwcs < 0.17using the standard extra_coords_3d produces the (probably) expect results

from ndcube.conftest import data_nd
data = data_nd((2, 3, 4))
cube = ndcube.NDCube(data, wcs_3d_l_lt_ln, mask=data>0, uncertainty=data)
coord0 = time.Time(['2000-01-01T12:00:00', '2000-01-02T12:00:00'], format='isot')
cube._extra_coords = ndcube.ExtraCoords.from_lookup_tables(('time', 'hello', 'bye'), (0, 1, 2), (coord0, list(range(3))*u.pix, list(range(4))*u.m))
cube._extra_coords._ndcube = cube
cube.axis_world_coords(wcs=cube.extra_coords)

(<Time object: scale='utc' format='isot' value=['2000-01-01T12:00:00.000' '2000-01-02T12:00:00.000']>, <Quantity [0., 1., 2.] pix>, <Quantity [0., 1., 2., 3.] m>)

whereas simple_extra_coords_3d produces

cube._extra_coords = ndcube.ExtraCoords.from_lookup_tables(('time', 'hello', 'bye'), (0, 1, 2), (list(range(2))*u.pix, list(range(3))*u.pix, list(range(4))*u.pix))
cube._extra_coords._ndcube = cube
cube.axis_world_coords(wcs=cube.extra_coords)

(<Quantity [0., 1., 2., 3.] pix>,)

so apparently for repeated occurrences of identical physical types gwcs 0.16.1 only keeps the last instance, while the dev version more consistently returns

(<Quantity [ 0.,  1.] pix>, <Quantity [0., 1., 2.] pix>, <Quantity [0., 1., 2., 3.] pix>)

I have not traced back where, how or why that happens, but arguably the dev version is working more correctly here, so we should perhaps refrain from using the extra_coords wcs in these functions more widely until gwcs >= 0.17 can be made a requirement.

[dhomeier]

I am now testing for the new functionality and have it xfail with gwcs 0.16, so marking this ready for a first review.

[Cadair]

My initial thought is I wonder if it would be better if this starts off as private API? I am not sure of it's utility to the general users?

[dhomeier]

Yes, at the start I was relatively indifferent whether it should be public or private, just thought some general uses might spring from it as it develops. The speed differences probably won't matter there, but the output can be much more compact than axis_world_coords.

[Cadair]

There are potentially many situations where we would be passed a HighLevelWCSWrapper, not just those which you list in the comment (mainly a wrapped SlicedLowLevelWCS). A more nuanced check would be to extract the low level object first and then dispatch on the type of that instead.

[dhomeier]

Not sure how to get an overview of what could potentially be passed there.

[Cadair]

I will have to get @DanRyanIrish to check me on this, but I think there is a subset of cases where this _bounding_box_to_points method generates invalid inputs to pixel_to_world.

[dhomeier]

I think that's from the suggestion in #413 (comment), but the link there must have pointed to an older revision.

[Cadair]

I am not sure I am following exactly what this block is doing.

[dhomeier]

The one below?
It takes the world coordinate at crpix[i] in the respective dimension, and creates a grid over the other dependent axes (will probably be usually only one) through that point.

[dhomeier]

For example in

ndcube/ndcube/tests/test_ndcube.py

Lines 468 to 476 in 1efb0b8

	coords = ndcube_3d_l_ra_pol.axis_world_coords_limits(max_size=10000)
	assert u.allclose(coords[0], [1.02e-09, 1.20e-09] * u.m)
	assert u.allclose(coords[1].ra, [26.1484, 333.4424] * u.deg)
	assert u.allclose(coords[1].dec, [61.8572, 89.98945] * u.deg)
	coords = ndcube_3d_l_ra_pol.axis_world_coords_limits(max_size=10)
	assert u.allclose(coords[0], [1.02e-09, 1.20e-09] * u.m)
	assert u.allclose(coords[1].ra, [26.1484, 333.4424] * u.deg)
	assert u.allclose(coords[1].dec, [61.8572, 80.429] * u.deg)

which is working on wcs_3d_l_ra_pol with CRPIX2=200 at CRVAL2=90.5 and CRPIX3 at CRVAL3=80.4, it is spanning a line through all declinations at RA=90.5 deg and all Right Ascensions at Dec=80.4 deg each.
This is covering the point closest to the pole, or very nearly, while the last call, which is skipping this block as it requires more than max_size=10 points, only returns the world coordinates at the bbox corners.
It's of course arguable which is really more useful as limit values...

[Cadair]

Is this identical to the logic at the end of axis_world_coords? if so should we move it into a shared method?

[dhomeier]

Much of the code is; I was just hesitant to overload axis_world_coords (like implementing the whole thing as axis_world_coords(limits=True)).

[stale]

This pull request has been automatically marked as stale because it has not had any activity for the past five months. It will be closed if no further activity occurs. If the ideas in this pull request are still worth implementing, please make sure you open an issue to keep track of that idea!

[stale]

This pull request has been automatically closed since there was no activity for a month since it was marked as stale. If the ideas in this pull request are still worth implementing, please make sure you open an issue to keep track of that idea!

[dhomeier]

Rebasing just to check the tests are still passing; not sure where we last stood re. design decisions.
Probably needs to be reopened first to update again from my branch.