From ecfdeb5fc7e31401d7719c2eb148e1b0a70fca06 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 28 Jul 2026 12:00:33 -0600 Subject: [PATCH] Weight `ctis.regrid` by the voxel volume to conserve the integral Multiply each input voxel by its volume (via `SpectralPositionalVectorArray.volume_cell`) before the two-stage conservative regridding and divide by the output voxel volume afterward, so `ctis.regrid` preserves the integral of the field and treats `values_input` as a density (such as a spectral radiance) rather than a per-voxel total. Update the test to check integral conservation and bump the named-arrays pin to ~=2.3. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- ctis/_regrid.py | 30 ++++++++++++++++++++---------- ctis/_regrid_test.py | 16 ++++++++-------- pyproject.toml | 2 +- 3 files changed, 29 insertions(+), 19 deletions(-) diff --git a/ctis/_regrid.py b/ctis/_regrid.py index 4d6560f..a7d7d68 100644 --- a/ctis/_regrid.py +++ b/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 diff --git a/ctis/_regrid_test.py b/ctis/_regrid_test.py index 18cf4d7..55b5711 100644 --- a/ctis/_regrid_test.py +++ b/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) diff --git a/pyproject.toml b/pyproject.toml index 9874bce..dfa2d13 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,7 +17,7 @@ classifiers = [ ] dependencies = [ "astropy", - "named-arrays~=2.2", + "named-arrays~=2.3", ] dynamic = ["version"]