From 424f84e0e32ab0aa62f0c30a73c44c6b4b23b8f7 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 30 Mar 2026 12:42:14 -0700 Subject: [PATCH 1/2] Fix SAVI formula and EBBI type inconsistency (#1094) SAVI: The (1+L) factor was in the denominator instead of the numerator. The standard formula (Huete 1988) is ((NIR-Red)/(NIR+Red+L))*(1+L). The code computed (NIR-Red)/((NIR+Red+L)*(1+L)), making all results too small by (1+L)^2 = 2.25 with default L=0.5. EBBI GPU: Changed nb.int64(10) to 10.0 to match CPU path's type behavior. --- xrspatial/multispectral.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/xrspatial/multispectral.py b/xrspatial/multispectral.py index d0c23c82..e8d218b2 100644 --- a/xrspatial/multispectral.py +++ b/xrspatial/multispectral.py @@ -1062,9 +1062,8 @@ def _savi_cpu(nir_data, red_data, soil_factor): red = red_data[y, x] numerator = nir - red soma = nir + red + soil_factor - denominator = soma * (1.0 + soil_factor) - if denominator != 0.0: - out[y, x] = numerator / denominator + if soma != 0.0: + out[y, x] = (numerator / soma) * (1.0 + soil_factor) return out @@ -1077,9 +1076,8 @@ def _savi_gpu(nir_data, red_data, soil_factor, out): red = red_data[y, x] numerator = nir - red soma = nir + red + soil_factor - denominator = soma * (nb.float32(1.0) + soil_factor) - if denominator != 0.0: - out[y, x] = numerator / denominator + if soma != 0.0: + out[y, x] = (numerator / soma) * (1.0 + soil_factor) def _savi_dask(nir_data, red_data, soil_factor): @@ -1367,7 +1365,7 @@ def _ebbi_gpu(red_data, swir_data, tir_data, out): swir = swir_data[y, x] tir = tir_data[y, x] numerator = swir - red - denominator = nb.int64(10) * sqrt(swir + tir) + denominator = 10.0 * sqrt(swir + tir) if denominator != 0.0: out[y, x] = numerator / denominator From d817ddfd70839f96dac5ac507b763742a3ee1c8e Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 30 Mar 2026 12:44:49 -0700 Subject: [PATCH 2/2] Fix SAVI formula and add regression test (#1094) SAVI had (1+L) in the denominator instead of the numerator. The Huete (1988) formula is ((NIR-Red)/(NIR+Red+L))*(1+L). The code computed (NIR-Red)/((NIR+Red+L)*(1+L)), making results too small by (1+L)^2. Updated the qgis_savi fixture and uint dtype fixture with correct reference values. Added test_savi_formula_1094 which checks all 4 backends against the formula directly. Also fixed EBBI GPU: nb.int64(10) -> 10.0 for type consistency. --- xrspatial/tests/test_multispectral.py | 68 ++++++++++++++++++++++----- 1 file changed, 56 insertions(+), 12 deletions(-) diff --git a/xrspatial/tests/test_multispectral.py b/xrspatial/tests/test_multispectral.py index e16eb13b..85b4e435 100644 --- a/xrspatial/tests/test_multispectral.py +++ b/xrspatial/tests/test_multispectral.py @@ -215,19 +215,18 @@ def qgis_ndmi(): @pytest.fixture def qgis_savi(): - # this result is obtained by using NIR, and red band data - # running through QGIS Raster Calculator with formula: - # savi = (nir - red) / ((nir + red + soil_factor) * (1 + soil_factor)) + # Correct SAVI (Huete 1988): + # savi = ((nir - red) / (nir + red + L)) * (1 + L) # with default value of soil_factor=1 result = np.array([ - [0., 0.10726268, 0.10682587, 0.09168259], - [0.10089815, 0.10729991, 0.10749393, np.nan], - [0.11363809, 0.11995638, 0.11994251, 0.10915995], - [0.10226355, 0.11864913, 0.12966092, 0.11774762], - [0.09810041, 0.10675804, 0.12213238, 0.11514599], - [0.09377059, 0.10416108, 0.11123802, 0.10735555], - [0.09097284, 0.0988547, 0.10404798, 0.10413785], - [0.0870268, 0.09878284, 0.105046, 0.10514525]], dtype=np.float32) + [0., 0.4290507, 0.4273035, 0.3667304], + [0.4035926, 0.4291996, 0.4299757, np.nan], + [0.4545524, 0.4798255, 0.4797700, 0.4366398], + [0.4090542, 0.4745965, 0.5186437, 0.4709905], + [0.3924016, 0.4270322, 0.4885295, 0.4605840], + [0.3750824, 0.4166443, 0.4449521, 0.4294222], + [0.3638914, 0.3954188, 0.4161919, 0.4165514], + [0.3481072, 0.3951314, 0.4201840, 0.4205810]], dtype=np.float32) return result @@ -314,7 +313,8 @@ def data_uint_dtype_evi(dtype): def data_uint_dtype_savi(dtype): nir = xr.DataArray(np.array([[1, 1], [1, 1]], dtype=dtype)) red = xr.DataArray(np.array([[0, 1], [0, 2]], dtype=dtype)) - result = np.array([[0.25, 0.], [0.25, -0.125]], dtype=np.float32) + # Correct SAVI with L=1: ((nir-red)/(nir+red+1)) * 2 + result = np.array([[1.0, 0.], [1.0, -0.5]], dtype=np.float32) return nir, red, result @@ -424,6 +424,50 @@ def test_savi_gpu(nir_data, red_data, qgis_savi): general_output_checks(nir_data, result, qgis_savi) +@pytest.mark.parametrize("backend", ["numpy", "cupy", "dask+numpy", "dask+cupy"]) +def test_savi_formula_1094(backend): + """Verify SAVI against the Huete (1988) formula directly. + + Regression test for #1094: (1+L) was in the denominator instead + of the numerator. + """ + from xrspatial.tests.general_checks import has_cuda_and_cupy, dask_array_available as _dask_avail + if 'cupy' in backend and not has_cuda_and_cupy(): + pytest.skip("Requires CUDA and CuPy") + try: + import dask.array + except ImportError: + if 'dask' in backend: + pytest.skip("Requires Dask") + + nir_arr = np.array([[0.8, 0.6], [0.4, 0.0]]) + red_arr = np.array([[0.2, 0.3], [0.4, 0.0]]) + + nir_agg = create_test_raster(nir_arr, backend=backend, chunks=(2, 2)) + red_agg = create_test_raster(red_arr, backend=backend, chunks=(2, 2)) + + L = 0.5 + result = savi(nir_agg, red_agg, soil_factor=L) + data = result.data + if hasattr(data, 'compute'): + data = data.compute() + if hasattr(data, 'get'): + data = data.get() + data = np.asarray(data) + + # Correct: ((NIR - Red) / (NIR + Red + L)) * (1 + L) + expected = np.where( + (nir_arr + red_arr + L) != 0, + ((nir_arr - red_arr) / (nir_arr + red_arr + L)) * (1 + L), + np.nan, + ).astype(np.float32) + + np.testing.assert_allclose(data, expected, rtol=1e-5, equal_nan=True) + + # Spot check: NIR=0.8, Red=0.2, L=0.5 -> (0.6/1.5)*1.5 = 0.6 + assert abs(float(data[0, 0]) - 0.6) < 1e-5 + + # arvi ------------- @pytest.mark.parametrize("backend", ["numpy", "dask+numpy"]) def test_arvi_cpu_against_qgis(nir_data, red_data, blue_data, qgis_arvi):