jupyter/R/stac/stac-introduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/brazil-data-cube/code-gallery/master/img/logo-bdc.png\" align=\"right\" width=\"64\"/>\n",
    "\n",
    "# <span style=\"color:#336699\">Introduction to the SpatioTemporal Asset Catalog (STAC) in R language</span>\n",
    "<hr style=\"border:2px solid #0077b9;\">\n",
    "\n",
    "<div style=\"text-align: left;\">\n",
    "    <a href=\"https://nbviewer.jupyter.org/github/brazil-data-cube/code-gallery/blob/master/jupyter/Python/stac/stac-introduction.ipynb\"><img src=\"https://raw.githubusercontent.com/jupyter/design/master/logos/Badges/nbviewer_badge.svg\" align=\"center\"/></a>\n",
    "</div>\n",
    "\n",
    "<br/>\n",
    "\n",
    "<div style=\"text-align: center;font-size: 90%;\">\n",
    "    Felipe Carvalho de Souza<sup><a href=\"https://orcid.org/0000-0002-5826-1700\"><i class=\"fab fa-lg fa-orcid\" style=\"color: #a6ce39\"></i></a></sup>, Felipe Menino Carlos<sup><a href=\"https://orcid.org/0000-0002-3334-4315\"><i class=\"fab fa-lg fa-orcid\" style=\"color: #a6ce39\"></i></a></sup>, Rennan Marujo<sup><a href=\"https://orcid.org/0000-0002-0082-9498\"><i class=\"fab fa-lg fa-orcid\" style=\"color: #a6ce39\"></i></a></sup>, Gilberto R. Queiroz<sup><a href=\"https://orcid.org/0000-0001-7534-0219\"><i class=\"fab fa-lg fa-orcid\" style=\"color: #a6ce39\"></i></a></sup>\n",
    "    <br/><br/>\n",
    "    Earth Observation and Geoinformatics Division, National Institute for Space Research (INPE)\n",
    "    <br/>\n",
    "    Avenida dos Astronautas, 1758, Jardim da Granja, São José dos Campos, SP 12227-010, Brazil\n",
    "    <br/><br/>\n",
    "    Contact: <a href=\"mailto:brazildatacube@inpe.br\">brazildatacube@inpe.br</a>\n",
    "    <br/><br/>\n",
    "    Last Update: February 13, 2023\n",
    "</div>\n",
    "\n",
    "<br/>\n",
    "\n",
    "<div style=\"text-align: justify;  margin-left: 25%; margin-right: 25%;\">\n",
    "<b>Abstract.</b> This Jupyter Notebook overviews how to use the STAC service to discover and access the data products from the <em>Brazil Data Cube</em> using <em>rstac</em> package.\n",
    "</div>\n",
    "\n",
    "<br/>\n",
    "<div style=\"text-align: justify;  margin-left: 25%; margin-right: 25%;font-size: 75%; border-style: solid; border-color: #0077b9; border-width: 1px; padding: 5px;\">\n",
    "    <b>This Jupyter Notebook is a supplement to the following paper:</b>\n",
    "    <div style=\"margin-left: 10px; margin-right: 10px\">\n",
    "    Zaglia, M.; Vinhas, L.; Queiroz, G. R.; Simões, R. <a href=\"http://urlib.net/rep/8JMKD3MGPDW34R/3UFEFD8\" target=\"_blank\">Catalogação de Metadados do Cubo de Dados do Brasil com o SpatioTemporal Asset Catalog</a>. In: Proceedings XX GEOINFO, November 11-13, 2019, São José dos Campos, SP, Brazil. p 280-285.\n",
    "    </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/brazil-data-cube/code-gallery/master/img/stac/stac.png?raw=true\" align=\"right\" width=\"66\"/>\n",
    "\n",
    "\n",
    "## <span style=\"color:#336699\">Introduction</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [**S**patio**T**emporal **A**sset **C**atalog (STAC)](https://stacspec.org/) is a specification created through several organizations' collaboration to increase satellite image search interoperability.\n",
    "\n",
    "The diagram depicted in the picture contains the most important concepts behind the STAC data model:\n",
    "\n",
    "<center>\n",
    "<img src=\"https://raw.githubusercontent.com/brazil-data-cube/code-gallery/master/img/stac/stac-concept.png\" width=\"480\" />\n",
    "<br/>\n",
    "<b>Figure 1</b> - STAC model.\n",
    "</center>\n",
    "\n",
    "The description of the concepts below are adapted from the [STAC Specification](https://github.com/radiantearth/stac-spec):\n",
    "\n",
    "- **Item**: a `STAC Item` is the atomic unit of metadata in STAC, providing links to the actual `assets` (including thumbnails) they represent. It is a `GeoJSON Feature` with additional fields for things like time, links to related entities, and mainly to the assets. According to the specification, this is the atomic unit that describes the data to be discovered in a `STAC Catalog` or `Collection`.\n",
    "\n",
    "- **Asset**: a `spatiotemporal asset` is any file representing information about the earth captured in a certain space and time.\n",
    "\n",
    "\n",
    "- **Catalog**: provides a structure to link various `STAC Items` together or even to other `STAC Catalogs` or `Collections`.\n",
    "\n",
    "\n",
    "- **Collection:** is a specialization of the `Catalog` that allows additional information about a spatio-temporal collection of data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:#336699\">Clients</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">\n",
    "\n",
    "\n",
    "The facilities provided by the STAC service can be used in any programming language that supports network communication via HTTP requests. This means that all modern languages can be used to access the operations of the service.\n",
    "\n",
    "So, in addition to the specification and implementation of the STAC service reference, the BDC also provides clients with different programming languages. These clients offer facilities that allow easy use of the STAC service in the programming languages they are implemented.\n",
    "\n",
    "- [stac.py - Python client](https://github.com/brazil-data-cube/stac.py);\n",
    "- [rstac - R Client](https://github.com/brazil-data-cube/rstac).\n",
    "\n",
    "This Jupyter Notebook will present how STAC can be used in R through the `rstac` client."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:#336699\">First step</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To run the examples in this Jupyter Notebook, you need to install the [rstac](https://github.com/brazil-data-cube/rstac) package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "install.packages(\"rstac\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are running this notebook on your local machine, consider installing the packages listed below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Remove '#' to install packages\n",
    "# install.packages(c(\"magrittr\", \"tibble\", \"dplyr\", \"raster\", \"tmap\"), dependencies = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#system(\"apt update && apt install ca-certificates\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load the `rstac` and `terra` packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(magrittr) # Package to use pipe operator %>%\n",
    "library(rstac)    # package rstac\n",
    "library(terra)    # package to manipulate rasters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-08-23T19:35:33.36541Z",
     "iopub.status.busy": "2022-08-23T19:35:33.363203Z",
     "iopub.status.idle": "2022-08-23T19:35:33.444735Z"
    }
   },
   "source": [
    "Access to Brazil Data Cube (BDC) data products requires an access key. If you do not have a BDC `Access Token`, please go to [Brazil Data Cube Explorer](https://brazildatacube.dpi.inpe.br/portal/explore) and create an account. Then go to the management panel and generate a new `Access Token`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "access_token <- \"change-me\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we will create a query object called `stac_obj` pointing to the service address, allowing us to communicate with the `STAC` service."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "stac_obj <- stac(\"https://brazildatacube.dpi.inpe.br/stac/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#336699\">Listing the available Data Products</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To list all the image collections and data cube collections, we will make a request using the `get_request()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# query to the data catalog\n",
    "#\n",
    "catalog <- stac_obj %>% get_request()\n",
    "\n",
    "print(catalog)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# query the available product collections\n",
    "#\n",
    "collections <- stac_obj %>%\n",
    "    collections() %>%\n",
    "    get_request()\n",
    "\n",
    "print(collections, n = 31)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/brazil-data-cube/code-gallery/master/img/stac/stac-catalog.png?raw=true\" align=\"right\" width=\"300\"/>\n",
    "\n",
    "## <span style=\"color:#336699\">Retrieving the Metadata of a Collection</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">\n",
    "\n",
    "The `collections()` function returns information about a given image or data cube collection identified by its name. In this example, we are retrieving information about the datacube collection `CB4-16D-2`:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "collection_info <- stac_obj %>% \n",
    "    collections(\"CB4-16D-2\") %>% \n",
    "    get_request()\n",
    "\n",
    "print(collection_info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/brazil-data-cube/code-gallery/master/img/stac/stac-item.png?raw=true\" align=\"right\" width=\"300\"/>\n",
    "\n",
    "## <span style=\"color:#336699\">Retrieving Items</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">\n",
    "\n",
    "The `items()` function returns a query given a bounding box (`bbox`) and a date range (`datetime`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "items <- stac_obj %>% \n",
    "    collections(\"CB4-16D-2\") %>% \n",
    "    items(datetime = \"2018-08-01/2019-07-31\",\n",
    "          bbox  = c(-45.9, -12.9, -45.4, -12.6),\n",
    "          limit = 20) %>%\n",
    "    get_request() %>%\n",
    "    items_sign(sign_bdc(access_token))\n",
    "\n",
    "print(items)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's view the available bands that we can query using the `items_assets()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "items_assets(items)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/brazil-data-cube/code-gallery/master/img/stac/stac-asset.png?raw=true\" align=\"right\" width=\"300\"/>\n",
    "\n",
    "## <span style=\"color:#336699\">Assets</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">\n",
    "\n",
    "Assets are links to images, thumbnails or specific metadata files and can be accessed through the `assets` property (of an item):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can view the assets url using `assets_url()` as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "assets_url(items, asset_names = c(\"BAND14\", \"BAND13\"))[1:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#336699\">Reading and viewing the images</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will read and view the images through the terra package. First, we'll filter the assets corresponding to the date `2019-07-28`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# listing the datetime of all items \n",
    "#\n",
    "items_datetime(items)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# filtering the assets by a datetime\n",
    "#\n",
    "item_filtered <- items_filter(items, filter_fn = function(item) item$properties[[\"datetime\"]] == \"2019-07-28T00:00:00\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(item_filtered)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's read the images of the filtered item:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "blue_url  <- assets_url(item_filtered, asset_names = \"BAND13\", append_gdalvsi = TRUE)\n",
    "green_url <- assets_url(item_filtered, asset_names = \"BAND14\", append_gdalvsi = TRUE)\n",
    "red_url   <- assets_url(item_filtered, asset_names = \"BAND15\", append_gdalvsi = TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# reading the first images of each band\n",
    "#\n",
    "blue_rast  <- terra::rast(blue_url)\n",
    "green_rast <- terra::rast(green_url)\n",
    "red_rast   <- terra::rast(red_url)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to crop the image with a row and column window. To do this, we will create a `bbox` that satisfies this extension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "proj_orig <- sf::st_crs(\"+proj=longlat +datum=WGS84\")\n",
    "#BDC proj4 string\n",
    "proj_dest <- sf::st_crs(\"+proj=aea +lat_0=-12 +lon_0=-54 +lat_1=-2 +lat_2=-22 +x_0=5000000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs\")\n",
    "\n",
    "pts <- tibble::tibble(\n",
    "    lon = c(-45.89957, -45.40046),\n",
    "    lat = c(-12.9142, -12.58579)\n",
    ")\n",
    "pts_sf <- sf::st_as_sf(pts, coords = c(\"lon\", \"lat\"), crs = proj_orig)\n",
    "pts_transf <- sf::st_transform(pts_sf, crs = proj_dest)\n",
    "\n",
    "lat_dest <- sf::st_coordinates(pts_transf)[, 2]\n",
    "lon_dest <- sf::st_coordinates(pts_transf)[, 1]\n",
    "\n",
    "cat(\"Reprojected longitude:\", lon_dest, \"\\nReprojected latitude:\", lat_dest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# defining a clipping length\n",
    "#\n",
    "transformed_bbox <- terra::ext(5865751, 5920212, 9884783, 9920060)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# cropping images from an extent\n",
    "#\n",
    "blue_rast_cropped  <- terra::crop(blue_rast, transformed_bbox)\n",
    "green_rast_cropped <- terra::crop(green_rast, transformed_bbox)\n",
    "red_rast_cropped   <- terra::crop(red_rast, transformed_bbox)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualizing each band separately"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# setting plot display options\n",
    "options(repr.plot.width = 16, repr.plot.height = 5)\n",
    "par(mfrow = c(1, 3))\n",
    "\n",
    "plot(blue_rast_cropped,  main = \"Blue Band\")\n",
    "plot(green_rast_cropped, main = \"Green Band\")\n",
    "plot(red_rast_cropped,   main = \"Red Band\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#336699\">Composite Image Viewing</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a stack of bands for our composite plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# creating a band composition\n",
    "#\n",
    "rgb <- c(red_rast_cropped, green_rast_cropped, blue_rast_cropped)\n",
    "\n",
    "#\n",
    "# rgb view of the created composition\n",
    "#\n",
    "plotRGB(rgb, r = 1, g = 2, b = 3, stretch = \"lin\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#336699\">Calculating the Normalized Difference Vegetation Index (NDVI)</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **N**ormalized **D**ifference **V**egetation **I**ndex (NDVI) is calculated using the **Red** and **Near Infrared** (NIR) spectral bands. This index is used to assess whether or not the observed target contains live green vegetation. It can be calculated using the following equation:\n",
    "\n",
    "$$\n",
    "NDVI = \\frac{(NIR - RED)}{(NIR + RED)}\n",
    "$$\n",
    "\n",
    "<center><b>Equation 1</b> - NDVI</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<div style=\"text-align: justify;  margin-left: 15%; margin-right: 15%; border-style: solid; border-color: #0077b9; border-width: 1px; padding: 5px;\">\n",
    "<b>Note:</b>  Brazil Data Cube already provides for this data cube: <em>NDVI</em> and <em>EVI</em> along with the spectral bands. In addition, it also provides quality indicators (<em>CLEAROB</em>, <em>PROVENANCE</em>, <em>CMASK</em>, <em>TOTALOB</em>).\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, the bands `BAND15` and `BAND16` will be loaded from the filtered `items`.\n",
    "\n",
    "> As can be seen in the metadata of the `items`, `BAND15` corresponds to the **red** wavelength and `BAND16` to the **near-infrared**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the url of the **Red** band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "red <- assets_url(item_filtered, asset_names = \"BAND15\", append_gdalvsi = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the url of the **Near Infrared** band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "nir <- assets_url(item_filtered, asset_names = \"BAND16\", append_gdalvsi = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-12-13T04:19:48.308459Z",
     "iopub.status.busy": "2021-12-13T04:19:48.306379Z",
     "iopub.status.idle": "2021-12-13T04:19:48.322038Z"
    }
   },
   "source": [
    "Reading the cropped scene from the **Red (BAND15)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "red_rast <- terra::crop(terra::rast(red), transformed_bbox)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "red_rast"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reading the cropped scene from the **Near Infrared (BAND16)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "nir_rast <- terra::crop(terra::rast(nir), transformed_bbox)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "nir_rast"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's view the data that was loaded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "plot(red_rast)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "plot(nir_rast)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's calculate the **NDVI**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ndvi <- (nir_rast - red_rast) / (nir_rast + red_rast)\n",
    "ndvi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "plot(ndvi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#336699\">Image Thresholding</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the simplest approaches to separate different values in images is thresholding. This process consists of labeling the data based on fixed values.\n",
    "\n",
    "Let's try to separate our data into groups according to their NDVI values. To get started, first, let's see what the histogram of the image looks like. To do this we will use the `hist` function from the **terra** package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# creating the histogram\n",
    "#\n",
    "terra::hist(\n",
    "    ndvi, \n",
    "    xlim   = c(0, 1), \n",
    "    breaks = 10,  \n",
    "    main   = \"NDVI Distribution\",\n",
    "    xlab   = \"NDVI\", \n",
    "    ylab   = \"Frequency\", \n",
    "    col    = \"wheat\", \n",
    "    xaxt   = \"n\"\n",
    ")\n",
    "\n",
    "#\n",
    "# defining the interval on the x-axis\n",
    "#\n",
    "axis(side = 1, at = seq(0, 1, 0.1), labels = seq(0, 1, 0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming that we can separate the `ndvi` image with thresholding, we assume for this specific case that:\n",
    "\n",
    "* All pixels with values below 0.2 are dark pixels;\n",
    "* All pixels with values from 0.2 to 0.45 are sparsely vegetated areas.\n",
    "* All pixels above 0.45 are heavily vegetated areas;\n",
    "\n",
    "We can perform this thresholding by selecting from the `ndvi` matrix all values belonging to a given interval and assigning an integer value to it. We assume the following integer values:\n",
    "\n",
    "* `1`: Dark pixels;\n",
    "* `2`: Little vegetation\n",
    "* `3`: Abundant vegetation.\n",
    "\n",
    "To get started, we will first create a copy of the `ndvi` image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "labelled_img <- ndvi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# defining a vector of labels from the thresholds\n",
    "#\n",
    "vector_labels <- c(\n",
    "    0, 0.2, 1,    # valor 1\n",
    "    0.2, 0.45, 2, # valor 2\n",
    "    0.45, 1, 3    # valor 3\n",
    ")   \n",
    "\n",
    "#\n",
    "# transforming the label vector into an matrix\n",
    "#\n",
    "matrix_labels <- matrix(vector_labels, ncol = 3, byrow = TRUE)\n",
    "\n",
    "#\n",
    "# image with added thresholds\n",
    "#\n",
    "image_labelled <- terra::classify(labelled_img, matrix_labels, include.lowest = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets now see *labels*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width = 18, repr.plot.height = 10)\n",
    "\n",
    "#\n",
    "# Creating Figure\n",
    "#\n",
    "par(mfrow = c(1, 2))\n",
    "\n",
    "#\n",
    "# Plot of NDVI data\n",
    "#\n",
    "plot(ndvi, main = \"NDVI\")\n",
    "\n",
    "#\n",
    "# Plot of the thresholding result\n",
    "#\n",
    "plot(image_labelled, col = c(\"#E4E538\", \"#EFB17B\", \"#00AF22\"), main = \"Imagem rotulada\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#336699\">Calculating the difference between images</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compare the NDVI for images from two dates and the same location. This can be used, for example, to check areas where crops have grown and areas that have lost vegetation.\n",
    "\n",
    "For this calculation, we will use the NDVI indices provided in the data cube from two items with the same location and different dates using STAC.\n",
    "\n",
    "The first image comprises pixels from September 30, 2018 to October 15, 2018 (`2018-09-30_2018-10-15`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# filtering item by datetime\n",
    "#\n",
    "items_first <- items_filter(items, filter_fn = function(item) item$properties[[\"datetime\"]] == \"2018-09-30T00:00:00\")\n",
    "\n",
    "#\n",
    "# get the url from NDVI index\n",
    "#\n",
    "ndvi_first <- assets_url(items_first, asset_names = \"NDVI\", append_gdalvsi = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second selected image comprises pixels from January 1, 2019 to January 16, 2019 (`2019-01-01_2019-01-16`), that is, three months after the first selected image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# filtering item by datetime\n",
    "#\n",
    "items_second <- items_filter(items, filter_fn = function(item) item$properties[[\"datetime\"]] == \"2019-01-01T00:00:00\")\n",
    "\n",
    "#\n",
    "# get the url from NDVI index\n",
    "#\n",
    "ndvi_second <- assets_url(items_second, asset_names = \"NDVI\", append_gdalvsi = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"text-align: justify;  margin-left: 15%; margin-right: 15%; border-style: solid; border-color: #0077b9; border-width: 1px; padding: 5px;\">\n",
    "    <b>Nota:</b> The NDVI index pre-computed by BDC ranges from <em>-10000</em> to <em>10000</em>, instead of <em>-1</em> to <em>1</em>,  as can be seen in the `item` metadata. This is due to the smaller volume required to store files that use integer (16-bit) values rather than floating-point (32-bit) values.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considering that these images are of an agricultural area and that crops are usually planted near August (first observation), six months before or after the first observation, crops are expected to be found, implying higher NDVI values (more vigorous vegetation). This will cause the NDVI band to show brighter values in these areas. Using the gray color map, the high value NDVI pixels will be more like white, while the low-value NDVI pixels will be closer to black.\n",
    "\n",
    "Based on this, let's visually compare both NDVI images:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reading the cropped scene for the first NDVI date:\n",
    "> The region will be the **same** already used in the previous examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ndvi_first_rast <- terra::crop(terra::rast(ndvi_first), transformed_bbox)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reading the cropped scene for the second NDVI date:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ndvi_second_rast <- terra::crop(terra::rast(ndvi_second), transformed_bbox)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#\n",
    "# plot setup\n",
    "#\n",
    "par(mfrow = c(1, 2))\n",
    "\n",
    "#\n",
    "# NDVI data plot (First scene)\n",
    "#\n",
    "plot(ndvi_first_rast, main = \"First scene - 2018-09-30\")\n",
    "\n",
    "#\n",
    "# NDVI data plot (Second scene)\n",
    "#\n",
    "plot(ndvi_second_rast, main = \"Second scene - 2019-01-01\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we want to see what has grown and what has been lost, we will subtract the newest image from the oldest and plot it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ndvi_diff <- ndvi_second_rast - ndvi_first_rast"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualizing the difference between the two scenes with the `tmap` package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# \n",
    "# loading package\n",
    "#\n",
    "library(tmap)\n",
    "\n",
    "#\n",
    "# plot the difference between the two images\n",
    "#\n",
    "tm_shape(ndvi_diff) +\n",
    "  tm_raster(style = \"pretty\", palette = c(\"-RdYlBu\"), legend.hist = TRUE, midpoint = NA) + \n",
    "   tm_layout(legend.outside = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen in the NDVI difference graph, the main changes in pixel values were found in agricultural areas, which was expected due to crop changes.\n",
    "\n",
    "The blue values indicate negative values, while the red values are positive. This means that there was a loss of vegetation for the blue areas, as a decreasing result in the NDVI value, which means that crops were harvested. Meanwhile, in the red areas, the NDVI value has increased due to more vigorous vegetation in the most recent dates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#336699\">References</span>\n",
    "<hr style=\"border:1px solid #0077b9;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- [Spatio Temporal Asset Catalog Specification](https://stacspec.org/)\n",
    "\n",
    "\n",
    "- [Brazil Data Cube R Client Library for STAC Service - GitHub Repository](https://github.com/brazil-data-cube/rstac)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}