diff --git a/changelog/857.bugfix.rst b/changelog/857.bugfix.rst new file mode 100644 index 000000000..f29c88e84 --- /dev/null +++ b/changelog/857.bugfix.rst @@ -0,0 +1 @@ +Fix translations in `~ndcube.wcs.wrappers.resampled_wcs.ResampledLowLevelWCS` between original and resampled pixel grids. diff --git a/ndcube/tests/test_ndcube_reproject_and_rebin.py b/ndcube/tests/test_ndcube_reproject_and_rebin.py index 5ac395ab8..2af66f6f6 100644 --- a/ndcube/tests/test_ndcube_reproject_and_rebin.py +++ b/ndcube/tests/test_ndcube_reproject_and_rebin.py @@ -143,11 +143,11 @@ def test_rebin(ndcube_3d_l_ln_lt_ectime, bin_shape): 49.99999902, 59.99999831, 69.99999731, 79.99999599], [9.99999999, 19.99999994, 29.99999979, 39.9999995, 49.99999902, 59.99999831, 69.99999731, 79.99999599]]) * u.arcsec - expected_Ty = np.array([[-14.99999996, -14.9999999, -14.99999981, -14.99999969, - -14.99999953, -14.99999934, -14.99999911, -14.99999885], - [-4.99999999, -4.99999998, -4.99999995, -4.9999999, - -4.99999985, -4.99999979, -4.99999971, -4.99999962]]) * u.arcsec - expected_spec = SpectralCoord([1.02e-09], unit=u.m) + expected_Ty = np.array([[-12.49999997, -12.49999993, -12.49999985, -12.49999975, + -12.49999962, -12.49999946, -12.49999926, -12.49999904], + [-2.5 , -2.49999999, -2.49999997, -2.49999995, + -2.49999993, -2.49999989, -2.49999986, -2.49999981]]) * u.arcsec + expected_spec = SpectralCoord([1.11e-09], unit=u.m) expected_time = Time([51544.00104167, 51544.00243056], format="mjd", scale="utc") expected_time.format = "fits" diff --git a/ndcube/wcs/wrappers/resampled_wcs.py b/ndcube/wcs/wrappers/resampled_wcs.py index 3f3c2836a..b24805e26 100644 --- a/ndcube/wcs/wrappers/resampled_wcs.py +++ b/ndcube/wcs/wrappers/resampled_wcs.py @@ -1,3 +1,4 @@ + import numpy as np from astropy.wcs.wcsapi.wrappers.base import BaseWCSWrapper @@ -20,9 +21,59 @@ class ResampledLowLevelWCS(BaseWCSWrapper): axis. If a scalar, the same factor is used for all axes. offset: `int` or `float` or iterable of the same - The location on the underlying pixel grid which corresponds - to zero on the top level pixel grid. If a scalar, the grid will be - shifted by the same amount in all dimensions. + The shift of the lower edge of the 0th pixel (i.e. the pixel + coordinate -0.5) of the resampled grid relative to the lower + edge of the 0th pixel in the original underlying pixel grid, + in units of original pixel widths. (See the schematic in the + Notes section for a graphical example.) If a scalar, the grid + will be shifted by the same amount in all dimensions. + + Notes + ----- + Below is a schematic of how ResampledLowLevelWCS works. The asterisks show the + corners of pixels in a grid before resampling, while the dashes and + pipes show the edges of the resampled pixels. The resampling along the + x-axis has been performed using a factor of 2 and offset of 1, respectively, + while the resampling of the y-axis uses a factor of 3 and offset of 2. + The right column and upper row of numbers along the side and bottom of the + grids denote the edges and centres of the original pixel grid in the original + pixel coordinates. The left column and lower row gives the same locations in + the pixel coordinates of the resampled grid. Note that the resampled pixels + have an (x, y) shape of (2, 3) relative to the original pixel grid. + Also note, the left/lower edge of the 0th pixel in the resampled grid (i.e. pixel + coord -0.5) is shifted relative to the left/lower edge of the original 0th pixel, + and that shift is given by the offset (+1 in the x-axis and +2 along the y-axis), + which is in units of original pixel widths. + + :: + + resampled original + factor=3 + offset=2 + + 0.5 4.5 *-----------*-----------*-----------*-----------* + | | + 2/6 4 | | + | | + 1/6 3.5 * * * * * + | | + 0 3 | | + | | + -1/3 2.5 * * * * * + | | + -2/6 2 | | + | | + -0.5 1.5 *-----------*-----------*-----------*-----------* + | | + -4/6 1 | | + | | + -5/6 0.5 * * * * * + | | + -1 0 | | + | | + -1-1/6 -0.5 * * * * * + -0.5 0 0.5 1 1.5 2 2.5 3 3.5 original pixel indices + -1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 resampled pixel indices: factor=2, offset=1 """ def __init__(self, wcs, factor, offset=0): @@ -42,13 +93,13 @@ def _top_to_underlying_pixels(self, top_pixels): # Convert user-facing pixel indices to the pixel grid of underlying WCS. factor = self._pad_dims(self._factor, top_pixels.ndim) offset = self._pad_dims(self._offset, top_pixels.ndim) - return top_pixels * factor + offset + return (top_pixels + 0.5) * factor - 0.5 + offset def _underlying_to_top_pixels(self, underlying_pixels): # Convert pixel indices of underlying pixel grid to user-facing grid. factor = self._pad_dims(self._factor, underlying_pixels.ndim) offset = self._pad_dims(self._offset, underlying_pixels.ndim) - return (underlying_pixels - offset) / factor + return (underlying_pixels + 0.5 - offset) / factor - 0.5 def _pad_dims(self, arr, ndim): # Pad array with trailing degenerate dimensions. diff --git a/ndcube/wcs/wrappers/tests/test_resampled_wcs.py b/ndcube/wcs/wrappers/tests/test_resampled_wcs.py index 7b7b945e4..4d5781c2d 100644 --- a/ndcube/wcs/wrappers/tests/test_resampled_wcs.py +++ b/ndcube/wcs/wrappers/tests/test_resampled_wcs.py @@ -7,6 +7,7 @@ from astropy import units as u from astropy.coordinates import SkyCoord from astropy.tests.helper import assert_quantity_allclose +from astropy.wcs import WCS from astropy.wcs.wcsapi import HighLevelWCSWrapper from ndcube.wcs.wrappers import ResampledLowLevelWCS @@ -25,8 +26,8 @@ def celestial_wcs(request): Array shape (Numpy order): (2, 15) Pixel Dim Axis Name Data size Bounds - 0 None 15 (-2.5, 12.5) - 1 None 2.33333 (0.3333333333333333, 2.3333333333333335) + 0 None 15 (-1.75, 13.25) + 1 None 2.33333 (0.0, 2.0) World Dim Axis Name Physical Type Units 0 Right Ascension pos.eq.ra deg @@ -47,8 +48,8 @@ def celestial_wcs(request): Array shape (Numpy order): (2.3333333333333335, 15.0) Pixel Dim Axis Name Data size Bounds - 0 None 15 (-2.5, 12.5) - 1 None 2.33333 (0.3333333333333333, 2.3333333333333335) + 0 None 15 (-1.75, 13.25) + 1 None 2.33333 (0.0, 2.0) World Dim Axis Name Physical Type Units 0 Right Ascension pos.eq.ra deg @@ -69,7 +70,8 @@ def test_2d(celestial_wcs): # Upsample along the first pixel dimension and downsample along the second # pixel dimension. - wcs = ResampledLowLevelWCS(celestial_wcs, [0.4, 3]) + factors = [0.4, 3] + wcs = ResampledLowLevelWCS(celestial_wcs, factors) # The following shouldn't change compared to the original WCS assert wcs.pixel_n_dim == 2 @@ -77,83 +79,76 @@ def test_2d(celestial_wcs): assert tuple(wcs.world_axis_physical_types) == ('pos.eq.ra', 'pos.eq.dec') assert tuple(wcs.world_axis_units) == ('deg', 'deg') assert tuple(wcs.pixel_axis_names) == ('', '') - assert tuple(wcs.world_axis_names) == ('Right Ascension', - 'Declination') + assert tuple(wcs.world_axis_names) == ('Right Ascension', 'Declination') assert_equal(wcs.axis_correlation_matrix, np.ones((2, 2))) # Shapes and bounds should be floating-point if needed assert_allclose(wcs.pixel_shape, (15, 7/3)) assert_allclose(wcs.array_shape, (7/3, 15)) - assert_allclose(wcs.pixel_bounds, ((-2.5, 12.5), (1/3, 7/3))) + assert_allclose(wcs.pixel_bounds, ((-1.75, 13.25), (0, 2))) - 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), [1, 2]) + under_pixel, over_pixel = (2.3, 1.2), (5.75+0.6*1.25, 0.2/1.5*0.5) + world = celestial_wcs.pixel_to_world_values(*under_pixel) + + assert_allclose(wcs.pixel_to_world_values(*over_pixel), world) + assert_allclose(wcs.array_index_to_world_values(*over_pixel[::-1]), world) + assert_allclose(wcs.world_to_pixel_values(*world), over_pixel) + assert_allclose(wcs.world_to_array_index_values(*world), + np.around(over_pixel[::-1]).astype(int)) 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]), - np.array([13.8, 14.4])) - assert_allclose(wcs.pixel_to_world_values(*pixel_array), world_array) - assert_allclose(wcs.array_index_to_world_values(*pixel_array[::-1]), world_array) - assert_allclose(wcs.world_to_pixel_values(*world_array), pixel_array) + under_pixel_array = (np.array([2.3, 2.4]), + np.array([4.3, 4.4])) + over_pixel_array = (np.array([over_pixel[0], 5.75+0.8*1.25]), + 1 + np.array([0.3, 0.4]) / 1.5 * 0.5) + world_array = celestial_wcs.pixel_to_world_values(*under_pixel_array) + + assert_allclose(wcs.pixel_to_world_values(*over_pixel_array), world_array) + assert_allclose(wcs.array_index_to_world_values(*over_pixel_array[::-1]), world_array) + assert_allclose(wcs.world_to_pixel_values(*world_array), over_pixel_array) assert_allclose(wcs.world_to_array_index_values(*world_array), - [[4, 4], [2, 2]]) + np.around(np.asarray(over_pixel_array)[::-1]).astype(int)) wcs_hl = HighLevelWCSWrapper(wcs) - celestial = wcs_hl.pixel_to_world(*pixel_scalar) + celestial = wcs_hl.pixel_to_world(*over_pixel) assert isinstance(celestial, SkyCoord) - assert_quantity_allclose(celestial.ra, world_scalar[0] * u.deg) - assert_quantity_allclose(celestial.dec, world_scalar[1] * u.deg) + assert_quantity_allclose(celestial.ra, world[0] * u.deg) + assert_quantity_allclose(celestial.dec, world[1] * u.deg) - celestial = wcs_hl.pixel_to_world(*pixel_array) + celestial = wcs_hl.pixel_to_world(*over_pixel_array) assert isinstance(celestial, SkyCoord) assert_quantity_allclose(celestial.ra, world_array[0] * u.deg) assert_quantity_allclose(celestial.dec, world_array[1] * u.deg) @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) - world_scalar = (4.8, 5.2) - 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]) - - -@pytest.mark.parametrize('celestial_wcs', - ['celestial_2d_ape14_wcs', 'celestial_2d_fitswcs'], + ['celestial_2d_ape14_wcs', + 'celestial_2d_fitswcs'], indirect=True) -def test_offset(celestial_wcs): +@pytest.mark.parametrize(('factor', 'offset', 'over_pixel'), + [(2, 0, (0.9, 1.9)), + (2, 1, (0.4, 1.4)), + ]) +def test_scalar_factor_and_offset(celestial_wcs, factor, offset, over_pixel): celestial_wcs.pixel_bounds = None - offset = 1 - factor = 2 wcs = ResampledLowLevelWCS(celestial_wcs, factor, offset=offset) - - pixel_scalar = (2.3, 4.3) - world_scalar = celestial_wcs.pixel_to_world_values(*[p * factor + offset - for p in pixel_scalar]) - - 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]) + # Define the pixel coord pre-resampled pixel grid corresponding to + # the same location in the resampled grid, as defined in the parameterization. + under_pixel = (2.3, 4.3) + # Get the corresponding world location using original WCS. + world = celestial_wcs.pixel_to_world_values(*under_pixel) + # Confirm resampled WCS maps to and from the same world coordinate to the + # pixel location in the resample pixel coords. + assert_allclose(wcs.pixel_to_world_values(*over_pixel), world) + assert_allclose(wcs.array_index_to_world_values(*over_pixel[::-1]), world) + assert_allclose(wcs.world_to_pixel_values(*world), over_pixel) + assert_allclose(wcs.world_to_array_index_values(*world), + np.around(over_pixel[::-1]).astype(int)) @pytest.mark.parametrize('celestial_wcs', @@ -183,3 +178,110 @@ def test_int_fraction_pixel_shape(celestial_wcs): assert wcs.pixel_shape == (18, 21) for dim in wcs.pixel_shape: assert isinstance(dim, numbers.Integral) + +@pytest.fixture +def four_five_wcs(): + wcs = WCS(naxis=2) + wcs.wcs.ctype = ["HPLN-TAN", "HPLT-TAN"] + wcs.wcs.cdelt = [2, 2] + wcs.wcs.crval = [0, 0] + wcs.wcs.crpix = [1, 1] + + wcs.wcs.set() + wcs.pixel_shape = [4, 5] + return wcs + +@pytest.mark.parametrize(('factor', 'offset', 'expected_over_pixels'), + [([2, 3], [0, 0], np.meshgrid(np.linspace(-0.5, 1.5, 4*2+1), np.linspace(-0.5, 1+1/6, 5*2+1))), + ([2, 3], [1, 2], np.meshgrid(np.linspace(-1, 1, 4*2+1), np.linspace(-1-1/6, 0.5, 5*2+1))), + ]) +def test_resampled_pixel_to_world_values(four_five_wcs, factor, offset, expected_over_pixels): + """ + Notes + ----- + Below shows schematics of the two test cases tested in this test of how ResampledLowLevelWCS. + The asterisks show the corners of pixels in a grid before resampling, while the dashes and + pipes show the edges the resampled pixels. In the first case, the resampled pixels are + obtained by applied resampling factors of 2 along the x-axis and 3 along the y-axis. No + offset is applied to either axis. In the second case, the same resampling factors have + been applied, but an offsets of 1 and 2 pixels have been applied to the x- and y-axes, + respectively. The right column/upper row of numbers along the side of/below the grids denote + the edges and centres of the original pixel grid in the original pixel coordinates. + The left column and lower row gives the same locations in the pixel coordinates of the + resampled grid. + + Test Case 1 + + :: + + resampled original + factor=3 + offset=0 + + 1+1/6 4.5* * * * * + | | | + 1 4 | | | + | | | + 5/6 3.5* * * * * + | | | + 4/6 3 | | | + | | | + | | | + 0.5 2.5*---------*---------*---------*---------* + | | | + 2/6 2 | | | + | | | + 1/6 1.5* * * * * + | | | + 0 1 | | | + | | | + -1/6 0.5* * * * * + | | | + -2/6 0 | | | + | | | + -0.5 -0.5*---------*---------*---------*---------* + -0.5 0 0.5 1 1.5 2 2.5 3 3.5 original pixel indices + -0.5 -0.25 0 0.25 0.5 0.75 1 1.25 1.5 resampled pixel indices: factor=2, offset=0 + + Test Case 2 + + :: + + resampled original + factor=3 + offset=2 + + 0.5 4.5 *-----------*-----------*-----------*-----------* + | | + 2/6 4 | | + | | + 1/6 3.5 * * * * * + | | + 0 3 | | + | | + -1/3 2.5 * * * * * + | | + -2/6 2 | | + | | + -0.5 1.5 *-----------*-----------*-----------*-----------* + | | + -4/6 1 | | + | | + -5/6 0.5 * * * * * + | | + -1 0 | | + | | + -1-1/6 -0.5 * * * * * + -0.5 0 0.5 1 1.5 2 2.5 3 3.5 original pixel indices + -1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 resampled pixel indices: factor=2, offset=1 + """ + wcs = four_five_wcs + # Get world values of original pixel grid. + under_pixels = np.meshgrid(np.arange(-0.5, 4, 0.5), np.arange(-0.5, 5, 0.5)) + expected_world = wcs.pixel_to_world_values(*under_pixels) + + # Resample WCS + new_wcs = ResampledLowLevelWCS(wcs, factor, offset) + # Get expected pixel coords in resampled WCS of same pixel locations as above. + output_world = new_wcs.pixel_to_world_values(*expected_over_pixels) + assert_allclose(np.asarray(output_world), np.asarray(expected_world), atol=1e-15) diff --git a/pytest.ini b/pytest.ini index cd1290262..049235b81 100644 --- a/pytest.ini +++ b/pytest.ini @@ -57,4 +57,4 @@ filterwarnings = # This is raised by the Windows and mac os build for visualization.rst ignore:FigureCanvasAgg is non-interactive, and thus cannot be shown:UserWarning # Oldestdeps from gWCS - ignore:pkg_resources is deprecated as an API:DeprecationWarning + ignore:pkg_resources is deprecated as an API