Fix Errors in computation in ResampledLowLevelWCS

changelog/857.bugfix.rst

@@ -0,0 +1 @@
+Fix translations in `~ndcube.wcs.wrappers.resampled_wcs.ResampledLowLevelWCS` between original and resampled pixel grids.

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"
 

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.

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,91 +70,85 @@ 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
     assert wcs.world_n_dim == 2
     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)

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