Skip to content

Fix NDCube.rebin for lookup-table extra coords - #982

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

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

Conversation

@nabobalis

Copy link
Copy Markdown
Member

PR Description

Split out of #967. None of this depends on N-D lookup tables; it all reproduces on main and 2.4 with 1-D tables.

NDCube.rebin had three problems with lookup-table extra coords:

  • They were resampled with offset=0, so each new pixel got the value at the first pixel of its bin, while the rebinned WCS gives the bin centre. rebin now passes (bin_shape - 1) / 2 for lookup tables. WCS-backed extra coords keep 0, because ResampledLowLevelWCS treats the offset as a pixel-edge shift. Nothing public tells the two apart, so this checks the private _lookup_tables.
  • ExtraCoords.resample put the new grids in an object array. When they all have the same length (every 1-D cube, for example) numpy makes a 2-D object array and np.interp raises TypeError. They are now a list.
  • TimeTableCoordinate.interpolate interpolated absolute MJD floats (about 0.6 µs resolution for current dates) and dropped reference_time. It now interpolates seconds from the first time and keeps reference_time.

test_rebin already expected bin-centre times, but np.allclose's default rtol=1e-5 is about half a day at MJD 51544, so the 30 s error passed. It now uses rtol=0, atol=1e-8.

The second commit drops the "New array grids must all be same shape" check from QuantityTableCoordinate.interpolate, so a mesh Quantity table works when its axes are rebinned to different lengths. It is separate in case you'd rather not take it.

For direct callers of TimeTableCoordinate.interpolate, or ExtraCoords.resample with kwargs, a numeric left/right is now seconds from table[0] rather than an absolute MJD, and left=right=np.nan gives a NaN Time instead of raising ValueError. rebin passes no kwargs.

import numpy as np
import astropy.units as u
from astropy.time import Time
from astropy.wcs import WCS

from ndcube import NDCube

wcs = WCS(naxis=2)
wcs.wcs.ctype = ["WAVE", "FREQ"]
wcs.wcs.cunit = ["m", "Hz"]

# Time table on array axis 1 at 60 s cadence, rebinned by 2 along that axis.
cube = NDCube(np.zeros((4, 6)), wcs=wcs)
cube.extra_coords.add("time", 1, Time("2026-01-01T00:00:00") + np.arange(6) * 60 * u.s)
r = cube.rebin((1, 2))
print(r.axis_world_coords(wcs=r.extra_coords)[0].isot)

# 1-D cube.
wave = WCS(naxis=1)
wave.wcs.ctype = ["WAVE"]
wave.wcs.cunit = ["m"]
cube1 = NDCube(np.zeros(6), wcs=wave)
cube1.extra_coords.add("wave2", 0, np.arange(6) * u.m)
r1 = cube1.rebin((2,))
print(r1.axis_world_coords(wcs=r1.extra_coords)[0])

main:

['2026-01-01T00:00:00.000' '2026-01-01T00:02:00.000'
 '2026-01-01T00:04:00.000']
Traceback (most recent call last):
  ...
TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'

this branch:

['2026-01-01T00:00:30.000' '2026-01-01T00:02:30.000'
 '2026-01-01T00:04:30.000']
[0.5 2.5 4.5] m

Not in this PR:

This adds one small conflict with #950 in ExtraCoords.resample; resolving it means keeping both changes.

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.

NDCube.rebin resampled lookup-table extra coords with offset 0, so they
were sampled at the first pixel of each bin (0, 2, 4) while the rebinned
WCS describes the bin centres (0.5, 2.5, 4.5). Pass a bin-centre offset
for lookup tables; WCS-backed extra coords keep the pixel-edge offset
that ResampledLowLevelWCS expects. test_rebin already expected the
bin-centre times, but np.allclose's default rtol on MJD values (about
half a day) hid the error, so tighten it.

ExtraCoords.resample also packed the new grids into an object array,
which numpy turns into a 2-D object array when all grids have the same
length, so rebinning a 1-D cube, or any cube whose axes all have the
same length after rebinning, raised "Cannot cast array data from
dtype('O')". Keep the grids in a list.

TimeTableCoordinate.interpolate interpolated absolute MJD floats,
quantising times to ~0.6 us, and dropped reference_time. Interpolate
second offsets from the first time and keep reference_time.

Split out of sunpy#967; none of this depends on N-D lookup tables.
QuantityTableCoordinate always holds a mesh of 1-D tables and
interpolates each table along its own grid, but interpolate still
required every grid to have the same shape. So NDCube.rebin raised
"New array grids must all be same shape" for a Quantity lookup table
spanning several axes whose lengths differ after rebinning. Drop the
check.

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