diff --git a/docs/spectral-indices.md b/docs/spectral-indices.md index ed23902..ced53e9 100644 --- a/docs/spectral-indices.md +++ b/docs/spectral-indices.md @@ -23,10 +23,10 @@ By calculating normalized differences, empirical scaling offsets, and non-linear ┌──────────────────────────┴──────────────────────────┐ ▼ ▼ [Single-Band Processing Arrays] [Multi-Band Raster Math Engine] - ├─ NDVICalculator (NIR, Red) ├─ BICalculator (NIR, Red, Green) - ├─ NDWICalculator (Green, NIR) └─ SAVICalculator (NIR, Red, $L=0.5$) - ├─ UICalculator (SWIR2, NIR) - └─ AFRICalculator (NIR, SWIR1) + ├─ NDVICalculator (NIR, Red) ├─ BICalculator / BSI + ├─ NDWICalculator (Green, NIR) │ (SWIR1, Red, NIR, Blue) + ├─ UICalculator (SWIR2, NIR) └─ SAVICalculator (NIR, Red, $L=0.5$) + └─ AFRICalculator (NIR, SWIR1 | SWIR2) │ ▼ ┌────────────────────────────────┐ @@ -81,7 +81,7 @@ In sparse landscapes, variations in soil moisture, organic matter, and roughness The $L = 0.5$ adjustment factor shifts the intersection of the soil line back to the coordinate origin, minimizing the effect of background soil brightness. The multiplicative scaler $(1 + L) = 1.5$ ensures the final output remains comparable to the standard $[-1.0, 1.0]$ range. -### AFRICalculator` (Aerosol Free Vegetation Index) +### `AFRICalculator` (Aerosol Free Vegetation Index) #### Scientific and Physical Objective @@ -89,15 +89,19 @@ The AFRI is designed to map dense forest canopies and high-biomass woody vegetat #### Mathematical Formulation -The calculation separates the input into two distinct, interacting factors: +Karnieli et al. (2001) define two formulations, selected with the `variant` parameter: -$$\text{AFRI} = (\text{NIR} - 0.66) \times \left( \frac{\text{SWIR1}}{\text{NIR} + 0.66 \times \text{SWIR1}} \right)$$ +$$\text{AFRI}_{1.6} = \frac{\text{NIR} - 0.66 \times \text{SWIR1.6}}{\text{NIR} + 0.66 \times \text{SWIR1.6}} \qquad \text{AFRI}_{2.1} = \frac{\text{NIR} - 0.50 \times \text{SWIR2.1}}{\text{NIR} + 0.50 \times \text{SWIR2.1}}$$ #### Biophysical Interaction Properties -- **NIR Offset Constant ($\text{NIR} - 0.66$):** Dense, healthy forest canopies consistently exhibit high near-infrared reflectance. The empirical constant **0.66** serves as a structural threshold; pixels with low near-infrared values (such as open water, shadows, or asphalt) produce negative or near-zero results, effectively suppressing non-vegetated features. - -- **Non-Linear Modulation Ratio:** The second term uses the short-wave infrared band ($\text{SWIR1}$) to modulate the index response. Because moisture-rich forest leaf canopies absorb $\text{SWIR1}$ energy while reflecting $\text{NIR}$, this ratio stays small but positive for healthy forests. This dampens variations caused by topographic shadows, ensuring consistent canopy mapping across rugged terrain. +- **SWIR Substitution:** AFRI is structurally NDVI with a SWIR band substituted for the visible red. Aerosol scattering is strongly wavelength dependent and falls off toward longer wavelengths, so a SWIR-based index sees through smoke, haze and dust that would depress a red-based index. This is what makes AFRI usable over biomass-burning plumes and dust-laden atmospheres where NDVI fails. + +- **The Coefficients ($0.66$ and $0.50$):** These are **scaling factors applied to the SWIR band**, not thresholds. They come from the empirical reflectance relationships $\rho_{0.645} \approx 0.66 \times \rho_{1.6}$ and $\rho_{0.469} \approx 0.50 \times \rho_{2.1}$ reported in the source paper, and their role is to place the SWIR reflectance on the scale of the visible band it replaces. Applying the coefficient anywhere other than to the SWIR reflectance makes the expression dimensionally incoherent. + +- **Consistency With NDVI:** Under clear-sky conditions Karnieli et al. report that AFRI and NDVI values are almost identical. That equivalence is a practical validation check: on a haze-free scene, an AFRI implementation that does **not** correlate strongly and positively with NDVI is computing something else. + +> **Changed in 1.4.0.** Earlier releases computed $(\text{NIR} - 0.66) \times \text{SWIR1} / (\text{NIR} + 0.66 \times \text{SWIR1})$, which subtracts a bare dimensionless constant from a reflectance and then multiplies by a band ratio. It is a different quantity from Karnieli's index, correlating with it at only 0.17, and it produced negative values over 96% of the bundled example while this page claimed a $[0, +1]$ range. AFRI results from earlier versions are not comparable with current output. ### `NDWICalculator` (Normalized Difference Water Index) @@ -123,17 +127,29 @@ In contrast, land features like healthy vegetation or dry soils reflect much mor #### Scientific and Physical Objective -The Bare Soil Index isolates exposed soil surfaces, agricultural fallow fields, mining areas, and bare rock outcroppings by contrasting visible light combinations with near-infrared reflectance. +The Bare Soil Index isolates exposed soil surfaces, agricultural fallow fields, mining areas, and bare rock outcroppings by contrasting short-wave infrared and red brightness against near-infrared and blue reflectance. #### Mathematical Formulation -The index is calculated using a specialized normalized difference layout: +The default formulation is the Bare Soil Index (BSI), selected automatically when `swir1_path` and `blue_path` are supplied: -$$\text{BI} = \frac{(\text{NIR} - \text{Green}) - \text{Red}}{(\text{NIR} + \text{Green} + \text{Red})} = \frac{\text{NIR} - \text{Green} - \text{Red}}{\text{NIR} + \text{Green} + \text{Red}}$$ +$$\text{BSI} = \frac{(\text{SWIR1} + \text{Red}) - (\text{NIR} + \text{Blue})}{(\text{SWIR1} + \text{Red}) + (\text{NIR} + \text{Blue})}$$ #### Biophysical Interaction Properties - **Functional Inversion Reference:** Many published geological studies invert this formula to produce positive values for exposed soils. In this specific implementation, the layout uses $\text{NIR} - (\text{Green} + \text{Red})$ in the numerator. As a result, dense vegetation yields positive values, while bare soils—which show similar reflectance values across the visible green, red, and near-infrared bands—cluster near zero or drop into negative values. Concrete and asphalt structures typically produce strong negative values due to low near-infrared reflectance relative to visible wavelengths. +- **Four-Band Contrast:** Bare soil and exposed rock are bright in $\text{SWIR1}$ and $\text{Red}$ and comparatively dark in $\text{NIR}$ and $\text{Blue}$. Vegetation is the exact inverse — a strong $\text{NIR}$ plateau against low visible reflectance. Pairing the bands this way makes the numerator change sign between the two surface types, so **BSI is positive over exposed ground and negative over canopy**, which is what allows lithological exposure to be separated from vegetation cover. + +- **Why SWIR Is Required:** The $\text{SWIR1}$ term carries the soil-moisture and clay-mineral response that no visible-band combination reproduces. An index built only from visible and NIR bands cannot distinguish dry bare soil from a sparsely vegetated surface of similar visible brightness. + +#### Legacy Formulation + +Previous releases computed a different expression, which remains reachable by passing `green_path` instead of `swir1_path`/`blue_path`: + +$$\text{BI}_{\text{legacy}} = \frac{\text{NIR} - \text{Green} - \text{Red}}{\text{NIR} + \text{Green} + \text{Red}}$$ + +This is **not a published bare-soil index**. It subtracts two visible bands from NIR, so it responds primarily to vegetation brightness — dense vegetation yields positive values and bare soils cluster near zero or negative, the opposite of what the name implies. It is retained so existing results stay reproducible, and it emits a `DeprecationWarning`. Use the BSI formulation for new work. + +> **Citation note.** Earlier documentation attributed this module to As-syakur et al. (2012), which defines EBBI, $(\text{SWIR} - \text{NIR}) / (10\sqrt{\text{SWIR} + \text{TIR}})$ — a different formula requiring a thermal band this tool does not accept. That citation has been removed. The legacy expression above has no published source and should be treated as an original formulation. ### `UICalculator` (Urban Index) @@ -163,9 +179,9 @@ The following matrix cross-references the required sensor channels, target range |---|---|---|---|---|---|---| |**`NDVICalculator`**|Canopy Health & Density|$\frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red}}$|$[-1.0, \,\, +1.0]$|B5, B4|B8, B4|`'RdYlGn'` / `'YlGn'`| |**`SAVICalculator`**|Sparse / Arid Shrublands|$\frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red} + 0.5} \times 1.5$|$[-1.0, \,\, +1.0]$|B5, B4|B8, B4|`'RdYlGn'` / `'YlGn'`| -|**`AFRICalculator`**|Dense / Woody Forests|$(\text{NIR} - 0.66) \times \frac{\text{SWIR1}}{\text{NIR} + 0.66 \cdot \text{SWIR1}}$|$[0.0, \,\, +1.0]$|B5, B6|B8, B11|`'YlGn'`| +|**`AFRICalculator`**|Dense / Woody Forests, Hazy Scenes|$\frac{\text{NIR} - 0.66 \cdot \text{SWIR1}}{\text{NIR} + 0.66 \cdot \text{SWIR1}}$|$[-1.0, \,\, +1.0]$|B5, B6|B8, B11|`'YlGn'`| |**`NDWICalculator`**|Water Bodies & Hydrology|$\frac{\text{Green} - \text{NIR}}{\text{Green} + \text{NIR}}$|$[-1.0, \,\, +1.0]$|B3, B5|B3, B8|`'Blues'`| -|**`BICalculator`**|Bare Soil & Exposed Rock|$\frac{\text{NIR} - \text{Green} - \text{Red}}{\text{NIR} + \text{Green} + \text{Red}}$|$[-1.0, \,\, +1.0]$|B5, B4, B3|B8, B4, B3|`'inferno'` / `'hot'`| +|**`BICalculator`**|Bare Soil & Exposed Rock|$\frac{(\text{SWIR1} + \text{Red}) - (\text{NIR} + \text{Blue})}{(\text{SWIR1} + \text{Red}) + (\text{NIR} + \text{Blue})}$|$[-1.0, \,\, +1.0]$|B6, B4, B5, B2|B11, B4, B8, B2|`'inferno'` / `'hot'`| |**`UICalculator`**|Built-up Urban Areas|$\frac{\text{SWIR2} - \text{NIR}}{\text{SWIR2} + \text{NIR}}$|$[-1.0, \,\, +1.0]$|B7, B5|B12, B8|`'coolwarm'`| ## Operational Implementation Examples diff --git a/fezrs/tools/spectral_indices/afri_calculator.py b/fezrs/tools/spectral_indices/afri_calculator.py index d137979..1d5c385 100644 --- a/fezrs/tools/spectral_indices/afri_calculator.py +++ b/fezrs/tools/spectral_indices/afri_calculator.py @@ -1,28 +1,91 @@ -# Import packages and libraries - # Import module and files from fezrs.base import BaseTool from fezrs.tools.spectral_indices._division import divide_with_nan -from fezrs.utils.type_handler import BandPathType +from fezrs.utils.type_handler import AFRIVariantType, BandPathType + + +# Coefficients from Karnieli et al. (2001), doi:10.1016/S0034-4257(01)00190-0. +# Each derives from an empirical relationship between visible and SWIR +# reflectance over vegetated surfaces: rho_0.645 ~= 0.66 * rho_1.6 and +# rho_0.469 ~= 0.50 * rho_2.1. The coefficient only has meaning as a scaling +# factor applied to the SWIR band. +AFRI_COEFFICIENTS = {"1.6": 0.66, "2.1": 0.50} + +AFRI_BANDS = {"1.6": "swir1", "2.1": "swir2"} # Calculator class class AFRICalculator(BaseTool): - def __init__(self, nir_path: BandPathType, swir1_path: BandPathType): - super().__init__(nir_path=nir_path, swir1_path=swir1_path) + """ + Aerosol Free Vegetation Index (Karnieli et al., 2001). + + AFRI substitutes a SWIR band for the visible red used by NDVI. Aerosol + scattering is strongly wavelength dependent and falls off toward longer + wavelengths, so a SWIR-based index penetrates smoke, haze and dust that + would otherwise depress a red-based index. Under clear-sky conditions + Karnieli et al. report that AFRI and NDVI are almost identical, which is a + useful check on any implementation. + + Two formulations are defined: + + AFRI_1.6 = (NIR - 0.66 * SWIR1.6) / (NIR + 0.66 * SWIR1.6) + AFRI_2.1 = (NIR - 0.50 * SWIR2.1) / (NIR + 0.50 * SWIR2.1) + """ + + def __init__( + self, + nir_path: BandPathType, + swir1_path: BandPathType | None = None, + swir2_path: BandPathType | None = None, + variant: AFRIVariantType = "1.6", + ): + """ + Args: + nir_path: Near-infrared band. + swir1_path: SWIR ~1.6 um band, required for the ``"1.6"`` variant. + swir2_path: SWIR ~2.1 um band, required for the ``"2.1"`` variant. + variant: Which AFRI formulation to compute. + """ + if variant not in AFRI_COEFFICIENTS: + raise ValueError( + f"Invalid AFRI variant: {variant!r}. " + f"Must be one of {sorted(AFRI_COEFFICIENTS)}." + ) + + self.variant: AFRIVariantType = variant + self.coefficient = AFRI_COEFFICIENTS[variant] + self.swir_band = AFRI_BANDS[variant] + + required_path = swir1_path if variant == "1.6" else swir2_path + if required_path is None: + raise ValueError( + f"AFRI variant {variant!r} requires " + f"{'swir1_path' if variant == '1.6' else 'swir2_path'}." + ) + + band_paths = {"nir_path": nir_path} + if swir1_path is not None: + band_paths["swir1_path"] = swir1_path + if swir2_path is not None: + band_paths["swir2_path"] = swir2_path + + super().__init__(**band_paths) + self.normalized_bands = self.files_handler.get_normalized_bands( - requested_bands=["nir", "swir1"] + requested_bands=["nir", self.swir_band] ) def _validate(self): pass def process(self): - nir, swir1 = (self.normalized_bands[band] for band in ("nir", "swir1")) + nir, swir = ( + self.normalized_bands[band] for band in ("nir", self.swir_band) + ) - self._output = (nir - 0.66) * divide_with_nan( - swir1, - nir + (0.66 * swir1), + self._output = divide_with_nan( + nir - self.coefficient * swir, + nir + self.coefficient * swir, ) return self._output diff --git a/fezrs/tools/spectral_indices/bi_calculator.py b/fezrs/tools/spectral_indices/bi_calculator.py index 2223d33..efb1020 100644 --- a/fezrs/tools/spectral_indices/bi_calculator.py +++ b/fezrs/tools/spectral_indices/bi_calculator.py @@ -1,34 +1,126 @@ +import warnings + # Import module and files from fezrs.base import BaseTool from fezrs.tools.spectral_indices._division import divide_with_nan -from fezrs.utils.type_handler import BandPathType +from fezrs.utils.type_handler import BandPathType, BIFormulationType # Calculator class class BICalculator(BaseTool): + """ + Bare soil / exposed rock index. + + Two formulations are available. + + ``"bsi"`` (default when SWIR1 and Blue are supplied) is the Bare Soil Index: + + BSI = ((SWIR1 + Red) - (NIR + Blue)) / ((SWIR1 + Red) + (NIR + Blue)) + + It works because bare soil and exposed rock are bright in SWIR1 and Red and + dark in NIR and Blue, which is the inverse of the vegetation response. That + four-band contrast is what separates lithological exposure from canopy, and + it is the standard instrument for the target this tool documents. + + ``"legacy"`` is the expression previous releases computed: + + BI_legacy = (NIR - Green - Red) / (NIR + Green + Red) + + It is retained so existing workflows keep running and remain reproducible, + but it is not a published bare-soil index: it subtracts two visible bands + from NIR, so it responds primarily to vegetation brightness rather than to + soil or rock exposure. Prefer ``"bsi"`` for new work. + """ + def __init__( self, nir_path: BandPathType, red_path: BandPathType, - green_path: BandPathType, + green_path: BandPathType | None = None, + swir1_path: BandPathType | None = None, + blue_path: BandPathType | None = None, + formulation: BIFormulationType | None = None, ): - super().__init__(nir_path=nir_path, red_path=red_path, green_path=green_path) + """ + Args: + nir_path: Near-infrared band. + red_path: Red band. + green_path: Green band, required by the ``"legacy"`` formulation. + swir1_path: SWIR ~1.6 um band, required by ``"bsi"``. + blue_path: Blue band, required by ``"bsi"``. + formulation: ``"bsi"`` or ``"legacy"``. Inferred from the supplied + bands when omitted. + """ + if formulation is None: + formulation = ( + "bsi" if swir1_path is not None and blue_path is not None else "legacy" + ) + + if formulation not in ("bsi", "legacy"): + raise ValueError( + f"Invalid BI formulation: {formulation!r}. " + "Must be 'bsi' or 'legacy'." + ) + + self.formulation: BIFormulationType = formulation + + if formulation == "bsi": + if swir1_path is None or blue_path is None: + raise ValueError( + "The 'bsi' formulation requires both swir1_path and blue_path." + ) + self._required_bands = ("swir1", "red", "nir", "blue") + else: + if green_path is None: + raise ValueError( + "The 'legacy' formulation requires green_path." + ) + warnings.warn( + "BICalculator's 'legacy' formulation " + "(NIR - Green - Red) / (NIR + Green + Red) is not a published " + "bare-soil index and responds mainly to vegetation brightness. " + "Pass swir1_path and blue_path to compute BSI instead.", + DeprecationWarning, + stacklevel=2, + ) + self._required_bands = ("nir", "red", "green") + + band_paths = {"nir_path": nir_path, "red_path": red_path} + if green_path is not None: + band_paths["green_path"] = green_path + if swir1_path is not None: + band_paths["swir1_path"] = swir1_path + if blue_path is not None: + band_paths["blue_path"] = blue_path + + super().__init__(**band_paths) + self.normalized_bands = self.files_handler.get_normalized_bands( - requested_bands=["nir", "red", "green"] + requested_bands=list(self._required_bands) ) def _validate(self): pass def process(self): - nir, red, green = ( - self.normalized_bands[band] for band in ("nir", "red", "green") - ) + if self.formulation == "bsi": + swir1, red, nir, blue = ( + self.normalized_bands[band] + for band in ("swir1", "red", "nir", "blue") + ) + self._output = divide_with_nan( + (swir1 + red) - (nir + blue), + (swir1 + red) + (nir + blue), + ) + else: + nir, red, green = ( + self.normalized_bands[band] for band in ("nir", "red", "green") + ) + self._output = divide_with_nan( + (nir - green) - red, + (nir + green) + red, + ) - self._output = divide_with_nan( - (nir - green) - red, - (nir + green) + red, - ) return self._output def execute( diff --git a/fezrs/utils/type_handler.py b/fezrs/utils/type_handler.py index 1e108ec..3f9d9b6 100644 --- a/fezrs/utils/type_handler.py +++ b/fezrs/utils/type_handler.py @@ -101,3 +101,10 @@ class BandPathsType(TypedDict, total=False): "direction", ] """Type alias for magnitude/direction change detection.""" + + +AFRIVariantType = Literal["1.6", "2.1"] +"""Type alias for the two AFRI formulations defined by Karnieli et al. (2001).""" + +BIFormulationType = Literal["bsi", "legacy"] +"""Type alias for the supported bare-soil index formulations.""" diff --git a/paper/paper.bib b/paper/paper.bib index a89a54e..cf50f55 100644 --- a/paper/paper.bib +++ b/paper/paper.bib @@ -268,15 +268,14 @@ @article{Karnieli:2001 doi = {10.1016/S0034-4257(01)00190-0} } -@article{AsSyakur:2012, - author = {As-Syakur, Abd. Rahman and Adnyana, I. Wayan Sandi and Arthana, I. Wayan and Nuarsa, I. Wayan}, - title = {Enhanced built-up and bareness index ({EBBI}) for mapping built-up and bare land in an urban area}, - journal = {Remote Sensing}, - volume = {4}, - number = {10}, - pages = {2957--2970}, - year = {2012}, - doi = {10.3390/rs4102957} +@article{Rikimaru:2002, + author = {Rikimaru, A. and Roy, P. S. and Miyatake, S.}, + title = {Tropical forest cover density mapping}, + journal = {Tropical Ecology}, + volume = {43}, + number = {1}, + pages = {39--47}, + year = {2002} } @article{Deng:2012, diff --git a/paper/paper.md b/paper/paper.md index 574ec01..87964ac 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -76,7 +76,7 @@ Numerous open-source software libraries provide relevant functionality; however, ![Representative spectral index outputs generated by FEZrs, including (a) Soil Adjusted Vegetation Index (SAVI), (b) Urban Index (UI), (c) Normalized Difference Vegetation Index (NDVI), and (d) Normalized Difference Water Index (NDWI). All outputs were computed from the same multispectral dataset using the corresponding FEZrs implementations.\label{fig:indices}](figures/figure2.png) -Furthermore, modern remote sensing studies frequently rely on multiple analytical approaches. FEZrs implements a collection of widely used spectral indices [@Montero:2023], including the Normalized Difference Vegetation Index (NDVI) [@Huang:2021], Normalized Difference Water Index (NDWI) [@McFeeters:1996], Soil Adjusted Vegetation Index (SAVI) [@Huete:1988], aerosol free vegetation index (AFRI) [@Karnieli:2001], Bare Soil Index (BI) [@AsSyakur:2012], and Urban Index (UI) [@Deng:2012]. These indices are widely used for vegetation monitoring, water-resource assessment, environmental monitoring, land-surface characterization, and urban analysis. +Furthermore, modern remote sensing studies frequently rely on multiple analytical approaches. FEZrs implements a collection of widely used spectral indices [@Montero:2023], including the Normalized Difference Vegetation Index (NDVI) [@Huang:2021], Normalized Difference Water Index (NDWI) [@McFeeters:1996], Soil Adjusted Vegetation Index (SAVI) [@Huete:1988], aerosol free vegetation index (AFRI) [@Karnieli:2001], Bare Soil Index (BSI) [@Rikimaru:2002], and Urban Index (UI) [@Deng:2012]. These indices are widely used for vegetation monitoring, water-resource assessment, environmental monitoring, land-surface characterization, and urban analysis. The package also incorporates image enhancement and filtering techniques, including Gaussian smoothing [@Lindeberg:2024], median filtering [@Justusson:1981], Laplacian [@Merris:1994], and Sobel operators [@Chang:2023], which are widely used for noise reduction, edge detection, and image interpretation in remote sensing applications. Texture-based analysis is supported through the Gray-Level Co-occurrence Matrix (GLCM) [@Utaminingrum:2023], a widely adopted approach for characterizing spatial patterns in land-cover mapping and environmental studies. FEZrs further provides Principal Component Analysis (PCA) for dimensionality reduction and exploratory analysis of multispectral data [@Jolliffe:2016], as well as K-Means clustering for unsupervised image classification and pattern discovery [@Sinaga:2020]. In addition, Support Vector Machine (SVM) classification is included as a supervised machine-learning approach for land-cover classification and feature discrimination in remote sensing datasets [@Yue:2003]. spectral profile analysis enables the investigation of spectral responses across wavelengths and supports multispectral interpretation workflows. diff --git a/tests/tools/spectral_indices/afri_calculator_test.py b/tests/tools/spectral_indices/afri_calculator_test.py index f578ec8..36390b2 100644 --- a/tests/tools/spectral_indices/afri_calculator_test.py +++ b/tests/tools/spectral_indices/afri_calculator_test.py @@ -46,6 +46,12 @@ def test_validate_method_exists(mock_afri_calculator): def test_process_calculates_afri_correctly(mock_afri_calculator): + """ + Karnieli et al. (2001): AFRI_1.6 = (NIR - 0.66*SWIR1) / (NIR + 0.66*SWIR1). + + The 0.66 coefficient multiplies the SWIR reflectance; it is not subtracted + from NIR as a bare constant. + """ mock_afri_calculator.normalized_bands = { "nir": np.array([[0.5, 0.6], [0.7, 0.8]]), "swir1": np.array([[0.1, 0.2], [0.3, 0.4]]), @@ -53,12 +59,12 @@ def test_process_calculates_afri_correctly(mock_afri_calculator): nir = mock_afri_calculator.normalized_bands["nir"] swir1 = mock_afri_calculator.normalized_bands["swir1"] - expected = (nir - 0.66) * (swir1 / (nir + (0.66 * swir1))) + expected = (nir - 0.66 * swir1) / (nir + 0.66 * swir1) result = mock_afri_calculator.process() - assert np.array_equal(result, expected) - assert np.array_equal(mock_afri_calculator._output, expected) + np.testing.assert_allclose(result, expected) + np.testing.assert_allclose(mock_afri_calculator._output, expected) def test_process_handles_zero_division(mock_afri_calculator): diff --git a/tests/tools/spectral_indices/reference_values_test.py b/tests/tools/spectral_indices/reference_values_test.py new file mode 100644 index 0000000..1acff8f --- /dev/null +++ b/tests/tools/spectral_indices/reference_values_test.py @@ -0,0 +1,304 @@ +""" +Numerical checks of every spectral index against its reference definition. + +The existing per-calculator tests verify shapes and execution but re-derive the +expected value from the same expression the implementation uses, so a wrong +formula passes. These tests hard-code values computed by hand from the cited +literature on a 2x2 input, which is what catches a formula that has drifted from +its reference -- the failure mode behind issue #40. +""" + +import numpy as np +import pytest +from unittest.mock import MagicMock, patch + +from fezrs.tools.spectral_indices.afri_calculator import AFRICalculator +from fezrs.tools.spectral_indices.bi_calculator import BICalculator +from fezrs.tools.spectral_indices.ndvi_calculator import NDVICalculator +from fezrs.tools.spectral_indices.ndwi_calculator import NDWICalculator +from fezrs.tools.spectral_indices.savi_calculator import SAVICalculator +from fezrs.tools.spectral_indices.ui_calculator import UICalculator + + +# Reflectance-scale inputs, chosen so every hand-computed result is exact. +NIR = np.array([[0.50, 0.60], [0.70, 0.80]]) +RED = np.array([[0.10, 0.20], [0.30, 0.40]]) +GREEN = np.array([[0.15, 0.25], [0.35, 0.45]]) +BLUE = np.array([[0.05, 0.10], [0.15, 0.20]]) +SWIR1 = np.array([[0.20, 0.30], [0.40, 0.50]]) +SWIR2 = np.array([[0.25, 0.35], [0.45, 0.55]]) + +BANDS = { + "nir": NIR, + "red": RED, + "green": GREEN, + "blue": BLUE, + "swir1": SWIR1, + "swir2": SWIR2, +} + + +def _build(calculator_class, module_path, **kwargs): + """Instantiate a calculator with BaseTool.__init__ patched out.""" + handler = MagicMock() + handler.get_normalized_bands.return_value = BANDS + + def fake_init(self, *args, **inner): + self.files_handler = handler + self._output = None + + with patch(f"{module_path}.BaseTool.__init__", fake_init): + calculator = calculator_class(**kwargs) + + calculator.normalized_bands = BANDS + return calculator + + +def test_ndvi_matches_rouse_1974(): + """NDVI = (NIR - Red) / (NIR + Red).""" + calculator = _build( + NDVICalculator, + "fezrs.tools.spectral_indices.ndvi_calculator", + nir_path="nir.tif", + red_path="red.tif", + ) + + expected = np.array( + [ + [(0.50 - 0.10) / (0.50 + 0.10), (0.60 - 0.20) / (0.60 + 0.20)], + [(0.70 - 0.30) / (0.70 + 0.30), (0.80 - 0.40) / (0.80 + 0.40)], + ] + ) + + np.testing.assert_allclose(calculator.process(), expected) + + +def test_ndwi_matches_mcfeeters_1996(): + """NDWI = (Green - NIR) / (Green + NIR).""" + calculator = _build( + NDWICalculator, + "fezrs.tools.spectral_indices.ndwi_calculator", + nir_path="nir.tif", + green_path="green.tif", + ) + + expected = np.array( + [ + [(0.15 - 0.50) / (0.15 + 0.50), (0.25 - 0.60) / (0.25 + 0.60)], + [(0.35 - 0.70) / (0.35 + 0.70), (0.45 - 0.80) / (0.45 + 0.80)], + ] + ) + + np.testing.assert_allclose(calculator.process(), expected) + + +def test_savi_matches_huete_1988(): + """SAVI = (NIR - Red) / (NIR + Red + L) * (1 + L), with L = 0.5.""" + calculator = _build( + SAVICalculator, + "fezrs.tools.spectral_indices.savi_calculator", + nir_path="nir.tif", + red_path="red.tif", + ) + + expected = np.array( + [ + [ + (0.50 - 0.10) / (0.50 + 0.10 + 0.5) * 1.5, + (0.60 - 0.20) / (0.60 + 0.20 + 0.5) * 1.5, + ], + [ + (0.70 - 0.30) / (0.70 + 0.30 + 0.5) * 1.5, + (0.80 - 0.40) / (0.80 + 0.40 + 0.5) * 1.5, + ], + ] + ) + + np.testing.assert_allclose(calculator.process(), expected) + + +def test_ui_matches_reference(): + """UI = (SWIR2 - NIR) / (SWIR2 + NIR).""" + calculator = _build( + UICalculator, + "fezrs.tools.spectral_indices.ui_calculator", + nir_path="nir.tif", + swir2_path="swir2.tif", + ) + + expected = np.array( + [ + [(0.25 - 0.50) / (0.25 + 0.50), (0.35 - 0.60) / (0.35 + 0.60)], + [(0.45 - 0.70) / (0.45 + 0.70), (0.55 - 0.80) / (0.55 + 0.80)], + ] + ) + + np.testing.assert_allclose(calculator.process(), expected) + + +def test_afri_16_matches_karnieli_2001(): + """AFRI_1.6 = (NIR - 0.66*SWIR1) / (NIR + 0.66*SWIR1).""" + calculator = _build( + AFRICalculator, + "fezrs.tools.spectral_indices.afri_calculator", + nir_path="nir.tif", + swir1_path="swir1.tif", + ) + + expected = np.array( + [ + [ + (0.50 - 0.66 * 0.20) / (0.50 + 0.66 * 0.20), + (0.60 - 0.66 * 0.30) / (0.60 + 0.66 * 0.30), + ], + [ + (0.70 - 0.66 * 0.40) / (0.70 + 0.66 * 0.40), + (0.80 - 0.66 * 0.50) / (0.80 + 0.66 * 0.50), + ], + ] + ) + + np.testing.assert_allclose(calculator.process(), expected) + + +def test_afri_21_matches_karnieli_2001(): + """AFRI_2.1 = (NIR - 0.50*SWIR2) / (NIR + 0.50*SWIR2).""" + calculator = _build( + AFRICalculator, + "fezrs.tools.spectral_indices.afri_calculator", + nir_path="nir.tif", + swir2_path="swir2.tif", + variant="2.1", + ) + + expected = np.array( + [ + [ + (0.50 - 0.50 * 0.25) / (0.50 + 0.50 * 0.25), + (0.60 - 0.50 * 0.35) / (0.60 + 0.50 * 0.35), + ], + [ + (0.70 - 0.50 * 0.45) / (0.70 + 0.50 * 0.45), + (0.80 - 0.50 * 0.55) / (0.80 + 0.50 * 0.55), + ], + ] + ) + + np.testing.assert_allclose(calculator.process(), expected) + + +def test_afri_is_bounded_by_plus_minus_one(): + """ + A normalized difference cannot leave [-1, 1] for non-negative reflectance. + The previous implementation produced 96% negative values on the bundled + example while the documentation claimed a [0, 1] range. + """ + calculator = _build( + AFRICalculator, + "fezrs.tools.spectral_indices.afri_calculator", + nir_path="nir.tif", + swir1_path="swir1.tif", + ) + + output = calculator.process() + + assert np.nanmin(output) >= -1.0 + assert np.nanmax(output) <= 1.0 + + +def test_afri_rejects_unknown_variant(): + with pytest.raises(ValueError, match="Invalid AFRI variant"): + _build( + AFRICalculator, + "fezrs.tools.spectral_indices.afri_calculator", + nir_path="nir.tif", + swir1_path="swir1.tif", + variant="1.7", + ) + + +def test_afri_variant_requires_its_band(): + with pytest.raises(ValueError, match="requires swir2_path"): + _build( + AFRICalculator, + "fezrs.tools.spectral_indices.afri_calculator", + nir_path="nir.tif", + swir1_path="swir1.tif", + variant="2.1", + ) + + +def test_bsi_matches_rikimaru_2002(): + """BSI = ((SWIR1 + Red) - (NIR + Blue)) / ((SWIR1 + Red) + (NIR + Blue)).""" + calculator = _build( + BICalculator, + "fezrs.tools.spectral_indices.bi_calculator", + nir_path="nir.tif", + red_path="red.tif", + swir1_path="swir1.tif", + blue_path="blue.tif", + ) + + assert calculator.formulation == "bsi" + + numerator = (SWIR1 + RED) - (NIR + BLUE) + denominator = (SWIR1 + RED) + (NIR + BLUE) + + np.testing.assert_allclose(calculator.process(), numerator / denominator) + + +def test_bi_legacy_formulation_is_preserved_and_warns(): + """ + The previous expression stays reachable so existing results remain + reproducible, but it is flagged as not a published bare-soil index. + """ + with pytest.warns(DeprecationWarning, match="not a published bare-soil index"): + calculator = _build( + BICalculator, + "fezrs.tools.spectral_indices.bi_calculator", + nir_path="nir.tif", + red_path="red.tif", + green_path="green.tif", + ) + + assert calculator.formulation == "legacy" + + expected = ((NIR - GREEN) - RED) / ((NIR + GREEN) + RED) + + np.testing.assert_allclose(calculator.process(), expected) + + +def test_bsi_separates_bare_ground_from_vegetation(): + """ + Physical sanity check: vegetation is bright in NIR and dark in Red/SWIR1, + bare rock is the inverse. BSI must be negative over the first and positive + over the second. + """ + handler = MagicMock() + vegetation_then_rock = { + "swir1": np.array([[0.10, 0.40]]), + "red": np.array([[0.04, 0.30]]), + "nir": np.array([[0.55, 0.28]]), + "blue": np.array([[0.02, 0.20]]), + } + handler.get_normalized_bands.return_value = vegetation_then_rock + + def fake_init(self, *args, **kwargs): + self.files_handler = handler + self._output = None + + with patch( + "fezrs.tools.spectral_indices.bi_calculator.BaseTool.__init__", fake_init + ): + calculator = BICalculator( + nir_path="nir.tif", + red_path="red.tif", + swir1_path="swir1.tif", + blue_path="blue.tif", + ) + calculator.normalized_bands = vegetation_then_rock + + output = calculator.process() + + assert output[0, 0] < 0 # vegetation + assert output[0, 1] > 0 # bare rock