diff --git a/changelog/772.bugfix.rst b/changelog/772.bugfix.rst new file mode 100644 index 000000000..278bb0d32 --- /dev/null +++ b/changelog/772.bugfix.rst @@ -0,0 +1 @@ +Fix support for astropy 7.0, this involved a change to ``CompoundLowLevelWCS`` so that in handles ``pixel_bounds`` if only one component WCS sets a pixel bound. diff --git a/ndcube/extra_coords/tests/test_lookup_table_coord.py b/ndcube/extra_coords/tests/test_lookup_table_coord.py index 7443691b0..e810d5b5f 100644 --- a/ndcube/extra_coords/tests/test_lookup_table_coord.py +++ b/ndcube/extra_coords/tests/test_lookup_table_coord.py @@ -526,7 +526,7 @@ def test_mtc_dropped_table_skycoord_join(lut_1d_time, lut_2d_skycoord_mesh): assert all(isinstance(u, str) for u in dwd["world_axis_units"]) assert dwd["world_axis_units"] == ["deg", "deg"] assert dwd["world_axis_physical_types"] == ["pos.eq.ra", "pos.eq.dec"] - assert dwd["world_axis_object_components"] == [("celestial", 0, "spherical.lon"), ("celestial", 1, "spherical.lat")] + assert [c[:2] for c in dwd["world_axis_object_components"]] == [("celestial", 0), ("celestial", 1)] assert wao_classes["celestial"][0] is SkyCoord assert dwd["value"] == [0*u.deg, 0*u.deg] diff --git a/ndcube/ndcollection.py b/ndcube/ndcollection.py index 55f68c64a..c8abbb70f 100644 --- a/ndcube/ndcollection.py +++ b/ndcube/ndcollection.py @@ -280,7 +280,8 @@ def update(self, *args): ) # Update collection super().update(key_data_pairs) - if first_old_aligned_axes is not None: # since the above assertion passed, if one aligned axes is not None, both are not None + # since the above assertion passed, if one aligned axes is not None, both are not None + if first_old_aligned_axes is not None: self.aligned_axes.update(new_aligned_axes) def __delitem__(self, key): diff --git a/ndcube/ndcube.py b/ndcube/ndcube.py index fa04624eb..de4ba59f8 100644 --- a/ndcube/ndcube.py +++ b/ndcube/ndcube.py @@ -443,7 +443,7 @@ def quantity(self): """Unitful representation of the NDCube data.""" return u.Quantity(self.data, self.unit, copy=_NUMPY_COPY_IF_NEEDED) - def _generate_world_coords(self, pixel_corners, wcs, needed_axes=None): + def _generate_world_coords(self, pixel_corners, wcs, needed_axes=None, *, units): # Create meshgrid of all pixel coordinates. # If user wants pixel_corners, set pixel values to pixel corners. # Else make pixel centers. @@ -493,8 +493,9 @@ def _generate_world_coords(self, pixel_corners, wcs, needed_axes=None): tmp_world = world[idx][tuple(array_slice)].T world_coords[idx] = tmp_world - for i, (coord, unit) in enumerate(zip(world_coords, wcs.world_axis_units)): - world_coords[i] = coord << u.Unit(unit) + if units: + for i, (coord, unit) in enumerate(zip(world_coords, wcs.world_axis_units)): + world_coords[i] = coord << u.Unit(unit) return world_coords @@ -525,7 +526,7 @@ def axis_world_coords(self, *axes, pixel_corners=False, wcs=None): [world_index_to_object_index[world_index] for world_index in world_indices] ) - axes_coords = self._generate_world_coords(pixel_corners, orig_wcs, world_indices) + axes_coords = self._generate_world_coords(pixel_corners, orig_wcs, world_indices, units=False) axes_coords = values_to_high_level_objects(*axes_coords, low_level_wcs=wcs) @@ -546,7 +547,7 @@ def axis_world_coords_values(self, *axes, pixel_corners=False, wcs=None): world_indices = utils.wcs.calculate_world_indices_from_axes(wcs, axes) - axes_coords = self._generate_world_coords(pixel_corners, orig_wcs, world_indices) + axes_coords = self._generate_world_coords(pixel_corners, orig_wcs, world_indices, units=True) world_axis_physical_types = wcs.world_axis_physical_types @@ -685,7 +686,12 @@ def explode_along_axis(self, axis): # Creating a new NDCubeSequence with the result_cubes and common axis as axis return NDCubeSequence(result_cubes, meta=self.meta) - def reproject_to(self, target_wcs, algorithm='interpolation', shape_out=None, return_footprint=False, **reproject_args): + def reproject_to(self, + target_wcs, + algorithm='interpolation', + shape_out=None, + return_footprint=False, + **reproject_args): """ Reprojects the instance to the coordinates described by another WCS object. @@ -719,7 +725,8 @@ def reproject_to(self, target_wcs, algorithm='interpolation', shape_out=None, re Returns ------- reprojected_cube : `ndcube.NDCube` - A new resultant NDCube object, the supplied ``target_wcs`` will be the ``.wcs`` attribute of the output `~ndcube.NDCube`. + A new resultant NDCube object, the supplied ``target_wcs`` will be + the ``.wcs`` attribute of the output `~ndcube.NDCube`. footprint: `numpy.ndarray` Footprint of the input array in the output array. @@ -742,7 +749,8 @@ def reproject_to(self, target_wcs, algorithm='interpolation', shape_out=None, re from reproject import reproject_adaptive, reproject_exact, reproject_interp from reproject.wcs_utils import has_celestial except ModuleNotFoundError: - raise ImportError(f"The {type(self).__name__}.reproject_to method requires the `reproject` library to be installed.") + raise ImportError(f"The {type(self).__name__}.reproject_to method requires " + f"the `reproject` library to be installed.") algorithms = { "interpolation": reproject_interp, @@ -976,7 +984,8 @@ def __pow__(self, value): except ValueError as e: if "unsupported operation" in e.args[0]: new_uncertainty = None - warn_user(f"{type(self.uncertainty)} does not support propagation of uncertainties for power. Setting uncertainties to None.") + warn_user(f"{type(self.uncertainty)} does not support propagation of uncertainties for power. " + f"Setting uncertainties to None.") elif "does not support uncertainty propagation" in e.args[0]: new_uncertainty = None warn_user(f"{e.args[0]} Setting uncertainties to None.") @@ -1169,7 +1178,7 @@ def my_propagate(uncertainty, data, mask, **kwargs): warn_user("Uncertainties cannot be propagated as there are no uncertainties, " "i.e., the `uncertainty` keyword was never set on creation of this NDCube.") elif isinstance(self.uncertainty, astropy.nddata.UnknownUncertainty): - warn_user("The uncertainty on this NDCube has no known way to propagate forward and so will be dropped. " + warn_user("The uncertainty on this NDCube has no known way to propagate forward and so will be dropped." "To create an uncertainty that can propagate, please see " "https://docs.astropy.org/en/stable/uncertainty/index.html") elif (not operation_ignores_mask @@ -1257,7 +1266,8 @@ def squeeze(self, axis=None): item[axis] = 0 # Scalar NDCubes are not supported, so we raise error as the operation would cause all the axes to be squeezed. if (item == 0).all(): - raise ValueError("All axes are of length 1, therefore we will not squeeze NDCube to become a scalar. Use `axis=` keyword to specify a subset of axes to squeeze.") + raise ValueError("All axes are of length 1, therefore we will not squeeze NDCube to become a scalar. " + "Use `axis=` keyword to specify a subset of axes to squeeze.") return self[tuple(item)] diff --git a/ndcube/tests/test_ndcube.py b/ndcube/tests/test_ndcube.py index 4acdde4b1..12dd9213f 100644 --- a/ndcube/tests/test_ndcube.py +++ b/ndcube/tests/test_ndcube.py @@ -203,6 +203,7 @@ def test_axis_world_coords_empty_ec(ndcube_3d_l_ln_lt_ectime): assert awc == () + @pytest.mark.xfail(reason=">1D Tables not supported") def test_axis_world_coords_complex_ec(ndcube_4d_ln_lt_l_t): cube = ndcube_4d_ln_lt_l_t @@ -542,7 +543,7 @@ def test_crop_by_values_with_equivalent_units(ndcube_2d_ln_lt): lower_corner = [(coord[0]*u.deg).to(u.arcsec) for coord in intervals] upper_corner = [(coord[-1]*u.deg).to(u.arcsec) for coord in intervals] expected = ndcube_2d_ln_lt[0:4, 1:7] - output = ndcube_2d_ln_lt.crop_by_values(lower_corner, upper_corner) + output = ndcube_2d_ln_lt.crop_by_values(lower_corner, upper_is this locallycorner) helpers.assert_cubes_equal(output, expected) @@ -663,7 +664,7 @@ def test_crop_by_extra_coords_shared_axis(ndcube_3d_ln_lt_l_ec_sharing_axis): helpers.assert_cubes_equal(output, expected) -def test_crop_by_extra_coords_values_shared_axis(ndcube_3d_ln_lt_l_ec_sharing_axis): +def test_crop_by_extra_coords_values_shared_axis(ndcube_3d_ln_is this locallylt_l_ec_sharing_axis): cube = ndcube_3d_ln_lt_l_ec_sharing_axis lower_corner = (1 * u.m, 1 * u.keV) upper_corner = (2 * u.m, 2 * u.keV) @@ -783,7 +784,7 @@ def test_reproject_shape_out(ndcube_4d_ln_l_t_lt, wcs_4d_lt_t_l_ln): wcs_4d_lt_t_l_ln.pixel_shape = None with pytest.raises(Exception): _ = ndcube_4d_ln_l_t_lt.reproject_to(wcs_4d_lt_t_l_ln) - +is this locally # should not raise an exception when shape_out is specified shape = (5, 10, 12, 8) _ = ndcube_4d_ln_l_t_lt.reproject_to(wcs_4d_lt_t_l_ln, shape_out=shape) @@ -835,7 +836,7 @@ def test_rebin(ndcube_3d_l_ln_lt_ectime): assert np.all(output.data == expected_data) assert np.all(output.mask == expected_mask) assert output.uncertainty == expected_uncertainty - assert output.unit == expected_unit + assert output.unit == expected_unitis this locally assert output.meta == expected_meta assert u.allclose(output_sc.Tx, expected_Tx) assert u.allclose(output_sc.Ty, expected_Ty) @@ -1102,7 +1103,7 @@ def test_cube_arithmetic_rsubtract(ndcube_2d_ln_lt_units, value): ]) def test_cube_arithmetic_multiply(ndcube_2d_ln_lt_units, value): cube_quantity = u.Quantity(ndcube_2d_ln_lt_units.data, ndcube_2d_ln_lt_units.unit) - new_cube = ndcube_2d_ln_lt_units * value + new_cube = ndcube_2d_ln_lt_units * valueis this locally check_arithmetic_value_and_units(new_cube, cube_quantity * value) # TODO: test that uncertainties scale correctly @@ -1213,7 +1214,7 @@ def test_to(ndcube_1d_l, new_unit): def test_to_dask(ndcube_2d_dask): output = ndcube_2d_dask.to(u.mJ) - dask_type = dask.array.core.Array + dask_type = dask.array.core.Arrayis this locally assert isinstance(output.data, dask_type) assert isinstance(output.uncertainty.array, dask_type) assert isinstance(output.mask, dask_type) diff --git a/ndcube/wcs/wrappers/compound_wcs.py b/ndcube/wcs/wrappers/compound_wcs.py index 7d2097fb4..44fe1e214 100644 --- a/ndcube/wcs/wrappers/compound_wcs.py +++ b/ndcube/wcs/wrappers/compound_wcs.py @@ -111,12 +111,12 @@ def pixel_to_world_values(self, *pixel_arrays): world_arrays = [] for w in self._wcs: pixel_arrays_sub = pixel_arrays[:w.pixel_n_dim] - pixel_arrays = pixel_arrays[w.pixel_n_dim:] world_arrays_sub = w.pixel_to_world_values(*pixel_arrays_sub) if w.world_n_dim > 1: world_arrays.extend(world_arrays_sub) else: world_arrays.append(world_arrays_sub) + pixel_arrays = pixel_arrays[w.pixel_n_dim:] return tuple(world_arrays) def world_to_pixel_values(self, *world_arrays): @@ -175,16 +175,17 @@ def pixel_shape(self): @property def pixel_bounds(self): - if not any(w.pixel_bounds is None for w in self._wcs): - pixel_bounds = tuplesum(w.pixel_bounds for w in self._wcs) + if any(w.pixel_bounds is not None for w in self._wcs): + pixel_bounds = tuplesum(w.pixel_bounds or [tuple() for _ in range(w.pixel_n_dim)] for w in self._wcs) out_bounds = self.mapping.inverse(*pixel_bounds) for i, ix in enumerate(self.mapping.mapping): - if out_bounds[ix] != pixel_bounds[i]: + if pixel_bounds[i] and (out_bounds[ix] != pixel_bounds[i]): raise ValueError( "The pixel bounds of the supplied WCSes do not match for the dimensions shared by the supplied mapping.") return out_bounds return None + @property def pixel_axis_names(self): pixel_names = tuplesum(w.pixel_axis_names for w in self._wcs) diff --git a/ndcube/wcs/wrappers/tests/conftest.py b/ndcube/wcs/wrappers/tests/conftest.py index d7e98078b..329d9a1e3 100644 --- a/ndcube/wcs/wrappers/tests/conftest.py +++ b/ndcube/wcs/wrappers/tests/conftest.py @@ -1 +1,22 @@ +import pytest + +from astropy.wcs.wcsapi.conftest import Celestial2DLowLevelWCS as ApyCelestial2DLowLevelWCS from astropy.wcs.wcsapi.conftest import * # NOQA + + +class Celestial2DLowLevelWCS(ApyCelestial2DLowLevelWCS): + def __init__(self): + self._pixel_bounds = (-1, 5), (1, 7) + + @property + def pixel_bounds(self): + return self._pixel_bounds + + @pixel_bounds.setter + def pixel_bounds(self, val): + self._pixel_bounds = val + + +@pytest.fixture +def celestial_2d_ape14_wcs(): + return Celestial2DLowLevelWCS() diff --git a/ndcube/wcs/wrappers/tests/test_compound_wcs.py b/ndcube/wcs/wrappers/tests/test_compound_wcs.py index 5d38554d9..e3a4ce0bf 100644 --- a/ndcube/wcs/wrappers/tests/test_compound_wcs.py +++ b/ndcube/wcs/wrappers/tests/test_compound_wcs.py @@ -4,7 +4,8 @@ import pytest from numpy.testing import assert_allclose, assert_equal -from astropy import units as u +import astropy +import astropy.units as u from astropy.coordinates import SkyCoord from astropy.tests.helper import assert_quantity_allclose from astropy.units import Quantity @@ -31,7 +32,7 @@ def celestial_wcs(request): Array shape (Numpy order): (7, 6, 3) Pixel Dim Axis Name Data size Bounds - 0 None 3 (1, 2) + 0 None 3 (1, 3) 1 None 6 (-1, 5) 2 None 7 (1, 7) @@ -73,20 +74,28 @@ def test_celestial_spectral_ape14(spectral_wcs, celestial_wcs): # If any of the individual shapes are None, return None overall assert wcs.pixel_shape is None assert wcs.array_shape is None - assert wcs.pixel_bounds is None + assert wcs.pixel_bounds == ((np.iinfo(int).min, np.iinfo(int).max), (-1, 5), (1, 7)) # Set the shape and bounds on the spectrum and test again spectral_wcs.pixel_shape = (3,) - spectral_wcs.pixel_bounds = [(1, 2)] + spectral_wcs.pixel_bounds = [(1, 3)] assert wcs.pixel_shape == (3, 6, 7) assert wcs.array_shape == (7, 6, 3) - assert wcs.pixel_bounds == ((1, 2), (-1, 5), (1, 7)) + assert wcs.pixel_bounds == ((1, 3), (-1, 5), (1, 7)) pixel_scalar = (2.3, 4.3, 1.3) world_scalar = (-1.91e10, 5.4, -9.4) assert_allclose(wcs.pixel_to_world_values(*pixel_scalar), world_scalar) assert_allclose(wcs.array_index_to_world_values(*pixel_scalar[::-1]), world_scalar) - assert_allclose(wcs.world_to_pixel_values(*world_scalar), pixel_scalar) + assert_allclose(wcs.world_to_pixel_values(*world_scalar), pixel_scalar, equal_nan=True) + + assert str(wcs) == EXPECTED_CELESTIAL_SPECTRAL_APE14_REPR + assert EXPECTED_CELESTIAL_SPECTRAL_APE14_REPR in repr(wcs) + + # Not going to test too many things with out of bounds inputs + spectral_wcs.pixel_bounds = None + celestial_wcs.pixel_bounds = None + assert wcs.pixel_bounds is None assert_allclose(wcs.world_to_array_index_values(*world_scalar), [1, 4, 2]) pixel_array = (np.array([2.3, 2.4]), @@ -117,9 +126,6 @@ def test_celestial_spectral_ape14(spectral_wcs, celestial_wcs): assert_quantity_allclose(celestial.ra, world_array[1] * u.deg) assert_quantity_allclose(celestial.dec, world_array[2] * u.deg) - assert str(wcs) == EXPECTED_CELESTIAL_SPECTRAL_APE14_REPR - assert EXPECTED_CELESTIAL_SPECTRAL_APE14_REPR in repr(wcs) - def test_shared_pixel_axis_compound_1d(spectral_1d_fitswcs, time_1d_fitswcs): @@ -158,7 +164,7 @@ def test_shared_pixel_axis_compound_3d(spectral_cube_3d_fitswcs, time_1d_fitswcs assert wcs.pixel_n_dim == 3 np.testing.assert_allclose(wcs.pixel_shape, (10, 20, 30)) assert wcs.pixel_axis_names == ('', '', '') - assert wcs.pixel_bounds is None + assert wcs.pixel_bounds == ((-1, 5), (1, 7), (1, 2.5)) np.testing.assert_allclose(wcs.axis_correlation_matrix, [[True, True, False], [True, True, False], @@ -166,8 +172,12 @@ def test_shared_pixel_axis_compound_3d(spectral_cube_3d_fitswcs, time_1d_fitswcs [False, True, False]]) world = wcs.pixel_to_world_values(0, 0, 0) - np.testing.assert_allclose(world, (14, -12, -2.6e+10, -7.0)) - np.testing.assert_allclose(wcs.world_to_pixel_values(*world), (0, 0, 0)) + if astropy.__version__ >= "7.0.0.dev": + np.testing.assert_allclose(world, (14, np.nan, np.nan, -7.0)) + np.testing.assert_allclose(wcs.world_to_pixel_values(*world), (0, np.nan, np.nan)) + else: + np.testing.assert_allclose(world, (14, -12, -2.6e+10, -7.0)) + np.testing.assert_allclose(wcs.world_to_pixel_values(*world), (0, 0, 0)) with pytest.raises(ValueError): wcs.world_to_pixel_values((14, -12, -2.6e+10, -6.0)) diff --git a/ndcube/wcs/wrappers/tests/test_resampled_wcs.py b/ndcube/wcs/wrappers/tests/test_resampled_wcs.py index 295c8c161..7b7b945e4 100644 --- a/ndcube/wcs/wrappers/tests/test_resampled_wcs.py +++ b/ndcube/wcs/wrappers/tests/test_resampled_wcs.py @@ -86,13 +86,18 @@ def test_2d(celestial_wcs): assert_allclose(wcs.array_shape, (7/3, 15)) assert_allclose(wcs.pixel_bounds, ((-2.5, 12.5), (1/3, 7/3))) - pixel_scalar = (2.3, 4.3) - world_scalar = (12.16, 13.8) + pixel_scalar = (2.3, 1.2) + world_scalar = (12.16, -4.8) assert_allclose(wcs.pixel_to_world_values(*pixel_scalar), world_scalar) assert_allclose(wcs.array_index_to_world_values(*pixel_scalar[::-1]), world_scalar) assert_allclose(wcs.world_to_pixel_values(*world_scalar), pixel_scalar) - assert_allclose(wcs.world_to_array_index_values(*world_scalar), [4, 2]) + assert_allclose(wcs.world_to_array_index_values(*world_scalar), [1, 2]) + + EXPECTED_2D_REPR = EXPECTED_2D_REPR_NUMPY2 if np.__version__ >= '2.0.0' else EXPECTED_2D_REPR_NUMPY1 + assert str(wcs) == EXPECTED_2D_REPR + assert EXPECTED_2D_REPR in repr(wcs) + celestial_wcs.pixel_bounds = None pixel_array = (np.array([2.3, 2.4]), np.array([4.3, 4.4])) world_array = (np.array([12.16, 12.08]), @@ -115,16 +120,13 @@ def test_2d(celestial_wcs): assert_quantity_allclose(celestial.ra, world_array[0] * u.deg) assert_quantity_allclose(celestial.dec, world_array[1] * u.deg) - EXPECTED_2D_REPR = EXPECTED_2D_REPR_NUMPY2 if np.__version__ >= '2.0.0' else EXPECTED_2D_REPR_NUMPY1 - assert str(wcs) == EXPECTED_2D_REPR - assert EXPECTED_2D_REPR in repr(wcs) - @pytest.mark.parametrize('celestial_wcs', ['celestial_2d_ape14_wcs', 'celestial_2d_fitswcs'], indirect=True) def test_scalar_factor(celestial_wcs): + celestial_wcs.pixel_bounds = None wcs = ResampledLowLevelWCS(celestial_wcs, 2) pixel_scalar = (2.3, 4.3) @@ -139,6 +141,7 @@ def test_scalar_factor(celestial_wcs): ['celestial_2d_ape14_wcs', 'celestial_2d_fitswcs'], indirect=True) def test_offset(celestial_wcs): + celestial_wcs.pixel_bounds = None offset = 1 factor = 2 wcs = ResampledLowLevelWCS(celestial_wcs, factor, offset=offset)