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

ndcube/ndcube.py

@@ -314,7 +314,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, **kwargs):
         """
         Returns WCS coordinate values of all pixels for all axes.
 
@@ -325,6 +325,10 @@ def axis_world_coords(self, *axes):
             `~ndcube.NDCube.world_axis_physical_types`
             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
         -------
@@ -341,6 +345,12 @@ 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
+        # Extract the kwarg `edges` if present.
+        # True signifies pixel_edges
+        # False signifies pixel_values
+        # Default value is False
+        edges = kwargs.get("edges",False)
+
         # Determine axis numbers of user supplied axes.
         if axes == ():
             int_axes = np.arange(n_dimensions)
@@ -381,24 +391,46 @@ 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)
+
+                    # Check for the edges arguments
+                    if edges:
+                        # Setting start and stop ranges for `pixel_edges`
+                        start = -0.5
+                        stop = cube_dimensions[axis] + 0.5
+                        quantity_list = [
+                        u.Quantity(np.zeros(cube_dimensions[dependent_axes[i]]+1),
+                                   unit=u.pix)] * n_dimensions
+                        quantity_list[axis] = u.Quantity(np.arange(start,stop), unit=u.pix)
+                    else:
+                        quantity_list = [
+                        u.Quantity(np.zeros(cube_dimensions[dependent_axes[i]]),
+                                   unit=u.pix)] * n_dimensions
+                        quantity_list[axis] = u.Quantity(np.arange(cube_dimensions[axis]), 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(
+                    # Construct orthogonal pixel index arrays for dependent axes.
+                    # Check for edges arguments
+                    if edges:
+                        quantity_list = [u.Quantity(np.zeros(tuple(
+                        [cube_dimensions[k]+1 for k in dependent_axes[i]])),
+                        unit=u.pix)] * n_dimensions
+
+                        dependent_pixel_quantities = np.meshgrid(
+                            *[np.arange(-0.5, cube_dimensions[k] + 0.5) * u.pix
+                              for k in dependent_axes[i]], indexing="ij")
+                    else:
+                        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.
-                    dependent_pixel_quantities = np.meshgrid(
-                        *[np.arange(cube_dimensions[k]) * u.pix
-                          for k in dependent_axes[i]], indexing="ij")
+
+                        dependent_pixel_quantities = np.meshgrid(
+                            *[np.arange(cube_dimensions[k]) * 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]
                 # Perform wcs translation

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/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