Skip to content

Fix crop world-object order for WCSes with non-contiguous components - #977

Open
nabobalis wants to merge 3 commits into
sunpy:mainfrom
nabobalis:crop-object-order
Open

nabobalis wants to merge 3 commits into
sunpy:mainfrom
nabobalis:crop-object-order

Conversation

@nabobalis

Copy link
Copy Markdown
Member

PR Description

crop takes one entry per world object, with None for objects that should not be cropped. _get_crop_item built that list by popping repeated component names, which keeps the last occurrence of each object; pixel_to_world and world_to_pixel use the first. So on WCSes with interleaved world objects (e.g. HPLN, WAVE, HPLT, or DKIST VISP), crop rejected points in pixel_to_world order, including partial ones, and the count error didn't say what it expected.

import numpy as np

from astropy.wcs import WCS

from ndcube import NDCube

w = WCS(naxis=3)
w.wcs.ctype = ["HPLN-TAN", "WAVE", "HPLT-TAN"]
w.wcs.cunit = ["deg", "m", "deg"]
w.wcs.cdelt = [0.1, 1e-10, 0.1]
w.wcs.crpix = [1, 1, 1]
w.wcs.crval = [0, 5e-7, 0]
w.wcs.dateobs = "2020-01-01"
w.wcs.set()
cube = NDCube(np.zeros((4, 5, 6)), wcs=w)
sky0, wl0 = cube.wcs.pixel_to_world(1, 1, 1)
sky1, wl1 = cube.wcs.pixel_to_world(3, 3, 2)

for label, points in {
    "[sky, wl]": ([sky0, wl0], [sky1, wl1]),
    "[sky, None]": ([sky0, None], [sky1, None]),
    "[None, wl]": ([None, wl0], [None, wl1]),
    "[wl, None]": ([wl0, None], [wl1, None]),
    "[wl, sky]": ([wl0, sky0], [wl1, sky1]),
    "[sky]": ([sky0], [sky1]),
}.items():
    try:
        result = cube.crop(*points).shape
    except Exception as e:
        result = f"{type(e).__name__}: {e}"
    print(f"{label:12} -> {result}")
points main this PR
[sky, wl] (pixel_to_world order) TypeError: <class '...SkyCoord'> of component 0 in point 0 is incompatible with WCS component spectral <class '...Quantity'>. (2, 3, 3)
[sky, None] same TypeError (2, 5, 3)
[None, wl] TypeError: <class '...SpectralCoord'> of component 1 in point 0 is incompatible with WCS component celestial <class '...SkyCoord'>. (4, 3, 6)
[wl, None] ValueError: Expected the following order of world arguments: SkyCoord (4, 3, 6)
[wl, sky] (old order) (2, 3, 3) (2, 3, 3)
[sky] ValueError: 1 components in point 0 do not match WCS with 2 components. same, followed by Each point must have one entry per world object (use None for a component that should not be cropped), in order: celestial (SkyCoord), spectral (Quantity).

Changes:

  • The object order is first appearance (utils.misc.unique_sorted, as in axis_world_coords).
  • If each non-None object is an instance of exactly one world object class (no two the same), objects are matched by class, as in astropy's world_to_pixel; otherwise by position.
  • The count and type errors list the expected objects in order (replaces Name expected world objects in crop component-count error #939).
  • get_crop_item_from_points sorts the pixel axes instead of using set order, which swapped some bounds on cubes with nine or more dimensions.

User-visible changes:

  • pixel_to_world order works, and so does Allow greater flexibility in crop bounds order #608's [spec, sky] example.
  • Old-order points still work when every world object has a distinct class (none a subclass of another), as on DKIST VISP; dkist's crop tests pass, two of them only because of the class matching.
  • Matching is all or nothing per point. With shared classes (FITS WAVE and LINEAR are both Quantity) or subclasses (SpectralCoord/Quantity), the point is positional in first-appearance order, so old-order [wl, sky, None] on HPLN, WAVE, HPLT, LINEAR now raises.
  • Some inputs main rejected with a TypeError now go to the wrong slot: a same-class object in the wrong slot (as with world_to_pixel), a plain Quantity for a gWCS SpectralCoord slot (moved to the only Quantity slot), and, with wcs=cube.extra_coords, a wrongly typed Quantity moved to a dummy pixel slot (I'll open a separate issue for the extra_coords cases).

Related:

  • Related to Allow greater flexibility in crop bounds order #608, not closing it: None placeholders are still required.
  • Conflicts with Typing NDCube #950, which renames points to sanitized_points. points[i] = point = [...] merges cleanly but must become sanitized_points[i] = ..., otherwise it raises 'tuple' object does not support item assignment; test_crop_non_contiguous_world_objects catches this.
  • Please backport to 2.4; it cherry-picks cleanly and the tests pass there.

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

When a crop point has the wrong number of entries, or an entry of the
wrong class, the error now lists the expected world objects in order,
so users migrating from a WCS with fewer world objects, or passing
objects in the wrong order, learn the fix from the error itself.
The component dedup in NDCube._get_crop_item popped entries while
enumerating, which kept the last occurrence of each world object. WCSes
whose object components are not adjacent (e.g. HPLN, WAVE, HPLT, or
DKIST VISP's lon, wl, lat, time, stokes) then rejected points given in
the order pixel_to_world returns, including partial points with None.
Keep the first occurrence instead.

As astropy's world_to_pixel does, also match objects to components by
class when each object in a point matches exactly one class and no two
objects match the same one (sunpy#608). Points in the old order therefore
keep working when the world objects have distinct classes, none a
subclass of another, as in dkist's VISP crop tests. Otherwise, for
example when objects share a class or their classes are related by
subclassing (Quantity and SpectralCoord), the whole point is matched by
position.

Also sort the pixel axes with input in get_crop_item_from_points:
iterating a set swapped the bounds of pixel axes (e.g. 3 and 8) on
cubes with nine or more dimensions.

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