Weight ctis.regrid by the voxel volume to conserve the integral

ctis/_regrid.py

@@ -28,12 +28,13 @@ def regrid(
     depend on spatial position (and vice versa), which is the usual case for a
     CTIS scene.
 
-    Both steps use :func:`named_arrays.regridding.regrid` with
-    ``method="conservative"``, so the total (the sum of ``values_input`` over
-    the resampled axes) is preserved. ``values_input`` is therefore treated as
-    an extensive quantity: a per-voxel total rather than a density. Multiply a
-    spectral radiance by its voxel volume before calling this function if the
-    integral of the radiance is what should be conserved.
+    Each input voxel is weighted by its volume (the wavelength bin width times
+    the solid angle subtended by each field pixel) before the conservative
+    regridding and divided by the output voxel volume afterward, so the
+    *integral* of the field is preserved and ``values_input`` is treated as a
+    density (such as a spectral radiance) rather than a per-voxel total. Both
+    steps use :func:`named_arrays.regridding.regrid` with
+    ``method="conservative"``.
 
     Parameters
     ----------
@@ -104,9 +105,8 @@ def regrid(
         )
 
         # Plot the input and output grids side by side, with wavelength
-        # represented by color. Coarsening conservatively sums voxels, so the
-        # output voxels are brighter; a shared normalization lets one colorbar
-        # describe both panels.
+        # represented by color. The resampling preserves the field density, so
+        # a shared normalization lets one colorbar describe both panels.
         vmax = values_output.max()
         wavelength_min = coordinates_input.wavelength.min()
         wavelength_max = coordinates_input.wavelength.max()
@@ -149,11 +149,17 @@ def regrid(
             ax[1].set_title("output grid")
     """
 
+    axis = (axis_wavelength, *axis_position)
+
+    # weight each input voxel by its volume so the conservative regridding
+    # preserves the integral of the field rather than the per-voxel sum.
+    values = values_input * coordinates_input.volume_cell(axis)
+
     # 1D conservative interpolation along the wavelength axis.
     values = na.regridding.regrid(
         coordinates_input=coordinates_input.wavelength,
         coordinates_output=coordinates_output.wavelength,
-        values_input=values_input,
+        values_input=values,
         axis_input=(axis_wavelength,),
         axis_output=(axis_wavelength,),
         method="conservative",
@@ -169,4 +175,8 @@ def regrid(
         method="conservative",
     )
 
+    # recover the field density on the output grid by dividing out the output
+    # voxel volume.
+    values = values / coordinates_output.volume_cell(axis)
+
     return values

ctis/_regrid_test.py

@@ -61,11 +61,11 @@ def test_regrid(
         axis_position[1]: num_output[axis_position[1]] - 1,
     }
 
-    # the input and output grids span the same volume, so the conservative
-    # resampling preserves the total to within the tolerance of the
-    # perturbation applied by the 2D conservative step.
-    assert np.isclose(
-        result.sum().ndarray,
-        values_input.sum().ndarray,
-        rtol=0.05,
-    )
+    # the input and output grids span the same volume, so the volume-weighted
+    # conservative resampling preserves the integral of the field (to within
+    # the tolerance of the perturbation applied by the 2D conservative step).
+    axis = (axis_wavelength, *axis_position)
+    integral_input = (values_input * coordinates_input.volume_cell(axis)).sum()
+    integral_output = (result * coordinates_output.volume_cell(axis)).sum()
+    ratio = float((integral_output / integral_input).ndarray)
+    assert np.isclose(ratio, 1, rtol=0.05)

pyproject.toml

@@ -17,7 +17,7 @@ classifiers = [
 ]
 dependencies = [
     "astropy",
-    "named-arrays~=2.2",
+    "named-arrays~=2.3",
 ]
 dynamic = ["version"]