Describe the bug
xrspatial.resample.resample(..., method='cubic') calls scipy.ndimage.map_coordinates(data, coords, order=3, mode='nearest') without prefilter=False. With order=3, scipy runs an internal spline_filter before sampling. The internal prefilter honors the mode argument since scipy 1.6, but the boundary transient of the IIR prefilter still depends on how the array is presented to it. A full unchunked array sees one global prefilter pass; each dask chunk sees its own prefilter pass on the (overlapped) chunk view. Even with overlap depth 10, the implicit per-block prefilter on the NaN-aware filled and weights arrays gives values that differ from the eager pass by ~1e-6 near chunk seams.
The current dask-vs-numpy parity test for cubic uses atol=1e-5, which hides the difference.
Expected behavior
A cubic resample of a smooth surface (e.g. a degree-2 polynomial) should match between the eager numpy path and the chunked dask+numpy path to within float64 round-off (~1e-10 atol).
To reproduce
import numpy as np
import xarray as xr
from xrspatial.resample import resample
y, x = np.mgrid[0:30, 0:30].astype(np.float64)
f = x*x + y*y - x*y
da_np = xr.DataArray(f, dims=('y', 'x'),
coords={'y': np.arange(30), 'x': np.arange(30)},
attrs={'res': (1.0, 1.0)})
import dask.array as da
da_dk = da_np.copy()
da_dk.data = da.from_array(da_np.data, chunks=(8, 8))
out_np = resample(da_np, scale_factor=0.5, method='cubic').values
out_dk = resample(da_dk, scale_factor=0.5, method='cubic').values
print(np.max(np.abs(out_np - out_dk))) # ~1e-6, want ~1e-10
Fix
Call scipy.ndimage.spline_filter(data, order=3, mode='nearest') once per block (and per array - filled and weights get separate prefilter passes for the NaN-aware case), then call map_coordinates(..., prefilter=False) on the prefiltered coefficients. Mirror the change in the cupy path with cupyx.scipy.ndimage.spline_filter. Skip the explicit prefilter when order < 2; bilinear order=1 prefilter is a no-op.
Why it matters
Chunk-seam parity for cubic resample is a correctness property users rely on when porting eager prototypes to chunked dask jobs. Tightening the parity bound from 1e-5 to 1e-10 in the test catches future regressions in prefilter handling.
Describe the bug
xrspatial.resample.resample(..., method='cubic')callsscipy.ndimage.map_coordinates(data, coords, order=3, mode='nearest')withoutprefilter=False. Withorder=3, scipy runs an internalspline_filterbefore sampling. The internal prefilter honors themodeargument since scipy 1.6, but the boundary transient of the IIR prefilter still depends on how the array is presented to it. A full unchunked array sees one global prefilter pass; each dask chunk sees its own prefilter pass on the (overlapped) chunk view. Even with overlap depth 10, the implicit per-block prefilter on the NaN-awarefilledandweightsarrays gives values that differ from the eager pass by ~1e-6 near chunk seams.The current dask-vs-numpy parity test for cubic uses
atol=1e-5, which hides the difference.Expected behavior
A cubic resample of a smooth surface (e.g. a degree-2 polynomial) should match between the eager numpy path and the chunked dask+numpy path to within float64 round-off (~1e-10 atol).
To reproduce
Fix
Call
scipy.ndimage.spline_filter(data, order=3, mode='nearest')once per block (and per array -filledandweightsget separate prefilter passes for the NaN-aware case), then callmap_coordinates(..., prefilter=False)on the prefiltered coefficients. Mirror the change in the cupy path withcupyx.scipy.ndimage.spline_filter. Skip the explicit prefilter whenorder < 2; bilinearorder=1prefilter is a no-op.Why it matters
Chunk-seam parity for cubic resample is a correctness property users rely on when porting eager prototypes to chunked dask jobs. Tightening the parity bound from 1e-5 to 1e-10 in the test catches future regressions in prefilter handling.