# -*- coding: utf-8 -*-
# Copyright 2022, SERTIT-ICube - France, https://sertit.unistra.fr/
# This file is part of eoreader project
# https://github.com/sertit/eoreader
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
Set of usual optical index.
**Note**: This is easier to manage index as raw functions in a file rather than stored in a class
"""
# Index not snake case
# pylint: disable=C0103
import inspect
import logging
import re
import sys
from functools import wraps
from typing import Callable
import numpy as np
import xarray as xr
from sertit import rasters
from eoreader.bands.bands import OpticalBandNames as obn
from eoreader.utils import EOREADER_NAME
LOGGER = logging.getLogger(EOREADER_NAME)
np.seterr(divide="ignore", invalid="ignore")
def _idx_fct(function: Callable) -> Callable:
"""
Decorator of index functions
"""
@wraps(function)
def _idx_fct_wrapper(bands: dict) -> xr.DataArray:
"""
Index functions wrapper
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
# WARNING: for performance issues, use numpy arrays here to speed up computation !
out_np = function({key: val.data for key, val in bands.items()})
# Take the first band as a template for xarray
first_xda = list(bands.values())[0]
out_xda = first_xda.copy(data=out_np)
out = rasters.set_metadata(out_xda, first_xda, new_name=str(function.__name__))
return out
return _idx_fct_wrapper
def _norm_diff(band_1: xr.DataArray, band_2: xr.DataArray) -> xr.DataArray:
"""
Get normalized difference index between band 1 and band 2:
(band_1 - band_2)/(band_1 + band_2)
Args:
band_1 (xr.DataArray): Band 1
band_2 (xr.DataArray): Band 2
Returns:
xr.DataArray: Normalized Difference between band 1 and band 2
"""
norm = (band_1 - band_2) / (band_1 + band_2)
return norm
[docs]@_idx_fct
def RGI(bands: dict) -> xr.DataArray:
"""
Relative Greenness Index: https://www.indexdatabase.de/db/i-single.php?id=326
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return bands[obn.RED] / bands[obn.GREEN]
[docs]@_idx_fct
def NDVI(bands: dict) -> xr.DataArray:
"""
Normalized Difference Vegetation Index: https://www.indexdatabase.de/db/i-single.php?id=59
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.NIR], bands[obn.RED])
[docs]@_idx_fct
def TCBRI(bands: dict) -> xr.DataArray:
"""
Tasseled Cap Brightness:
https://en.wikipedia.org/wiki/Tasseled_cap_transformation
https://www.indexdatabase.de/db/r-single.php?id=723
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return (
0.3037 * bands[obn.BLUE]
+ 0.2793 * bands[obn.GREEN]
+ 0.4743 * bands[obn.RED]
+ 0.5585 * bands[obn.NIR]
+ 0.5082 * bands[obn.SWIR_1]
+ 0.1863 * bands[obn.SWIR_2]
)
[docs]@_idx_fct
def TCGRE(bands: dict) -> xr.DataArray:
"""
Tasseled Cap Greenness:
https://en.wikipedia.org/wiki/Tasseled_cap_transformation
https://www.indexdatabase.de/db/r-single.php?id=723
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return (
-0.2848 * bands[obn.BLUE]
- 0.2435 * bands[obn.GREEN]
- 0.5436 * bands[obn.RED]
+ 0.7243 * bands[obn.NIR]
+ 0.0840 * bands[obn.SWIR_1]
- 0.1800 * bands[obn.SWIR_2]
)
[docs]@_idx_fct
def TCWET(bands: dict) -> xr.DataArray:
"""
Tasseled Cap Wetness:
https://en.wikipedia.org/wiki/Tasseled_cap_transformation
https://www.indexdatabase.de/db/r-single.php?id=723
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return (
0.1509 * bands[obn.BLUE]
+ 0.1973 * bands[obn.GREEN]
+ 0.3279 * bands[obn.RED]
+ 0.3406 * bands[obn.NIR]
- 0.7112 * bands[obn.SWIR_1]
- 0.4572 * bands[obn.SWIR_2]
)
[docs]@_idx_fct
def NDRE2(bands: dict) -> xr.DataArray:
"""
Normalized Difference Red-Edge: https://www.indexdatabase.de/db/i-single.php?id=223
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.NIR], bands[obn.VRE_1])
[docs]@_idx_fct
def NDRE3(bands: dict) -> xr.DataArray:
"""
Normalized Difference Red-Edge: https://www.indexdatabase.de/db/i-single.php?id=223
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.NIR], bands[obn.VRE_2])
[docs]@_idx_fct
def GLI(bands: dict) -> xr.DataArray:
"""
Green leaf index: https://www.indexdatabase.de/db/i-single.php?id=375
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return (2 * bands[obn.GREEN] - bands[obn.RED] - bands[obn.BLUE]) / (
2 * bands[obn.GREEN] + bands[obn.RED] + bands[obn.BLUE]
)
[docs]@_idx_fct
def GNDVI(bands: dict) -> xr.DataArray:
"""
Green NDVI: https://www.indexdatabase.de/db/i-single.php?id=401
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.NIR], bands[obn.GREEN])
[docs]@_idx_fct
def RI(bands: dict) -> xr.DataArray:
"""
Normalized Difference RED/GREEN Redness Index: https://www.indexdatabase.de/db/i-single.php?id=74
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.VRE_1], +bands[obn.GREEN])
[docs]@_idx_fct
def NDGRI(bands: dict) -> xr.DataArray:
"""
Normalized Difference GREEN/RED Index: https://www.indexdatabase.de/db/i-single.php?id=390
Also known as NDGR.
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.GREEN], bands[obn.RED])
[docs]@_idx_fct
def CIG(bands: dict) -> xr.DataArray:
"""
Chlorophyll Index Green: https://www.indexdatabase.de/db/i-single.php?id=128
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return (bands[obn.NIR] / bands[obn.GREEN]) - 1
[docs]@_idx_fct
def NDMI(bands: dict) -> xr.DataArray:
"""
Normalized Difference Moisture Index: https://www.indexdatabase.de/db/i-single.php?id=56
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.NIR], bands[obn.SWIR_1])
[docs]@_idx_fct
def DSWI(bands: dict) -> xr.DataArray:
"""
Disease water stress index: https://www.indexdatabase.de/db/i-single.php?id=106
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return (bands[obn.NIR] + bands[obn.GREEN]) / (bands[obn.SWIR_1] + bands[obn.RED])
[docs]@_idx_fct
def SRSWIR(bands: dict) -> xr.DataArray:
"""
Simple Ratio SWIR_1/SWIR_2 Clay Minerals: https://www.indexdatabase.de/db/i-single.php?id=204
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return bands[obn.SWIR_1] / bands[obn.SWIR_2]
[docs]@_idx_fct
def RDI(bands: dict) -> xr.DataArray:
"""
Ratio Drought Index: https://www.indexdatabase.de/db/i-single.php?id=71
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return bands[obn.SWIR_2] / bands[obn.NARROW_NIR]
[docs]@_idx_fct
def NDWI(bands: dict) -> xr.DataArray:
"""
Normalized Difference Water Index (GREEN Version):
https://pro.arcgis.com/fr/pro-app/2.7/arcpy/image-analyst/ndwi.htm
NDWI = (Green - NIR) / (Green + NIR)
For the SWIR version, see the NDMI.
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.GREEN], bands[obn.NIR])
[docs]@_idx_fct
def BAI(bands: dict) -> xr.DataArray:
"""
Burn Area Index: https://www.harrisgeospatial.com/docs/BackgroundBurnIndices.html
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return 1.0 / ((0.1 - bands[obn.RED]) ** 2 + (0.06 - bands[obn.NIR]) ** 2)
[docs]@_idx_fct
def BAIS2(bands: dict) -> xr.DataArray:
"""
Burn Area Index for Sentinel-2:
https://www.researchgate.net/publication/323964124_BAIS2_Burned_Area_Index_for_Sentinel-2
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
# (1-((B06*B07*B8A)/B04)**0.5)*((B12-B8A)/((B12+B8A)**0.5)+1);
a = (
(bands[obn.VRE_2] * bands[obn.VRE_3] * bands[obn.NARROW_NIR]) / bands[obn.RED]
) ** 0.5
b = (bands[obn.SWIR_2] - bands[obn.NARROW_NIR]) / (
(bands[obn.SWIR_2] + bands[obn.NARROW_NIR]) ** 0.5
)
return (1 - a) * (1 + b)
[docs]@_idx_fct
def NBR(bands: dict) -> xr.DataArray:
"""
Normalized Burn Ratio: https://www.indexdatabase.de/db/i-single.php?id=53
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.NARROW_NIR], bands[obn.SWIR_2])
[docs]@_idx_fct
def MNDWI(bands: dict) -> xr.DataArray:
"""
Modified Normalised Difference Water Index : https://wiki.orfeo-toolbox.org/index.php/MNDWI
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.GREEN], bands[obn.SWIR_1])
[docs]@_idx_fct
def AWEInsh(bands: dict) -> xr.DataArray:
"""
Automated Water Extraction Index not shadow: Feyisa et al. (2014)
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return 4 * (bands[obn.GREEN] - bands[obn.SWIR_1]) - (
0.25 * bands[obn.NIR] + 2.75 * bands[obn.SWIR_2]
)
[docs]@_idx_fct
def AWEIsh(bands: dict) -> xr.DataArray:
"""
Automated Water Extraction Index shadow: Feyisa et al. (2014)
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return (
bands[obn.BLUE]
+ 2.5 * bands[obn.GREEN]
- 1.5 * (bands[obn.NIR] + bands[obn.SWIR_1])
- 0.25 * bands[obn.SWIR_2]
)
[docs]@_idx_fct
def WI(bands: dict) -> xr.DataArray:
"""
Water Index (2015): Fisher et al. (2016)
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return (
1.7204
+ 171 * bands[obn.GREEN]
+ 3 * bands[obn.RED]
- 70 * bands[obn.NIR]
- 45 * bands[obn.SWIR_1]
- 71 * bands[obn.SWIR_2]
)
[docs]@_idx_fct
def AFRI_1_6(bands: dict) -> xr.DataArray:
"""
Aerosol free vegetation index 1600: https://www.indexdatabase.de/db/i-single.php?id=393
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.NIR], 0.66 * bands[obn.SWIR_1])
[docs]@_idx_fct
def AFRI_2_1(bands: dict) -> xr.DataArray:
"""
Aerosol free vegetation index 2100: https://www.indexdatabase.de/db/i-single.php?id=395
.. WARNING::
There is an error in the formula, go see the papers to get the right one (0.56 instead of 0.5):
https://core.ac.uk/download/pdf/130673386.pdf
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.NIR], 0.5 * bands[obn.SWIR_2])
[docs]@_idx_fct
def BSI(bands: dict) -> xr.DataArray:
"""
Barren Soil Index:
Rikimaru et al., 2002. Tropical forest cover density mapping.
http://tropecol.com/pdf/open/PDF_43_1/43104.pdf
BSI = ((RED+SWIR) – (NIR+BLUE)) / ((RED+SWIR) + (NIR+BLUE))
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(
bands[obn.RED] + bands[obn.SWIR_1], bands[obn.NIR] + bands[obn.BLUE]
)
# WorldView index (without the ones with SWIR)
# https://resources.maxar.com/optical-imagery/multispectral-reference-guide
[docs]@_idx_fct
def WV_WI(bands: dict) -> xr.DataArray:
"""
WorldView-Water (WV-WI)
Useful for detecting standing, flowing water, or shadow in VNIR imagery
:code:`WV_WI = ((B8-B1)/(B8+B1))`
https://resources.maxar.com/optical-imagery/multispectral-reference-guide
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.WV], bands[obn.CA])
[docs]@_idx_fct
def WV_VI(bands: dict) -> xr.DataArray:
"""
WorldView – Vegetation (WV-VI)
Useful for detecting vegetation and assessing vegetation health
:code:`WV_VI = ((B8-B5)/(B8+B5))`
https://resources.maxar.com/optical-imagery/multispectral-reference-guide
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.WV], bands[obn.RED])
[docs]@_idx_fct
def WV_SI(bands: dict) -> xr.DataArray:
"""
WorldView – Soil (WV-SI)
Useful for detecting and differentiating exposed soil
:code:`WV_SI = ((B4-B3)/(B4+B3))`
https://resources.maxar.com/optical-imagery/multispectral-reference-guide
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.YELLOW], bands[obn.GREEN])
[docs]@_idx_fct
def WV_BI(bands: dict) -> xr.DataArray:
"""
WorldView – Built-up (WV-BI)
Useful for detecting impervious surfaces such as buildings and roads
:code:`WV_BI = ((B6-B1)/(B6+B1))`
https://resources.maxar.com/optical-imagery/multispectral-reference-guide
Args:
bands (dict): Bands as {band_name: xr.DataArray}
Returns:
xr.DataArray: Computed index
"""
return _norm_diff(bands[obn.VRE_1], bands[obn.CA])
[docs]def get_all_index_names() -> list:
"""
Get all index names contained in this file
.. code-block:: python
>>> from eoreader.bands import index
>>> index.get_all_index_names()
['AFRI_1_6', 'AFRI_2_1', 'AWEInsh', 'AWEIsh', 'BAI', ..., 'WI']
Returns:
list: Index names
"""
return [idx_fct.__name__ for idx_fct in get_all_index()]
[docs]def get_all_index() -> list:
"""
Get all index functions contained in this file
.. code-block:: python
>>> from eoreader.bands import index
>>> index.get_all_index()
[<function AFRI_1_6 at 0x00000118FFFB51E0>, ..., <function WI at 0x00000118FFFB5158>]
Returns:
list: Index functions
"""
idx = []
functions = inspect.getmembers(sys.modules[__name__], predicate=inspect.isfunction)
for (name, fct) in functions:
# Do not gather this fct nor da.true_divide
if name[0].isupper():
idx.append(fct)
return idx
[docs]def get_needed_bands(index: Callable) -> list:
"""
Gather all the needed bands for the specified index function
.. code-block:: python
>>> index.get_needed_bands(NDVI)
[<OpticalBandNames.NIR: 'NIR'>, <OpticalBandNames.RED: 'RED'>]
Returns:
list: Needed bands for the index function
"""
# Get source code from this fct
code = inspect.getsource(index)
# Parse band's signature
b_regex = r"obn\.\w+"
return [getattr(obn, b.split(".")[-1]) for b in re.findall(b_regex, code)]
[docs]def get_all_needed_bands() -> dict:
"""
Gather all the needed bands for all index functions
.. code-block:: python
>>> index.get_all_needed_bands()
{
<function AFRI_1_6 at 0x00000261F6FF36A8>: [<OpticalBandNames.NIR: 'NIR'>, <OpticalBandNames.SWIR_2: 'SWIR_2'>],
...
<function WI at 0x00000261F6FF3620>: [<OpticalBandNames.NIR: 'NIR'>, <OpticalBandNames.SWIR_1: 'SWIR_1'>]
}
>>> # Or written in a more readable fashion:
>>> {idx.__name__: [band.value for band in bands] for idx, bands in index.get_all_needed_bands().items()}
{
'AFRI_1_6': ['NIR', 'SWIR_2'],
...,
'WI': ['NIR', 'SWIR_1']
}
Returns:
dict: Needed bands for all index functions
"""
needed_bands = {}
# Get all function from this file
functions = inspect.getmembers(sys.modules[__name__], predicate=inspect.isfunction)
for (name, function) in functions:
# Do not gather this fct nor da.true_divide
if name[0].isupper():
needed_bands[function] = get_needed_bands(function)
return needed_bands
NEEDED_BANDS = get_all_needed_bands()