diff --git a/xrspatial/resample.py b/xrspatial/resample.py index 33ff20c2..f7b9b623 100644 --- a/xrspatial/resample.py +++ b/xrspatial/resample.py @@ -39,8 +39,9 @@ # Overlap depth (input pixels) each interpolation kernel needs from # neighbouring chunks when processing dask arrays. Cubic requires extra # depth because the B-spline prefilter is a global IIR filter whose -# boundary transient decays as ~0.268^n; depth 10 brings it to machine -# epsilon. +# boundary transient decays as ~0.268^n; depth 10 brings the boundary +# transient to ~1e-7 (sufficient for float32 output and well under +# typical interpolation error). _INTERP_DEPTH = {'nearest': 1, 'bilinear': 1, 'cubic': 10} # Approximate working-set size per output cell for the eager backends: @@ -191,7 +192,11 @@ def _nan_aware_interp_np(data, out_h, out_w, order): z_data = _scipy_map_coords(filled, coords, order=order, mode='nearest') z_wt = _scipy_map_coords(weights, coords, order=order, mode='nearest') - result = np.where(z_wt > 0.01, + # Gate on majority weight: an output pixel is valid only when more + # than half of the resampling kernel weight came from valid input + # pixels. This rejects pixels lit only by cubic-kernel sidelobes + # leaking small positive weight from a single neighbour. + result = np.where(z_wt > 0.5, z_data / np.maximum(z_wt, 1e-10), np.nan) return result.reshape(out_h, out_w) @@ -223,7 +228,8 @@ def _nan_aware_interp_cupy(data, out_h, out_w, order): z_data = _cupy_map_coords(filled, coords, order=order, mode='nearest') z_wt = _cupy_map_coords(weights, coords, order=order, mode='nearest') - result = cupy.where(z_wt > 0.01, + # Majority-weight gate (see _nan_aware_interp_np for rationale). + result = cupy.where(z_wt > 0.5, z_data / cupy.maximum(z_wt, 1e-10), cupy.nan) return result.reshape(out_h, out_w) @@ -625,7 +631,8 @@ def _interp_block_np(block, global_in_h, global_in_w, weights = (~mask).astype(np.float64) z_data = _scipy_map_coords(filled, coords, order=order, mode='nearest') z_wt = _scipy_map_coords(weights, coords, order=order, mode='nearest') - result = np.where(z_wt > 0.01, + # Majority-weight gate (see _nan_aware_interp_np for rationale). + result = np.where(z_wt > 0.5, z_data / np.maximum(z_wt, 1e-10), np.nan) return result.reshape(target_h, target_w).astype(np.float32) @@ -667,7 +674,8 @@ def _interp_block_cupy(block, global_in_h, global_in_w, weights = (~mask).astype(cupy.float64) z_data = _cupy_map_coords(filled, coords, order=order, mode='nearest') z_wt = _cupy_map_coords(weights, coords, order=order, mode='nearest') - result = cupy.where(z_wt > 0.01, + # Majority-weight gate (see _nan_aware_interp_np for rationale). + result = cupy.where(z_wt > 0.5, z_data / cupy.maximum(z_wt, 1e-10), cupy.nan) return result.reshape(target_h, target_w).astype(cupy.float32) diff --git a/xrspatial/tests/test_resample.py b/xrspatial/tests/test_resample.py index 862f2dd9..a3f9333f 100644 --- a/xrspatial/tests/test_resample.py +++ b/xrspatial/tests/test_resample.py @@ -682,6 +682,12 @@ def test_target_resolution_list_form(self, grid_8x8): out_tuple.values, out_list.values, atol=1e-6 ) + def test_target_resolution_matches_scale_factor(self, grid_8x8): + # target_resolution=(2.0, 2.0) should match scale_factor=0.5. + a = resample(grid_8x8, target_resolution=(2.0, 2.0), method='nearest') + b = resample(grid_8x8, scale_factor=0.5, method='nearest') + np.testing.assert_allclose(a.values, b.values) + # --------------------------------------------------------------------------- # Memory guard (#1295) @@ -735,6 +741,136 @@ def test_dask_path_skips_guard(self, grid_4x4): assert out.shape == (400, 400) +# --------------------------------------------------------------------------- +# Coordinate orientation (T-3) +# --------------------------------------------------------------------------- + +class TestCoordinateOrientation: + """Resampling must preserve the y-axis orientation of the input.""" + + def test_descending_y_preserved(self): + # 8x8 grid with strong vertical asymmetry: row 0 = 100, row 7 = 0. + data = np.zeros((8, 8), dtype=np.float32) + for r in range(8): + data[r, :] = (7 - r) * (100.0 / 7.0) + # Descending y: 10, 9, ..., 3. x ascending. + y = np.arange(10, 2, -1, dtype=np.float64) + x = np.arange(8, dtype=np.float64) + agg = xr.DataArray( + data, dims=('y', 'x'), coords={'y': y, 'x': x}, + attrs={'res': (1.0, 1.0)}, + ) + + out = resample(agg, scale_factor=0.5, method='average') + + # Output should have 4 y-coords, still descending. + assert out.shape == (4, 4) + assert np.all(np.diff(out.y.values) < 0), ( + "output y must remain descending" + ) + + # Cell-centered y-coords for a 4-pixel output spanning the same + # extent as input y in [2.5, 10.5]: centers at 9.5, 7.5, 5.5, 3.5. + np.testing.assert_allclose( + out.y.values, [9.5, 7.5, 5.5, 3.5], atol=1e-6 + ) + + # Top row of output should average the (high-value) top of the + # input, bottom row should average the (low-value) bottom. No + # vertical flip. + assert out.values[0, 0] > out.values[-1, 0], ( + "top of output should be greater than bottom (no flip)" + ) + # Block-average of rows 0 and 1 = mean of 100 and ~85.7 = ~92.86. + np.testing.assert_allclose(out.values[0, 0], 92.857143, atol=1e-4) + # Bottom block averages rows 6 and 7 = mean of ~14.29 and 0 = ~7.14. + np.testing.assert_allclose(out.values[-1, 0], 7.142857, atol=1e-4) + + +# --------------------------------------------------------------------------- +# Inf input behavior (T-6) +# --------------------------------------------------------------------------- + +class TestInfHandling: + """Lock in the current behavior: infs propagate, they are not NaN-coerced. + + The NaN-aware interpolation path treats infs the same as any other + finite value (they pass through ``np.isnan`` checks), so the + interpolation kernel includes them in its weighted sum and the + aggregate kernel includes them in min / max / mean. Output values + that depend on an inf input cell come back as +/- inf. + """ + + @pytest.fixture + def grid_with_infs(self): + data = np.ones((6, 6), dtype=np.float32) + data[1, 1] = np.inf + data[4, 4] = -np.inf + return create_test_raster(data, attrs={'res': (1.0, 1.0)}) + + def test_bilinear_propagates_inf(self, grid_with_infs): + out = resample(grid_with_infs, scale_factor=0.5, method='bilinear') + # +inf and -inf both reach the output, no NaN coercion. + assert np.isposinf(out.values).any() + assert np.isneginf(out.values).any() + assert not np.isnan(out.values).any() + + def test_average_propagates_inf(self, grid_with_infs): + out = resample(grid_with_infs, scale_factor=0.5, method='average') + assert np.isposinf(out.values).any() + assert np.isneginf(out.values).any() + assert not np.isnan(out.values).any() + + def test_nearest_passes_inf_through(self, grid_with_infs): + out = resample(grid_with_infs, scale_factor=0.5, method='nearest') + # Nearest only picks one input cell per output pixel, so at + # least one infinite input survives in the output. + assert np.isinf(out.values).any() + + +# --------------------------------------------------------------------------- +# 4D input rejection message (T-7) +# --------------------------------------------------------------------------- + +class TestDimensionRejection: + def test_4d_input_error_message(self): + # 3D (band, y, x) is now valid; 4D should still be rejected. + data = np.zeros((2, 2, 4, 4), dtype=np.float32) + agg = xr.DataArray( + data, dims=('time', 'band', 'y', 'x'), + coords={ + 'time': [0, 1], + 'band': [1, 2], + 'y': np.arange(4, dtype=np.float64), + 'x': np.arange(4, dtype=np.float64), + }, + attrs={'res': (1.0, 1.0)}, + ) + with pytest.raises(ValueError, match="must be 2D"): + resample(agg, scale_factor=0.5, method='nearest') + + +# --------------------------------------------------------------------------- +# Bit reproducibility (T-9) +# --------------------------------------------------------------------------- + +class TestReproducibility: + """Repeat calls with identical inputs must produce identical outputs.""" + + @pytest.fixture + def random_grid(self): + data = np.random.RandomState(1471).rand(20, 20).astype(np.float32) + return create_test_raster(data, attrs={'res': (1.0, 1.0)}) + + @pytest.mark.parametrize('method', ['bilinear', 'cubic', 'average']) + def test_bit_identical_repeat(self, random_grid, method): + out1 = resample(random_grid, scale_factor=0.5, method=method) + out2 = resample(random_grid, scale_factor=0.5, method=method) + assert np.array_equal(out1.values, out2.values), ( + f"{method} resample is not bit-identical across runs" + ) + + # --------------------------------------------------------------------------- # Inlined dask aggregate kernel (#1463) # --------------------------------------------------------------------------- @@ -943,20 +1079,3 @@ def test_explicit_nodata_overrides_attr(self): # Without override the attr would say -999 (no match) and -1 would # leak through; with override the top-left output pixel is NaN. assert np.isnan(out.values[0, 0]) - - -# --------------------------------------------------------------------------- -# target_resolution as 2-tuple (issue #1466) -# --------------------------------------------------------------------------- - -class TestTargetResolutionTuple: - def test_tuple_resolution_independent_axes(self, grid_8x8): - # 8x8 grid with res=(1, 1) -> target (2, 4) -> output (4, 2). - out = resample(grid_8x8, target_resolution=(2.0, 4.0)) - assert out.shape == (4, 2) - - def test_tuple_resolution_matches_scale_factor(self, grid_8x8): - # target_resolution=(2.0, 2.0) should match scale_factor=0.5. - a = resample(grid_8x8, target_resolution=(2.0, 2.0), method='nearest') - b = resample(grid_8x8, scale_factor=0.5, method='nearest') - np.testing.assert_allclose(a.values, b.values)