{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DYndcZst_kdr"
   },
   "outputs": [],
   "source": [
    "# Copyright (c) TorchGeo Contributors. All rights reserved.\n",
    "# Licensed under the MIT License."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZKIkyiLScf9P"
   },
   "source": [
    "# Spectral Indices\n",
    "\n",
    "_Written by: Isaac A. Corley_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PevsPoE4cY0j"
   },
   "source": [
    "In this tutorial, we demonstrate how to use TorchGeo's functions and transforms to compute popular indices used in remote sensing and provide examples of how to utilize them for analyzing raw imagery or simply for visualization purposes. Some common indices and their formulas can be found at the following links:\n",
    "\n",
    "- [Index Database](https://www.indexdatabase.de/db/i.php)\n",
    "- [Awesome Spectral Indices](https://github.com/awesome-spectral-indices/awesome-spectral-indices)\n",
    "\n",
    "It's recommended to run this notebook on Google Colab if you don't have your own GPU. Click the \"Open in Colab\" button above to get started."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fsOYw-p2ccka"
   },
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VqdMMzvacOF8"
   },
   "source": [
    "Install TorchGeo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "wOwsb8KT_uXR",
    "outputId": "7a7ca2ff-a9a5-444f-99c8-61954b82c797"
   },
   "outputs": [],
   "source": [
    "%pip install torchgeo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "u2f5_f4X_-vV"
   },
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cvPMr76K_9uk"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import tempfile\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from torchgeo.datasets import EuroSAT100\n",
    "from torchgeo.transforms import AppendNDBI, AppendNDVI, AppendNDWI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Uh2IpthodK1R"
   },
   "source": [
    "## Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vLW4KhtEwx-e"
   },
   "source": [
    "We will use the [EuroSAT](https://torchgeo.readthedocs.io/en/stable/api/datasets.html#eurosat) dataset throughout this tutorial. Specifically, a subset containing only 100 images. EuroSAT contains 13-channel multispectral imagery captured by the Sentinel-2 satellite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "_seqhOz-Cw9c",
    "outputId": "99187c2c-9f4f-4ca5-ec10-6635feaaf064"
   },
   "outputs": [],
   "source": [
    "root = os.path.join(tempfile.gettempdir(), 'eurosat100')\n",
    "ds = EuroSAT100(root, download=True)\n",
    "sample = ds[21]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DLdI6L4A08vu"
   },
   "source": [
    "## True Color (RGB) Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nkb51TZpPylY"
   },
   "source": [
    "We can plot a true color image consisting of the first 3 channels (RGB) to visualize the sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "tymUtqGAQCQl",
    "outputId": "a5df2986-cc56-410a-ff73-4fd7caec2b80"
   },
   "outputs": [],
   "source": [
    "ds.plot(sample)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jWmSuZfI1gwN"
   },
   "source": [
    "## Normalized Difference Vegetation Index (NDVI)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yGTZzk6_QH1a"
   },
   "source": [
    "Below we use TorchGeo's `indices.AppendNDVI` to compute the [Normalized Difference Vegetation Index (NDVI)](https://gisgeography.com/ndvi-normalized-difference-vegetation-index/). NDVI is useful for measuring the presence of vegetation and vegetation health. It can be calculated using the Near Infrared (NIR) and Red bands using the formula below, resulting in a value between [-1, 1] where low NDVI values represents no or unhealthy vegetation and high NDVI values represents healthy vegetation. Here we use a diverging red, yellow, green colormap representing -1, 0, and 1, respectively.\n",
    "\n",
    "$$\\text{NDVI} = \\frac{\\text{NIR} - \\text{R}}{\\text{NIR} + \\text{R}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "e9Aob95YQCQn",
    "outputId": "0f12ca7e-d8ed-4327-8a41-0674d85e1df1"
   },
   "outputs": [],
   "source": [
    "# NDVI is appended to channel dimension (dim=0)\n",
    "index = AppendNDVI(index_nir=7, index_red=3)\n",
    "image = sample['image']\n",
    "image = index(image)[0]\n",
    "\n",
    "# Normalize from [-1, 1] -> [0, 1] for visualization\n",
    "image[-1] = (image[-1] + 1) / 2\n",
    "\n",
    "plt.imshow(image[-1], cmap='RdYlGn')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "igNq-9m91nDt"
   },
   "source": [
    "## Normalized Difference Water Index (NDWI)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fDFQm4vOZEUH"
   },
   "source": [
    "Below we use TorchGeo's `indices.AppendNDWI` to compute the [Normalized Difference Water Index (NDWI)](https://custom-scripts.sentinel-hub.com/custom-scripts/sentinel-2/ndwi/). NDWI is useful for measuring the presence of water content in water bodies. It can be calculated using the Green and Near Infrared (NIR) bands using the formula below, resulting in a value between [-1, 1] where low NDWI values represents no water and high NDWI values represents water bodies. Here we use a diverging brown, white, blue-green colormap representing -1, 0, and 1, respectively.\n",
    "\n",
    "$$\\text{NDWI} = \\frac{\\text{G} - \\text{NIR}}{\\text{G} + \\text{NIR}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "H8CPnPD9QCQp",
    "outputId": "81b48a45-ab72-496b-cc4c-283b8a239396"
   },
   "outputs": [],
   "source": [
    "# NDWI is appended to channel dimension (dim=0)\n",
    "index = AppendNDWI(index_green=2, index_nir=7)\n",
    "image = index(image)[0]\n",
    "\n",
    "# Normalize from [-1, 1] -> [0, 1] for visualization\n",
    "image[-1] = (image[-1] + 1) / 2\n",
    "\n",
    "plt.imshow(image[-1], cmap='BrBG')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oKEyz9TP2OK_"
   },
   "source": [
    "## Normalized Difference Built-up Index (NDBI)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ap3adJ06bXul"
   },
   "source": [
    "Below we use TorchGeo's `indices.AppendNDBI` to compute the [Normalized Difference Built-up Index (NDBI)](https://www.linkedin.com/pulse/ndvi-ndbi-ndwi-calculation-using-landsat-7-8-tek-bahadur-kshetri/). NDBI is useful for measuring the presence of urban buildings. It can be calculated using the Short-wave Infrared (SWIR)  and Near Infrared (NIR) bands using the formula below, resulting in a value between [-1, 1] where low NDBI values represents no urban land and high NDBI values represents urban land. Here we use a terrain colormap with blue, green-yellow, and brown representing -1, 0, and 1, respectively.\n",
    "\n",
    "$$\\text{NDBI} = \\frac{\\text{SWIR} - \\text{NIR}}{\\text{SWIR} + \\text{NIR}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "F607CaDoQCQq",
    "outputId": "f18b2337-05e1-48d7-bc50-d13c536e0177"
   },
   "outputs": [],
   "source": [
    "# NDBI is appended to channel dimension (dim=0)\n",
    "index = AppendNDBI(index_swir=11, index_nir=7)\n",
    "image = index(image)[0]\n",
    "\n",
    "# Normalize from [-1, 1] -> [0, 1] for visualization\n",
    "image[-1] = (image[-1] + 1) / 2\n",
    "\n",
    "plt.imshow(image[-1], cmap='terrain')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  }
 ],
 "metadata": {
  "execution": {
   "timeout": 1200
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}