PR for getting pixel_edges along with pixel_values in axis_world_coords

changelog/174.feature.rst

@@ -0,0 +1 @@
+- Added a feature to get the pixel_edges from axis_world_coords

docs/ndcube.rst

@@ -434,6 +434,26 @@ is in fact the default.::
              1.49986193e+00]] deg>,
    <Quantity [1.02e-09, 1.04e-09, 1.06e-09, 1.08e-09, 1.10e-09] m>)
 
+By default `~ndcube.NDCube.axis_world_coords` returns the coordinates at the
+center of each pixel. However, the pixel edges can be obtained by setting
+the ``edges`` kwarg to True.
+
+For example,
+  >>> my_cube.axis_world_coords(edges=True)
+  (<Quantity [[0.40006761, 0.40002193, 0.39997624, 0.39993054, 0.39988484],
+            [0.80001604, 0.80000081, 0.79998558, 0.79997035, 0.79995511],
+            [1.19998396, 1.19999919, 1.20001442, 1.20002965, 1.20004489],
+            [1.59993239, 1.59997807, 1.60002376, 1.60006946, 1.60011516]] deg>,
+   <Quantity [[-0.24994347,  0.24998788,  0.74995729,  1.24988864,
+              1.74970582],
+            [-0.24995565,  0.25000006,  0.74999384,  1.24994955,
+              1.74979108],
+            [-0.24995565,  0.25000006,  0.74999384,  1.24994955,
+              1.74979108],
+            [-0.24994347,  0.24998788,  0.74995729,  1.24988864,
+              1.74970582]] deg>,
+   <Quantity [1.01e-09, 1.03e-09, 1.05e-09, 1.07e-09, 1.09e-09, 1.11e-09] m>)
+
 As stated previously, `~ndcube.NDCube` is only written
 to handle single arrays described by single WCS instances.  For cases
 where data is made up of multiple arrays, each described by different

ndcube/ndcube.py

@@ -11,6 +11,7 @@
 from ndcube import utils
 from ndcube.ndcube_sequence import NDCubeSequence
 from ndcube.utils.wcs import wcs_ivoa_mapping
+from ndcube.utils.cube import _pixel_centers_or_edges, _get_dimension_for_pixel
 from ndcube.mixins import NDCubeSlicingMixin, NDCubePlotMixin
 
 
@@ -314,7 +315,7 @@ def world_to_pixel(self, *quantity_axis_list):
             result.append(u.Quantity(world_to_pixel[index], unit=u.pix))
         return result[::-1]
 
