Subtract blue in the ARVI denominator (#3747)

[brendancol]

Closes #3747

What changed

• _arvi_cpu and _arvi_gpu in xrspatial/multispectral.py added blue in the denominator. ARVI (Kaufman & Tanré 1992) uses rb = red - gamma*(blue - red), which is 2*red - blue for gamma = 1, so the denominator is now nir + 2*red - blue. The numerator was already right.
• The arvi docstring gives the formula, cites the paper, and has updated example output.
• The test fixtures had the old formula baked in. Their expected values are recomputed from the closed form in NumPy float64, not from the kernel. A new test checks the case from the issue on all four backends: nir=0.5, red=0.1, blue=0.05 gives 0.5385, where it used to give 0.4667. It also has a blue > 2*red cell. The qgis_arvi fixture is renamed reference_arvi, since its values no longer come from QGIS.

Behaviour change

Every arvi() result changes. On the Sentinel-2 sample in the docstring the values are about 1.8 to 1.9 times what they were. Anyone comparing against ARVI rasters written by an earlier release will see a difference.

The output is no longer limited to [-1, 1]. Where blue > 2*red (water, haze, snow), rb = 2*red - blue is negative, so the result goes above 1 and gets very large as nir + rb nears zero. For example nir=0.30, red=0.05, blue=0.12 gives 1.143, where the old formula gave 0.615. This is what the published formula does, so the kernel does not clamp or mask it. The docstring now says so.

The zero-denominator case also moves. NaN now comes back where nir + 2*red - blue == 0, so the uint fixture cell (nir=1, red=0, blue=1) is NaN where it used to be 1.

Backends

numpy, cupy, dask+numpy and dask+cupy are all fixed. The two dask paths wrap the CPU and GPU kernels, so the two-line change covers them. I ran the GPU tests locally on a CUDA box.

Not in this PR

examples/user_guide/images/remote_sensing_veg_preview.png has an ARVI panel rendered with the old formula. The notebook stores no outputs and its text still holds, so I left the image alone. It can be regenerated in a follow-up.

Test plan

• pytest xrspatial/tests/test_multispectral.py -k arvi fails before the fix (7 failures) and passes after (11 passed), including test_arvi_gpu[cupy] and test_arvi_gpu[dask+cupy]
• pytest xrspatial/tests/test_multispectral.py xrspatial/tests/test_accessor.py xrspatial/tests/test_validation.py (364 passed)
• Repro script from the issue prints 0.5384615

[brendancol]

PR Review: Subtract blue in the ARVI denominator (#3747)

Reviewed from the issue-3747 worktree at 90749a9. I read both changed files in full and ran the ARVI tests on numpy, dask+numpy, cupy and dask+cupy.

Blockers (must fix before merge)

None.

Suggestions (should fix, not blocking)

• xrspatial/multispectral.py:164-167: the docstring says zero-denominator cells are NaN, but it does not say the output is no longer bounded. With the corrected formula, rb = 2*red - blue goes negative wherever blue > 2*red (water, haze, snow). The result is then above 1, and it gets very large as nir + rb approaches zero. Probe on the PR head:

  nir	red	blue	old result	new result
  0.30	0.05	0.12	0.615	1.143
  0.02	0.03	0.10	0.333	-3.0
  0.10	0.04	0.18	0.556	-2.68e7

  The last row has a denominator of exactly zero in real arithmetic. In float32 it lands on a tiny non-zero value, so the denominator != 0.0 guard at lines 113 and 128 does not catch it and the cell comes back as -2.68e7, not NaN. The old formula could not do this for non-negative reflectance because its denominator was a sum of positive terms. This is what the published formula does, so I would not clamp or mask in the kernel. But the docstring should say it, and the PR's "Behaviour change" section should too, since anyone plotting ARVI with a fixed colour range over water will see it.

• xrspatial/tests/test_multispectral.py: nothing covers the blue > 2*red region. Every fixture cell has blue < red. A small test with one such cell (e.g. the 1.143 row above) would record that values above 1 are expected.

Nits (optional improvements)

• xrspatial/tests/test_multispectral.py:116: the fixture is still called qgis_arvi, and test_arvi_cpu_against_qgis (line 478) still says QGIS, but the values now come from NumPy. The comment explains this, though the names will mislead the next reader. Consider reference_arvi / test_arvi_cpu_against_reference.
• xrspatial/tests/test_multispectral.py:497: test_arvi_matches_kaufman_tanre_3747 runs on numpy only. The four backends are covered through the fixture tests, so this is not a gap in practice. It would still be cheap to run the hand-computed value on dask+numpy as well.

What looks good

• The formula matches Kaufman & Tanré (1992) and the ESA SNAP spec linked in the issue: rb = red - gamma*(blue - red), gamma = 1, arvi = (nir - rb) / (nir + rb). I checked the expansion by hand and against the issue's repro (0.5384615).
• Both kernels are changed (_arvi_cpu line 112, _arvi_gpu line 127). The dask and dask+cupy paths wrap them, so all four backends agree. test_arvi_gpu[cupy] and test_arvi_gpu[dask+cupy] ran for real here and passed.
• The new fixture values were computed from the closed form in float64, independently of the kernel. I recomputed them and they match to float32 precision.
• The uint fixture now includes a true zero-denominator cell, so the NaN guard has a test it did not have before.
• The docstring example output matches what arvi returns on the Sentinel-2 sample.
• The diff is two kernel lines plus docs and tests. No dispatch, dtype or metadata handling changed.

Checklist

• Algorithm matches reference/paper
• All implemented backends produce consistent results
• NaN handling is correct (NaN inputs propagate; exact zero denominator gives NaN)
• Edge cases are covered by tests (missing: blue > 2*red)
• Dask chunk boundaries handled correctly (per-cell op, no overlap needed)
• No premature materialization or unnecessary copies
• Benchmark exists (benchmarks/benchmarks/multispectral.py::Arvi)
• README feature matrix updated (not applicable, no new function or backend change)
• Docstrings present and accurate (missing: unbounded-output note)

[brendancol]

PR Review: follow-up pass on 1042f73

Second pass after the follow-up commit. I re-read both files and re-ran the ARVI tests on all four backends (11 passed, with the cupy and dask+cupy cases running on a real GPU).

Disposition of the first review

• Unbounded output not documented (xrspatial/multispectral.py:169-171): fixed. The docstring now says the result is not limited to [-1, 1] and explains the blue > 2*red case. The PR body's "Behaviour change" section has the same note with a worked value.
• No test for blue > 2*red: fixed. test_arvi_matches_kaufman_tanre_3747 has a second cell (nir=0.30, red=0.05, blue=0.12) that expects 0.32 / 0.28 = 1.143.
• qgis_arvi naming: fixed. The fixture is reference_arvi and the CPU test is test_arvi_cpu_against_reference. No references to the old names are left.
• Hand-computed test on numpy only: fixed. It is parametrized over numpy, dask+numpy, cupy and dask+cupy, with the GPU cases behind cuda_and_cupy_available.

Blockers (must fix before merge)

None.

Suggestions (should fix, not blocking)

None.

Nits (optional improvements)

None new.

Still open, by design

• examples/user_guide/images/remote_sensing_veg_preview.png was rendered with the old formula. The PR body calls this out and leaves it for a follow-up. The notebook itself stores no outputs, so nothing else there is stale.
• The kernel's denominator != 0.0 guard only catches exact float32 zeros. A cell like nir=0.10, red=0.04, blue=0.18 has a zero denominator on paper but returns about -2.7e7 in float32. I am not asking for a change here: an epsilon threshold would be a judgement call that the paper and SNAP do not make, and the docstring now warns about large values near nir + rb = 0.

Checklist

• Algorithm matches reference/paper
• All implemented backends produce consistent results
• NaN handling is correct
• Edge cases are covered by tests
• Dask chunk boundaries handled correctly
• No premature materialization or unnecessary copies
• Benchmark exists or is not needed
• README feature matrix updated (not applicable)
• Docstrings present and accurate