Skip to content

Fix slicing and transposed world coords of multi-axis extra coords - #980

Open
nabobalis wants to merge 3 commits into
sunpy:mainfrom
nabobalis:extra-coords-axes
Open

nabobalis wants to merge 3 commits into
sunpy:mainfrom
nabobalis:extra-coords-axes

Conversation

@nabobalis

Copy link
Copy Markdown
Member

PR Description

Split out of #941 and #967. It makes two independent fixes, one commit each, for extra coords that span more than one array axis.

Slicing (ExtraCoords._getitem_lookup_tables): table axes given as a list were treated as one axis, so any slice raised TypeError. ASDF stores the axes as a list, so this hit every such table loaded from ASDF. Integer indexing also kept the dropped axis ((0, 1) became (-1, 0)), which gave the wrong length or NaN on a cube with more axes than the table.

Transposed output (_generate_world_coords): correlated world coords were transposed with a blind .T, which is only right when the WCS pixel inputs are in ascending cube-axis order. That is true for the cube's WCS, but not for an ExtraCoords WCS, whose inputs follow ExtraCoords.mapping. So axis_world_coords(wcs=cube.extra_coords) transposed a 2-D table on (0, 1), and also WCS-backed extra coords whose mapping reorders correlated axes. The transpose now follows the mapping. For the cube's own WCS it is still the full reversal, so ordinary WCS output is unchanged, and the existing axis_world_coords tests guard that.

import numpy as np, astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
from ndcube import NDCube

lon = np.arange(12.).reshape(3, 4)
for axes in [(0, 1), [0, 1]]:
    cube = NDCube(np.zeros((3, 4, 5)), wcs=WCS(naxis=3))
    cube.extra_coords.add(("lon", "lat"), axes, SkyCoord(lon * u.deg, lon * 0 * u.deg), mesh=False)
    for item in [np.s_[:], np.s_[1]]:
        try:
            sliced = cube[item]
            out = sliced.axis_world_coords(wcs=sliced.extra_coords)[0].ra.deg
            print(axes, item, out.shape, "ok" if np.array_equal(out, lon[item]) else "WRONG")
        except Exception as e:
            print(axes, item, type(e).__name__, e)

main:

(0, 1) slice(None, None, None) (4, 3) WRONG
(0, 1) 1 (5,) WRONG
[0, 1] slice(None, None, None) TypeError list indices must be integers or slices, not list
[0, 1] 1 TypeError list indices must be integers or slices, not list

this branch:

(0, 1) slice(None, None, None) (3, 4) ok
(0, 1) 1 (4,) ok
[0, 1] slice(None, None, None) (3, 4) ok
[0, 1] 1 (4,) ok

Behaviour change: a table is now matched to its array axes in the order given. Code that used #342's workaround on a square table, add(..., (0, 1), table.T), got the intended array because the two transposes cancelled. It now gets the transpose. The fix is to drop the .T or list the axes descending, and the changelog says so. Non-square tables used this way gave NaN before and still do.

Known interaction: rebinning a non-mesh 2-D SkyCoord extra coord leaves a 1-D table declared on two axes. On main, axis_world_coords(wcs=cube.extra_coords) then returns wrong (diagonal) values; with this PR it raises IndexError. A follow-up after the rebin fix will add a meshgrid in ExtraCoords.resample.

#950 edits the same new_lut_axes line, so expect a small conflict.

Related to #342

AI Assistance Disclosure

AI tools were used for:

  • Code generation (e.g., when writing an implementation or fixing a bug)
  • Test/benchmark generation
  • Documentation (including examples)
  • Research and understanding
  • No AI tools were used

Regardless of AI use, the human contributor remains fully responsible for correctness, design choices, licensing compatibility, and long-term maintainability.

K

ExtraCoords._getitem_lookup_tables treated any axes that were not a
tuple as a single axis, and kept every axis of a table after integer
slicing. Tables added with a list of axes, and every multi-axis lookup
table loaded from ASDF (which stores the axes as a list), therefore
raised "list indices must be integers or slices" when the cube was
sliced. Integer indexing also left stale axes on the table (cube[1]
turned (0, 1) into (-1, 0)), so on a cube with more axes than the table
the extra coord was read along the wrong array axis, giving wrong or
NaN values.

Normalise the axes to a tuple and drop integer-sliced axes.
_generate_world_coords transposed each block of correlated world
coordinates with .T, which is only right when the WCS pixel inputs are
in ascending cube pixel order. That holds for the cube's own WCS, but an
ExtraCoords WCS takes its inputs in the order of its mapping, which for
lookup tables follows the order the array axes were given in. So a 2-D
SkyCoord table added on array axes (0, 1) came out of axis_world_coords
transposed, as did WCS-backed extra coords whose mapping reorders
correlated pixel axes.

Order each block by descending cube pixel axis instead. For the cube's
own WCS this is the same reversal as .T, so its output is unchanged.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant