{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background\n",
    "\n",
    "Analysis ready data (ARD) is data that has been processed to improve ease-of-use for an end user. This includes the process of applying corrections, for example, removing the atmospheric effects that appear in imagery. In this notebook, we work with Landsat data and will comparing unprocessed (Level-1) data with processed (Level-2, or analysis ready) data. \n",
    "\n",
    "When working with optical data, a common analysis method is to calculate various band indices. These are typically constructed as the ratio between the sum of two bands and their difference:\n",
    "$$\\text{Index} = \\frac{\\text{Band 1} - \\text{Band 2}}{\\text{Band 1} + \\text{Band 2}}.$$\n",
    "In particular, band indices are commonly used to extract information on vegetation coverage or classify water pixels. \n",
    "\n",
    "Given the level of processing required, it's worth understanding whether band-index analyses can be carried out with Level-1 data as opposed to Level-2 data. In this notebook, we calculate key band indicies for Landsat Level-1 and Level-2 products and compare them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import datacube\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from utils.utils import lat_lon_to_epsg, three_band_image\n",
    "from utils.dc_display_map import display_map\n",
    "\n",
    "from ipywidgets import interact, interactive\n",
    "from IPython.display import clear_output, display, HTML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Available Landsat products\n",
    "The `list_products` method in the Datacube class displays the names and details of all available products. In the below cell we will query what Landsat 8 Products are currently indexed in our instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Connect to a datacube\n",
    "dc = datacube.Datacube(app='Level 1 and Level 2 Comparison')\n",
    "\n",
    "# List metadata for all Landsat NBAR and NBART products available in DEA\n",
    "dc_products = dc.list_products()\n",
    "display_columns = ['name', 'description', 'product_type']\n",
    "dc_products[dc_products['name'].str.contains(\"ls8|ls8\")][display_columns].set_index('name')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Level-1 data has product name `ls8_level1_usgs`. Level-2 data has product name `ls8_usgs_sr_scene`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualise and load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The coordinates below correspond to an area containing Lake Rweru,\n",
    "# which is on the boarder of Rwanda and Burundi\n",
    "latitude = (-2.2, -2.45)\n",
    "longitude  = (30.1, 30.35)\n",
    "\n",
    "# Date range to load over\n",
    "date_range = (\"2018-01-01\", \"2018-02-01\")\n",
    "\n",
    "# Display the map before loading the data\n",
    "display_map(latitude, longitude)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the EPSG of a WGS UTM coordinate reference system that is appropriate for our data\n",
    "EPSG = lat_lon_to_epsg(latitude[1], longitude[1])\n",
    "\n",
    "# Load the Level-1 data\n",
    "data_cube_level1 = dc.load(\n",
    "    product='ls8_level1_usgs',\n",
    "    x=longitude,\n",
    "    y=latitude, \n",
    "    output_crs='epsg:' + EPSG,\n",
    "    resolution=(-30, 30),\n",
    "    time = date_range,\n",
    "    measurements = ['nir', 'red', 'swir1', 'green', 'blue']\n",
    ")\n",
    "print(data_cube_level1)\n",
    "\n",
    "# Load the Level-2 data\n",
    "data_cube_level2 = dc.load(\n",
    "    product='ls8_usgs_sr_scene',\n",
    "    x=longitude,\n",
    "    y=latitude, \n",
    "    output_crs='epsg:' + EPSG,\n",
    "    resolution=(-30, 30),\n",
    "    time = date_range,\n",
    "    measurements = ['nir', 'red', 'swir1', 'green', 'blue']\n",
    ")\n",
    "print(data_cube_level2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate various band indices\n",
    "\n",
    "For the demonstrative purpose of this notebook, we calculate three indices:\n",
    "\n",
    "### NDVI\n",
    "The normalised difference vegetation index (NDVI) indicates the presence of green vegetation. It is calculated as\n",
    "$$\\text{NDVI} = \\frac{\\text{NIR} - \\text{Red}}{\\text{NIR} + \\text{Red}}$$\n",
    "High values indicate dense vegetation.\n",
    "\n",
    "### NDWI\n",
    "The normalised difference water index (NDWI) indicates the presence of water. It is calculated as\n",
    "$$\\text{NDWI} = \\frac{\\text{Green} - \\text{NIR}}{\\text{Green} + \\text{NIR}}$$\n",
    "High values indicate water.\n",
    "\n",
    "### MNDWI\n",
    "The modified normalised difference water index (MNDWI) indicates the presence of water. It is calculated as\n",
    "$$\\text{MNDWI} = \\frac{\\text{Green} - \\text{SWIR}}{\\text{Green} + \\text{SWIR}}$$\n",
    "High values indicate water. This index provides enhanced water classification relative to the NDWI."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following cell, we define a function to calculate and add each index to a given data cube object. The band indices are stored as individual data varaiables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_band_indices(dc):\n",
    "    dc['ndvi'] = (dc.nir - dc.red)/(dc.nir + dc.red)\n",
    "    dc['ndwi'] = (dc.green - dc.nir)/(dc.green + dc.nir)\n",
    "    dc['mndwi'] = (dc.green - dc.swir1)/(dc.nir + dc.swir1)\n",
    "    return(dc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dc_level1 = calculate_band_indices(data_cube_level1)\n",
    "dc_level2 = calculate_band_indices(data_cube_level2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot indices for same time stamp\n",
    "\n",
    "We now show the index values for the same time step. Band indices are designed to occur on a scale of -1 to 1, so we keep these as the bounds for the colour maps in both images. Consequently, the intensity of the colour can be directly compared.\n",
    "\n",
    "Run the following commands to see the available time-steps for the Level-1 and Level-2 data cube objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(dc_level1.time)\n",
    "print(dc_level2.time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by stating which time steps to use for plotting. Choose index values that give the same times, noting that the times may not perfectly match between the Level-1 and Level-2 data cubes.\n",
    "\n",
    "The choices below correspond to the same time-stamp in the Level-1 and Level-2 data sets, which is an image with minimal cloud coverage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dc_level1_plotindex = 4\n",
    "dc_level2_plotindex = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDVI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(18,6))\n",
    "fig.suptitle(\"Comparison of NDVI between Level-1 (LHS) and Level-2 (RHS)\")\n",
    "\n",
    "plt.subplot(121)\n",
    "dc_level1.ndvi.isel(time=dc_level1_plotindex).plot(cmap='RdYlGn', vmin=-1, vmax=1)\n",
    "\n",
    "plt.subplot(122)\n",
    "dc_level2.ndvi.isel(time=dc_level2_plotindex).plot(cmap='RdYlGn', vmin=-1, vmax=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDWI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig2 = plt.figure(figsize=(18,6))\n",
    "fig2.suptitle(\"Comparison of NDWI between Level-1 (LHS) and Level-2 (RHS)\")\n",
    "\n",
    "plt.subplot(121)\n",
    "dc_level1.ndwi.isel(time=dc_level1_plotindex).plot(cmap='GnBu', vmin=-1, vmax=1)\n",
    "\n",
    "plt.subplot(122)\n",
    "dc_level2.ndwi.isel(time=dc_level2_plotindex).plot(cmap='GnBu', vmin=-1, vmax=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MNDWI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig3 = plt.figure(figsize=(18,6))\n",
    "fig3.suptitle(\"Comparison of MNDWI between Level-1 (LHS) and Level-2 (RHS)\")\n",
    "\n",
    "plt.subplot(121)\n",
    "dc_level1.mndwi.isel(time=dc_level1_plotindex).plot(cmap='RdBu', vmin=-1, vmax=1)\n",
    "\n",
    "plt.subplot(122)\n",
    "dc_level2.mndwi.isel(time=dc_level2_plotindex).plot(cmap='RdBu', vmin=-1, vmax=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot histograms of difference in index values between Level-1 and Level-2\n",
    "\n",
    "We can assess the distribution of the difference in the index values by taking a histogram of the difference in the index values. A value of 0 difference would indicate that the value of the index is the same in the Level 1 and Level 2 data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDVI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_diff = dc_level2.ndvi.isel(time=dc_level2_plotindex) - dc_level1.ndvi.isel(time=dc_level1_plotindex)\n",
    "\n",
    "fig4 = plt.figure(figsize=(6,6))\n",
    "ndvi_diff.plot.hist(bins=100)\n",
    "plt.title('Histogram of difference in NDVI (Level-2 - Level-1)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDWI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndwi_diff = dc_level2.ndwi.isel(time=dc_level2_plotindex) - dc_level1.ndwi.isel(time=dc_level1_plotindex)\n",
    "\n",
    "fig5 = plt.figure(figsize=(6,6))\n",
    "ndwi_diff.plot.hist(bins=100)\n",
    "plt.title('Histogram of difference in NDWI (Level-2 - Level-1)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MNDWI\n",
    "\n",
    "We note that there are some extreme values in the Level 1 data due to the presence of clouds, so we don't plot the difference values that are smaller than -2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mndwi_diff = dc_level2.mndwi.isel(time=dc_level2_plotindex) - dc_level1.mndwi.isel(time=dc_level1_plotindex)\n",
    "filtered_mndwi_diff = mndwi_diff.where(mndwi_diff > -2)\n",
    "\n",
    "fig6 = plt.figure(figsize=(6,6))\n",
    "filtered_mndwi_diff.plot.hist(bins=100)\n",
    "plt.title('Histogram of difference in MNDWI (Level-2 - Level-1)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "The visual inspection and histograms reveal that the band indices calculated for Level-1 and Level-2 data differ for the vast majority of pixels. Clouds also have spurious index values in the Level-1 data, which might also affect analyses. We conclude that you could not necessarily use the index values from Level-1 in place of those from Level-2."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}