diff --git a/changelog/982.bugfix.1.rst b/changelog/982.bugfix.1.rst new file mode 100644 index 000000000..427c819be --- /dev/null +++ b/changelog/982.bugfix.1.rst @@ -0,0 +1 @@ +Interpolating a `~ndcube.extra_coords.TimeTableCoordinate`, e.g. when rebinning, now keeps sub-microsecond precision and the coordinate's ``reference_time``. diff --git a/changelog/982.bugfix.2.rst b/changelog/982.bugfix.2.rst new file mode 100644 index 000000000..7ae4ecb89 --- /dev/null +++ b/changelog/982.bugfix.2.rst @@ -0,0 +1 @@ +~ndcube.extra_coords.QuantityTableCoordinate.interpolate now accepts grids of different lengths, so ~ndcube.NDCube.rebin no longer raises "New array grids must all be same shape" when a ~astropy.units.Quantity lookup table spans several axes that are rebinned to different lengths. diff --git a/changelog/982.bugfix.rst b/changelog/982.bugfix.rst new file mode 100644 index 000000000..db494b058 --- /dev/null +++ b/changelog/982.bugfix.rst @@ -0,0 +1 @@ +~ndcube.NDCube.rebin now samples lookup-table extra coords at the centre of each new pixel, matching the rebinned WCS, instead of at the first pixel of each bin. It also no longer raises a TypeError when every axis of the rebinned cube has the same length, which is always the case for 1-D cubes. diff --git a/ndcube/extra_coords/extra_coords.py b/ndcube/extra_coords/extra_coords.py index bcbf6724d..020d2ac50 100644 --- a/ndcube/extra_coords/extra_coords.py +++ b/ndcube/extra_coords/extra_coords.py @@ -493,12 +493,11 @@ def resample(self, factor, offset=0, ndcube=None, **kwargs): x = np.arange(c, d+f, f) x = x[x <= d-1] new_grids.append(x) - new_grids = np.array(new_grids, dtype=object) for array_axes, coord in self._lookup_tables: if np.isscalar(array_axes): new_coord = coord.interpolate(new_grids[array_axes], **kwargs) else: - new_coord = coord.interpolate(*new_grids[np.asarray(array_axes)], **kwargs) + new_coord = coord.interpolate(*[new_grids[axis] for axis in array_axes], **kwargs) new_ec.add(coord.names, array_axes, new_coord, physical_types=coord.physical_types) return new_ec diff --git a/ndcube/extra_coords/table_coord.py b/ndcube/extra_coords/table_coord.py index 0e18b2d3b..b72f877f9 100644 --- a/ndcube/extra_coords/table_coord.py +++ b/ndcube/extra_coords/table_coord.py @@ -448,9 +448,8 @@ def interpolate(self, *new_array_grids, **kwargs): new_array_grids: array-like The array index values at which the the new values of the coords are desired. An array grid must be provided as a separate arg - for each array dimension and corresponding elements in all arrays - represent a single location in the pixel grid. Therefore, array grids - must all have the same shape. + for each array dimension. Each table is interpolated independently + along its own grid, so the grids can have different lengths. kwargs All remaining kwargs are passed to underlying interpolation function. @@ -468,8 +467,6 @@ def interpolate(self, *new_array_grids, **kwargs): if len(new_array_grids) != ndim: raise ValueError( f"A new array grid must be given for each array axis/table, i.e. {ndim}") - if any(new_grid.shape != new_array_grids[0].shape for new_grid in new_array_grids): - raise ValueError("New array grids must all be same shape.") # Build array grids for non-interpolated table. old_array_grids = tuple(np.arange(d) for d in self.shape) # Iterate through tables and interpolate each. @@ -796,12 +793,13 @@ def interpolate(self, new_array_grids, **kwargs): raise ValueError("Cannot interpolate a scalar TimeTableCoordinate.") # Build pixel grids for current TimeTableCoord. old_array_grids = np.arange(len(self.table)) - # Interpolate using MJD format and convert back to a Time object. - new_table = np.interp(new_array_grids, old_array_grids, self.table.mjd, **kwargs) - new_table = Time(new_table, scale=self.table.scale, format="mjd") - new_table.format = self.table.format + # Interpolate offsets from the first time; absolute MJD floats lose sub-microsecond precision. + origin = self.table[0] + deltas = (self.table - origin).to_value(u.s) + new_table = origin + np.interp(new_array_grids, old_array_grids, deltas, **kwargs) * u.s # Rebuild new TimeTableCoord and return. - new_coord = type(self)(new_table, names=self.names, physical_types=self.physical_types) + new_coord = type(self)(new_table, names=self.names, physical_types=self.physical_types, + reference_time=self.reference_time) new_coord._dropped_world_dimensions = self._dropped_world_dimensions return new_coord diff --git a/ndcube/extra_coords/tests/test_lookup_table_coord.py b/ndcube/extra_coords/tests/test_lookup_table_coord.py index b4e86878f..fda24db6d 100644 --- a/ndcube/extra_coords/tests/test_lookup_table_coord.py +++ b/ndcube/extra_coords/tests/test_lookup_table_coord.py @@ -658,6 +658,15 @@ def test_time_interpolate(lut_1d_time): assert_lutc_ancilliary_data_same(output, lutc) +def test_time_interpolate_keeps_precision_and_reference_time(): + reference = Time("2020-01-01T00:00:00", scale="utc") + times = Time("2026-01-01T00:00:00", scale="utc") + np.arange(4) * 10 * u.ns + coord = TimeTableCoordinate(times, reference_time=reference) + new = coord.interpolate(np.arange(4.)) + assert new.reference_time == reference + assert u.allclose((new.table - times).to(u.ns), 0 * u.ns, atol=1 * u.ns) + + def test_skycoord_interpolate_no_mesh(lut_2d_skycoord_no_mesh): lutc = lut_2d_skycoord_no_mesh new_array_grids = np.meshgrid(np.arange(0.5, 2), np.arange(0, 3)) @@ -696,10 +705,6 @@ def test_quantity_interpolate_errors(): qtc.interpolate(np.ones(1)) assert "A new array grid must be given for each array axis" in str(ei) - with pytest.raises(ValueError) as ei: - qtc.interpolate(np.ones(1), np.ones(2)) - assert "New array grids must all be same shape." in str(ei) - def test_skycoord_interpolate_error(lut_2d_skycoord_mesh): sctc = SkyCoordTableCoordinate(SkyCoord(1 * u.deg, 2 * u.deg)) diff --git a/ndcube/ndcube.py b/ndcube/ndcube.py index 87ece2c04..f86b6dc1c 100644 --- a/ndcube/ndcube.py +++ b/ndcube/ndcube.py @@ -1392,7 +1392,9 @@ def my_propagate(uncertainty, data, mask, **kwargs): new_cube._global_coords = self._global_coords # Reconstitute extra coords if not self.extra_coords.is_empty: - new_cube._extra_coords = self.extra_coords.resample(bin_shape, ndcube=new_cube) + # Lookup tables are sampled at bin centres; WCS wrappers take pixel-edge offsets. + offset = (bin_shape - 1) / 2 if self.extra_coords._lookup_tables else 0 + new_cube._extra_coords = self.extra_coords.resample(bin_shape, offset=offset, ndcube=new_cube) return new_cube diff --git a/ndcube/tests/test_ndcube_reproject_and_rebin.py b/ndcube/tests/test_ndcube_reproject_and_rebin.py index b9bef439a..71f89e9ee 100644 --- a/ndcube/tests/test_ndcube_reproject_and_rebin.py +++ b/ndcube/tests/test_ndcube_reproject_and_rebin.py @@ -157,7 +157,7 @@ def test_rebin(ndcube_3d_l_ln_lt_ectime, bin_shape): assert u.allclose(output_spec, expected_spec) assert output_time.scale == expected_time.scale assert output_time.format == expected_time.format - assert np.allclose(output_time.mjd, expected_time.mjd) + assert np.allclose(output_time.mjd, expected_time.mjd, rtol=0, atol=1e-8) def test_rebin_dask(ndcube_2d_dask): @@ -325,3 +325,18 @@ def test_rebin_specutils(): spec = Spectrum(flux=y, spectral_axis=x, bin_specification='centers', mask=x > 2000*u.nm) output = spec.rebin((10,), operation=np.sum, operation_ignores_mask=False) assert output.shape == (400,) + + +def test_rebin_1d_extra_coords(ndcube_1d_l): + ndcube_1d_l.extra_coords.add("wave", 0, ndcube_1d_l.axis_world_coords_values()[0]) + output = ndcube_1d_l.rebin((2,)) + assert u.allclose(output.axis_world_coords_values(wcs=output.extra_coords)[0], + output.axis_world_coords_values()[0]) + + +def test_rebin_mesh_extra_coords(ndcube_2d_ln_lt): + ndcube_2d_ln_lt.extra_coords.add(("a", "b"), (0, 1), (np.arange(10) * u.m, np.arange(12) * u.m)) + output = ndcube_2d_ln_lt.rebin((5, 4)) + a, b = output.axis_world_coords(wcs=output.extra_coords) + assert u.allclose(a, [2, 7] * u.m) + assert u.allclose(b, [1.5, 5.5, 9.5] * u.m)