-    def axis_world_coords(self, *axes):
+    def axis_world_coords(self, *axes, edges=False):
         """
         Returns WCS coordinate values of all pixels for all axes.
 
@@ -326,6 +327,11 @@ def axis_world_coords(self, *axes):
             of axes for which real world coordinates are desired.
             axes=None implies all axes will be returned.
 
+        edges: `bool`
+            The edges argument helps in returning `pixel_edges`
+            instead of `pixel_values`. Default value is False,
+            which returns `pixel_values`. True return `pixel_edges`
+
         Returns
         -------
         axes_coords: `list` of `astropy.units.Quantity`
@@ -341,6 +347,7 @@ def axis_world_coords(self, *axes):
         cube_dimensions = np.array(self.dimensions.value, dtype=int)
         n_dimensions = cube_dimensions.size
         world_axis_types = self.world_axis_physical_types
+
         # Determine axis numbers of user supplied axes.
         if axes == ():
             int_axes = np.arange(n_dimensions)
@@ -381,23 +388,24 @@ def axis_world_coords(self, *axes):
                 if n_dependent_axes[i] == 1:
                     # Construct pixel quantities in each dimension letting
                     # other dimensions all have 0 pixel value.
-                    quantity_list = [
-                        u.Quantity(np.zeros(cube_dimensions[dependent_axes[i]]),
-                                   unit=u.pix)] * n_dimensions
                     # Replace array in quantity list corresponding to current axis with
                     # np.arange array.
-                    quantity_list[axis] = u.Quantity(np.arange(cube_dimensions[axis]), unit=u.pix)
+                        quantity_list = [
+                        u.Quantity(np.zeros(_get_dimension_for_pixel(cube_dimensions[dependent_axes[i]], edges)),
+                                   unit=u.pix)] * n_dimensions
+                        quantity_list[axis] = u.Quantity(_pixel_centers_or_edges(cube_dimensions[axis], edges), unit=u.pix)
                 else:
                     # If the axis is dependent on another, perform
                     # translations on all dependent axes.
                     # Construct pixel quantities in each dimension letting
                     # other dimensions all have 0 pixel value.
-                    quantity_list = [u.Quantity(np.zeros(tuple(
-                        [cube_dimensions[k] for k in dependent_axes[i]])),
-                        unit=u.pix)] * n_dimensions
                     # Construct orthogonal pixel index arrays for dependent axes.
+                    quantity_list = [u.Quantity(np.zeros(tuple(
+                        [_get_dimension_for_pixel(cube_dimensions[k], edges) for k in dependent_axes[i]])),
+                         unit=u.pix)] * n_dimensions
+
                     dependent_pixel_quantities = np.meshgrid(
-                        *[np.arange(cube_dimensions[k]) * u.pix
+                        *[_pixel_centers_or_edges(cube_dimensions[k], edges) * u.pix
                           for k in dependent_axes[i]], indexing="ij")
                     for k, axis in enumerate(dependent_axes[i]):
                         quantity_list[axis] = dependent_pixel_quantities[k]

ndcube/tests/test_ndcube.py

@@ -893,6 +893,15 @@ def test_all_world_coords_with_input(test_input, expected):
         np.testing.assert_allclose(all_coords[i].value, expected[i].value)
         assert all_coords[i].unit == expected[i].unit
 
+@pytest.mark.parametrize("test_input,expected", [
+    ((cubem, [2]), u.Quantity([1.01e-09, 1.03e-09, 1.05e-09, 1.07e-09, 1.09e-09], unit=u.m)),
+    ((cubem, ['em']), u.Quantity([1.01e-09, 1.03e-09, 1.05e-09, 1.07e-09, 1.09e-09], unit=u.m))
+    ])
+def test_all_world_coords_with_input_and_kwargs(test_input, expected):
+    all_coords = test_input[0].axis_world_coords(*test_input[1], **{"edges":True})
+    for i in range(len(all_coords)):
+        np.testing.assert_allclose(all_coords[i].value, expected[i].value)
+        assert all_coords[i].unit == expected[i].unit
 
 @pytest.mark.parametrize("test_input,expected", [
     (cubem, (u.Quantity([[0.60002173, 0.59999127, 0.5999608],

ndcube/tests/test_utils_cube.py

@@ -137,3 +137,20 @@ def test_convert_extra_coords_dict_to_input_format_error():
     with pytest.raises(KeyError):
         utils.cube.convert_extra_coords_dict_to_input_format(
             {"time": {"not axis": 0, "value": []}}, missing_axes_none)
+
+@pytest.mark.parametrize("test_input, expected",[
+    ((5, False),np.asarray([0, 1, 2, 3, 4])),
+    ((6, True), np.asarray([-0.5,  0.5,  1.5,  2.5,  3.5,  4.5,  5.5]))
+])
+def test_pixel_centers_or_edges(test_input, expected):
+    output = utils.cube._pixel_centers_or_edges(*test_input)
+    assert isinstance(output, np.ndarray)
+    np.testing.assert_allclose(output, expected)
+
+@pytest.mark.parametrize("test_input, expected",[
+    ((5, False), 5),
+    ((6, True), 7)
+])
+def test_get_dimension_for_pixel(test_input, expected):
+    output = utils.cube._get_dimension_for_pixel(*test_input)
+    assert output == expected

ndcube/utils/cube.py

@@ -157,3 +157,41 @@ def get_axis_number_from_axis_name(axis_name, world_axis_physical_types):
                          "a physical axis type. {0} not in any of {1}".format(
                              axis_name, world_axis_physical_types))
     return axis_index[0]
+
+def _pixel_centers_or_edges(axis_length, edges):
+    """
+    Returns a range of pixel_values or pixel_edges
+
+    Parameters
+    ----------
+    axis_length: `int`
+        The length of the axis
+
+    edges: `bool`
+        Boolean to signify whether pixel_edge or pixel_value requested
+        False stands for pixel_value, while True stands for pixel_edge
+
+    Returns
+    -------
+    `np.ndarray`
+        The axis_values for the given input
+    """
+    if edges is False:
+        axis_values = np.arange(axis_length)
+    else:
+        axis_values = np.arange(-0.5, axis_length+0.5)
+    return axis_values
+
+def _get_dimension_for_pixel(axis_length, edges):
+    """
+    Returns the dimensions for the given edges.
+
+    Parameters
+    ----------
+    axis_length : `int`
+        The length of the axis
+    edges : `bool`
+        Boolean to signify whether pixel_edge or pixel_value requested
+        False stands for pixel_value, while True stands for pixel_edge
+    """
+    return axis_length+1 if edges else axis_length

ndcube/utils/wcs.py

@@ -209,7 +209,7 @@ def _wcs_slicer(wcs, missing_axes, item):
             new_wcs = wcs.slice(item_)
             for i, it in enumerate(item_checked):
                 if isinstance(it, int):
-                    missing_axis[i] = True
+                    missing_axes[i] = True
     else:
         raise NotImplementedError("Slicing FITS-WCS by {0} not supported.".format(type(item)))
     # returning the reverse list of missing axis as in the item here was reverse of