{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# **Spectral Indices**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[](https://colab.research.google.com/github/kenoz/SITS_utils/blob/main/docs/source/tutorials/colab_sits_ex03.ipynb)" ] }, { "cell_type": "markdown", "metadata": { "id": "K2FIkbDYrq9l" }, "source": [ "---\n", "\n", "Using **Banc d'Arguin** as an example, we aim to **calculate one or more spectral indices**. To avoid reinventing the wheel, `sits` relies on the [spyndex](https://github.com/awesome-spectral-indices/spyndex) package. This package leverages spectral indices from the [Awesome Spectral Indices](https://github.com/awesome-spectral-indices/awesome-spectral-indices) list and provides an expression evaluation method compatible with Python object classes that support overloaded operators (e.g., `numpy.ndarray`, `pandas.Series`, `xarray.DataArray`).\n", "\n", "---\n", "\n", "## 1. Installation of SITS package and its depedencies\n", "\n", "First, install `sits` package with [pip](https://pypi.org/project/SITS/). We also need some other packages for displaying data." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "id": "xoL6NstiVcp9" }, "outputs": [], "source": [ "# SITS package\n", "!pip install -q --upgrade sits\n", "\n", "# other packages\n", "!pip install -q \"dask[dataframe]\"" ] }, { "cell_type": "markdown", "metadata": { "id": "OYAb0l91zqj4" }, "source": [ "Now we can import `sits` and some other libraries." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "executionInfo": { "elapsed": 8932, "status": "ok", "timestamp": 1737993800099, "user": { "displayName": "Kenji Ose", "userId": "04036279333112259084" }, "user_tz": -60 }, "id": "neqafgGHWIkj" }, "outputs": [], "source": [ "import os\n", "# sits lib\n", "from sits import sits, export\n", "# geospatial libs\n", "import geopandas as gpd\n", "import pandas as pd\n", "import spyndex\n", "# date format\n", "from datetime import datetime\n", "# ignore warnings messages\n", "import warnings\n", "warnings.filterwarnings('ignore')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Exploring spectral indices\n", "\n", "Scientific literature offers a wide array of **spectral indices** designed to reveal specific characteristics of objects or phenomena on the Earth's surface. The [**Awesome Spectral Indices**](https://awesome-ee-spectral-indices.readthedocs.io/en/latest/) (ASI) project, initiated in late 2021, catalogs these indices—both well-established and emerging—for **optical and radar** sensors. Each index is linked to a **standardized formula**, enabling consistent computation across imagery from various satellites. Furthermore, every index is accompanied by a scientific reference, ensuring transparency and traceability.\n", "\n", "The [**Spyndex**](https://spyndex.readthedocs.io/en/latest/) package provides seamless access to this library. It offers robust tools to **retrieve metadata** and efficiently **compute spectral indices**, streamlining their use in geospatial analyses.\n", "\n", "Reference:\n", "- Montero, D., Aybar, C., Mahecha, M.D. et al. A standardized catalogue of spectral indices to advance the use of remote sensing in Earth system research. *Sci Data* **10**, 197 (2023). https://doi.org/10.1038/s41597-023-02096-0\n", "\n", "### 2.1. Which indices? Which Spectral bands?\n", "\n", "Before loading satellite time series, we need to know which spectral indices we would like to calculate so that we can request the relevant spectral bands. First we list with Spyndex all the available indices." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "number of indices. 245\n", "---\n", "available indices:\n", "['AFRI1600', 'AFRI2100', 'ANDWI', 'ARI', 'ARI2', 'ARVI', 'ATSAVI', 'AVI', 'AWEInsh', 'AWEIsh', 'BAI', 'BAIM', 'BAIS2', 'BCC', 'BI', 'BITM', 'BIXS', 'BLFEI', 'BNDVI', 'BRBA', 'BWDRVI', 'BaI', 'CCI', 'CIG', 'CIRE', 'CSI', 'CSIT', 'CVI', 'DBI', 'DBSI', 'DPDD', 'DSI', 'DSWI1', 'DSWI2', 'DSWI3', 'DSWI4', 'DSWI5', 'DVI', 'DVIplus', 'DpRVIHH', 'DpRVIVV', 'EBBI', 'EBI', 'EMBI', 'ENDVI', 'EVI', 'EVI2', 'EVIv', 'ExG', 'ExGR', 'ExR', 'FAI', 'FCVI', 'GARI', 'GBNDVI', 'GCC', 'GDVI', 'GEMI', 'GLI', 'GM1', 'GM2', 'GNDVI', 'GOSAVI', 'GRNDVI', 'GRVI', 'GSAVI', 'GVMI', 'IAVI', 'IBI', 'IKAW', 'IPVI', 'IRECI', 'LSWI', 'MBI', 'MBWI', 'MCARI', 'MCARI1', 'MCARI2', 'MCARI705', 'MCARIOSAVI', 'MCARIOSAVI705', 'MGRVI', 'MIRBI', 'MLSWI26', 'MLSWI27', 'MNDVI', 'MNDWI', 'MNLI', 'MRBVI', 'MSAVI', 'MSI', 'MSR', 'MSR705', 'MTCI', 'MTVI1', 'MTVI2', 'MuWIR', 'NBAI', 'NBLI', 'NBLIOLI', 'NBR', 'NBR2', 'NBRSWIR', 'NBRT1', 'NBRT2', 'NBRT3', 'NBRplus', 'NBSIMS', 'NBUI', 'ND705', 'NDBI', 'NDBaI', 'NDCI', 'NDDI', 'NDGI', 'NDGlaI', 'NDII', 'NDISIb', 'NDISIg', 'NDISImndwi', 'NDISIndwi', 'NDISIr', 'NDMI', 'NDPI', 'NDPolI', 'NDPonI', 'NDREI', 'NDSI', 'NDSII', 'NDSIWV', 'NDSInw', 'NDSWIR', 'NDSaII', 'NDSoI', 'NDTI', 'NDVI', 'NDVI705', 'NDVIMNDWI', 'NDVIT', 'NDWI', 'NDWIns', 'NDYI', 'NGRDI', 'NHFD', 'NIRv', 'NIRvH2', 'NIRvP', 'NLI', 'NMDI', 'NRFIg', 'NRFIr', 'NSDS', 'NSDSI1', 'NSDSI2', 'NSDSI3', 'NSTv1', 'NSTv2', 'NWI', 'NormG', 'NormNIR', 'NormR', 'OCVI', 'OSAVI', 'OSI', 'PI', 'PISI', 'PSRI', 'QpRVI', 'RCC', 'RDVI', 'REDSI', 'RENDVI', 'RFDI', 'RGBVI', 'RGRI', 'RI', 'RI4XS', 'RNDVI', 'RVI', 'S2REP', 'S2WI', 'S3', 'SARVI', 'SAVI', 'SAVI2', 'SAVIT', 'SEVI', 'SI', 'SIPI', 'SLAVI', 'SR', 'SR2', 'SR3', 'SR555', 'SR705', 'SWI', 'SWM', 'SeLI', 'TCARI', 'TCARIOSAVI', 'TCARIOSAVI705', 'TCI', 'TDVI', 'TGI', 'TRRVI', 'TSAVI', 'TTVI', 'TVI', 'TWI', 'TriVI', 'UI', 'VARI', 'VARI700', 'VDDPI', 'VHVVD', 'VHVVP', 'VHVVR', 'VI6T', 'VI700', 'VIBI', 'VIG', 'VVVHD', 'VVVHR', 'VVVHS', 'VgNIRBI', 'VrNIRBI', 'WDRVI', 'WDVI', 'WI1', 'WI2', 'WI2015', 'WRI', 'bNIRv', 'kEVI', 'kIPVI', 'kNDVI', 'kRVI', 'kVARI', 'mND705', 'mSR705', 'sNIRvLSWI', 'sNIRvNDPI', 'sNIRvNDVILSWIP', 'sNIRvNDVILSWIS', 'sNIRvSWIR']\n" ] } ], "source": [ "idx_list = spyndex.indices\n", "\n", "print(f\"number of indices. {len(idx_list)}\")\n", "print(F\"---\\navailable indices:\")\n", "print(idx_list)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To get details about an index, we can use the following spyndex attributes:\n", "- `spyndex.indices.[*index short name*]`: index's metadata\n", "- `spyndex.indices.[*index short name*].short_name`: Short name of the Spectral Index.\n", "- `spyndex.indices.[*index short name*].long_name`: Long name of the Spectral Index.\n", "- `spyndex.indices.[*index short name*].application_domain`: Required bands and parameters for the Spectral Index computation.\n", "- `spyndex.indices.[*index short name*].bands`: Required bands and parameters for the Spectral Index computation.\n", "- `spyndex.indices.[*index short name*].formula`: Formula (as expression) of the Spectral Index.\n", "- `spyndex.indices.[*index short name*].reference`: URL to the reference/DOI of the Spectral Index.\n", "- `spyndex.indices.[*index short name*].platforms`: Platforms with the required bands for the Spectral Index computation.\n", "\n", "For example, we explore here the quite unknown **NDVI**:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "metadata:\n", "NDVI: Normalized Difference Vegetation Index\n", " * Application Domain: vegetation\n", " * Bands/Parameters: ['N', 'R']\n", " * Formula: (N-R)/(N+R)\n", " * Reference: https://ntrs.nasa.gov/citations/19740022614\n", " \n", "---\n", " short name: NDVI\n", "---\n", " long name:Normalized Difference Vegetation Index\n", "---\n", " application domain: vegetation\n", "---\n", " bands: ['N', 'R']\n", "---\n", " formula: (N-R)/(N+R)\n", "---\n", " reference: https://ntrs.nasa.gov/citations/19740022614\n", "---\n", " platform: ['Sentinel-2', 'Landsat-OLI', 'Landsat-TM', 'Landsat-ETM+', 'MODIS', 'Planet-Fusion']\n" ] } ], "source": [ "print(f\"metadata:\\n{spyndex.indices.NDVI}\")\n", "\n", "print(f\"---\\n short name: {spyndex.indices.NDVI.short_name}\")\n", "print(f\"---\\n long name:{spyndex.indices.NDVI.long_name}\")\n", "print(f\"---\\n application domain: {spyndex.indices.NDVI.application_domain}\")\n", "print(f\"---\\n bands: {spyndex.indices.NDVI.bands}\")\n", "print(f\"---\\n formula: {spyndex.indices.NDVI.formula}\")\n", "print(f\"---\\n reference: {spyndex.indices.NDVI.reference}\")\n", "print(f\"---\\n platform: {spyndex.indices.NDVI.platforms}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will compute three spectral indices:\n", "\n", "- NDVI, \n", "- NDWI, and \n", "- BI\n", "\n", "All of which include accessible metadata, just like NDVI." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "NDWI metadata:\n", "NDWI: Normalized Difference Water Index\n", " * Application Domain: water\n", " * Bands/Parameters: ['G', 'N']\n", " * Formula: (G-N)/(G+N)\n", " * Reference: https://doi.org/10.1080/01431169608948714\n", " \n", "EVI metadata:\n", "BI: Bare Soil Index\n", " * Application Domain: soil\n", " * Bands/Parameters: ['S1', 'R', 'N', 'B']\n", " * Formula: ((S1+R)-(N+B))/((S1+R)+(N+B))\n", " * Reference: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.465.8749&rep=rep1&type=pdf\n", " \n" ] } ], "source": [ "print(f\"NDWI metadata:\\n{spyndex.indices.NDWI}\")\n", "print(f\"EVI metadata:\\n{spyndex.indices.BI}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2. Selection of spectral bands\n", "\n", "We search for all necessary spectral bands to calculate the indices. As the band names are standardized, we need to map them to the Sentinel-2 nomenclature." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "We need to request the following bands:\n", "{'R', 'N', 'S1', 'G', 'B'}\n" ] } ], "source": [ "common_bands = set()\n", "common_bands.update(spyndex.indices.NDVI.bands)\n", "common_bands.update(spyndex.indices.NDWI.bands)\n", "common_bands.update(spyndex.indices.BI.bands)\n", "\n", "print(\"We need to request the following bands:\")\n", "print(common_bands)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Spyndex allows us to explore detailed information about spectral bands, including their corresponding names across different naming standards." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "---\n", "Short name of the Band: R\n", " - Description/Name of the Band: Red\n", " - Common name of the Band according to the EO Ext. Specification for STAC.: red\n", " - Minimum wavelength of the spectral range of the band (nm): 620\n", " - Maximum wavelength of the spectral range of the band (nm): 690\n", " - Sentinel-2 band name: B4\n", "---\n", "Short name of the Band: N\n", " - Description/Name of the Band: Near-Infrared (NIR)\n", " - Common name of the Band according to the EO Ext. Specification for STAC.: nir\n", " - Minimum wavelength of the spectral range of the band (nm): 760\n", " - Maximum wavelength of the spectral range of the band (nm): 900\n", " - Sentinel-2 band name: B8\n", "---\n", "Short name of the Band: S1\n", " - Description/Name of the Band: Short-wave Infrared (SWIR) 1\n", " - Common name of the Band according to the EO Ext. Specification for STAC.: swir16\n", " - Minimum wavelength of the spectral range of the band (nm): 1550\n", " - Maximum wavelength of the spectral range of the band (nm): 1750\n", " - Sentinel-2 band name: B11\n", "---\n", "Short name of the Band: G\n", " - Description/Name of the Band: Green\n", " - Common name of the Band according to the EO Ext. Specification for STAC.: green\n", " - Minimum wavelength of the spectral range of the band (nm): 510\n", " - Maximum wavelength of the spectral range of the band (nm): 600\n", " - Sentinel-2 band name: B3\n", "---\n", "Short name of the Band: B\n", " - Description/Name of the Band: Blue\n", " - Common name of the Band according to the EO Ext. Specification for STAC.: blue\n", " - Minimum wavelength of the spectral range of the band (nm): 450\n", " - Maximum wavelength of the spectral range of the band (nm): 530\n", " - Sentinel-2 band name: B2\n" ] } ], "source": [ "for band in common_bands:\n", " print(f\"---\\nShort name of the Band: {spyndex.bands[band].short_name}\")\n", " print(f\" - Description/Name of the Band: {spyndex.bands[band].long_name}\")\n", " print(f\" - Common name of the Band according to the EO Ext. Specification for STAC.: {spyndex.bands[band].common_name}\")\n", " print(f\" - Minimum wavelength of the spectral range of the band (nm): {spyndex.bands[band].min_wavelength}\")\n", " print(f\" - Maximum wavelength of the spectral range of the band (nm): {spyndex.bands[band].max_wavelength}\")\n", " print(f\" - Sentinel-2 band name: {spyndex.bands[band].sentinel2a.band}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Unfortunately, the Sentinel-2 band names don’t align perfectly with those used in the STAC catalog of the Microsoft Planetary Computer, so we need to map them manually." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [] }, "outputs": [], "source": [ "band_mapping = {'B': 'B02', 'G': 'B03', 'R': 'B04', 'N': 'B08', 'S1': 'B11'}" ] }, { "cell_type": "markdown", "metadata": { "id": "o4n__qxr0Kqi" }, "source": [ "## 3. Loading the sentinel-2 time series\n", "\n", "### 3.1. Data loading\n", "\n", "The geojson vector file describing the position of the sandbank is stored in the [Github repository](https://github.com/kenoz/SITS_utils). We download it into our current workspace. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "executionInfo": { "elapsed": 218, "status": "ok", "timestamp": 1737993800315, "user": { "displayName": "Kenji Ose", "userId": "04036279333112259084" }, "user_tz": -60 }, "id": "soIORtc8bGDd" }, "outputs": [], "source": [ "!mkdir -p test_data\n", "![ ! -f test_data/banc_arguin.geojson ] && wget https://raw.githubusercontent.com/kenoz/SITS_utils/refs/heads/main/sits/data/banc_arguin.geojson -P test_data" ] }, { "cell_type": "markdown", "metadata": { "id": "n0jbqcwQ1ppL" }, "source": [ "We load the vector file, named `banc_arguin.geojson`, as a geoDataFrame object with the `sits` method: `sits.Vec2gdf()`." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 81 }, "executionInfo": { "elapsed": 697, "status": "ok", "timestamp": 1737993801010, "user": { "displayName": "Kenji Ose", "userId": "04036279333112259084" }, "user_tz": -60 }, "id": "7MWDWpr7nxB3", "outputId": "cd939689-bad7-4b42-d469-a1a34f6a69ef" }, "outputs": [], "source": [ "data_dir = 'test_data'\n", "\n", "# Load vector\n", "v_arguin = sits.Vec2gdf(os.path.join(data_dir, 'banc_arguin.geojson'))\n", "\n", "# Define polygon bounding box\n", "v_arguin.set_bbox('gdf')\n", "\n", "# CRS management\n", "bbox_4326 = list(v_arguin.bbox.iloc[0]['geometry'].bounds)\n", "bbox_3035 = list(v_arguin.bbox.to_crs(3035).iloc[0]['geometry'].bounds)" ] }, { "cell_type": "markdown", "metadata": { "id": "ipDQq-9jqWjo" }, "source": [ "### 3.2. Loading and preprocessing of a Satellite Image Time-Series (SITS)\n", "\n", "In this example, we have only one area (one polygon) to process. We use the class `sits.StacAttack` to request and preprocess the data. The request consists in retrieving Sentinel-2 images (level 2A) acquired from January 1, 2016 to January 1, 2025 with cloud cover less than 10%. Then we build a 4 bands geo-datacube ('B03', 'B04', 'B08' and 'SCL') in EPSG:3035 with a 20m spatial resolution." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "executionInfo": { "elapsed": 7130, "status": "ok", "timestamp": 1737993825234, "user": { "displayName": "Kenji Ose", "userId": "04036279333112259084" }, "user_tz": -60 }, "id": "N8bn6ail-46F" }, "outputs": [ { "data": { "text/html": [ "
<xarray.Dataset> Size: 277MB\n",
"Dimensions: (y: 517, x: 413, time: 108)\n",
"Coordinates:\n",
" * y (y) float64 4kB 2.459e+06 2.459e+06 ... 2.448e+06 2.448e+06\n",
" * x (x) float64 3kB 3.426e+06 3.426e+06 ... 3.435e+06 3.435e+06\n",
" * time (time) datetime64[ns] 864B 2016-03-18T11:11:02.030000 ... 20...\n",
" spatial_ref int64 8B 0\n",
"Data variables:\n",
" B02 (time, y, x) uint16 46MB dask.array<chunksize=(1, 517, 413), meta=np.ndarray>\n",
" B03 (time, y, x) uint16 46MB dask.array<chunksize=(1, 517, 413), meta=np.ndarray>\n",
" B04 (time, y, x) uint16 46MB dask.array<chunksize=(1, 517, 413), meta=np.ndarray>\n",
" B08 (time, y, x) uint16 46MB dask.array<chunksize=(1, 517, 413), meta=np.ndarray>\n",
" B11 (time, y, x) uint16 46MB dask.array<chunksize=(1, 517, 413), meta=np.ndarray>\n",
" SCL (time, y, x) uint16 46MB dask.array<chunksize=(1, 517, 413), meta=np.ndarray><xarray.Dataset> Size: 184MB\n",
"Dimensions: (y: 517, x: 413, time: 108)\n",
"Coordinates:\n",
" * y (y) float64 4kB 2.459e+06 2.459e+06 ... 2.448e+06 2.448e+06\n",
" * x (x) float64 3kB 3.426e+06 3.426e+06 ... 3.435e+06 3.435e+06\n",
" * time (time) datetime64[ns] 864B 2016-03-18T11:11:02.030000 ... 20...\n",
" spatial_ref int64 8B 0\n",
"Data variables:\n",
" NDVI (time, y, x) float64 184MB dask.array<chunksize=(1, 517, 413), meta=np.ndarray><xarray.Dataset> Size: 553MB\n",
"Dimensions: (y: 517, x: 413, time: 108)\n",
"Coordinates:\n",
" * y (y) float64 4kB 2.459e+06 2.459e+06 ... 2.448e+06 2.448e+06\n",
" * x (x) float64 3kB 3.426e+06 3.426e+06 ... 3.435e+06 3.435e+06\n",
" * time (time) datetime64[ns] 864B 2016-03-18T11:11:02.030000 ... 20...\n",
"Data variables:\n",
" NDVI (time, y, x) float64 184MB dask.array<chunksize=(1, 517, 413), meta=np.ndarray>\n",
" NDWI (time, y, x) float64 184MB dask.array<chunksize=(1, 517, 413), meta=np.ndarray>\n",
" BI (time, y, x) float64 184MB dask.array<chunksize=(1, 517, 413), meta=np.ndarray>\n",
" spatial_ref int64 8B 0