{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To start, we install Rasterio, a Python module for interacting with gridded spatial data\n",
    "!pip install rasterio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reading and manipulating tiled GeoTIFF datasets\n",
    "\n",
    "This notebook shows how to perform simple calculations with a GeoTIFF dataset using XArray and Dask. We load and rescale a Landsat 8 image and compute NDVI (Normalized difference vegetation index). This can be used to distinguish green vegetation from areas of bare land or water.\n",
    "\n",
    "We'll use an image of the Denver, USA area taken in July 2018."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![RGB image](https://landsat-pds.s3.amazonaws.com/c1/L8/033/033/LC08_L1TP_033033_20180706_20180717_01_T1/LC08_L1TP_033033_20180706_20180717_01_T1_thumb_small.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download data\n",
    "\n",
    "First, we download the dataset. We are using an image from the cloud-hosted [Landsat 8 public dataset](https://landsatonaws.com/) and each band is available as a separate GeoTIFF file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import json\n",
    "import rasterio\n",
    "import requests\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nir_filename = 'https://landsat-pds.s3.amazonaws.com/c1/L8/033/033/LC08_L1TP_033033_20180706_20180717_01_T1/LC08_L1TP_033033_20180706_20180717_01_T1_B5.TIF'\n",
    "red_filename = 'https://landsat-pds.s3.amazonaws.com/c1/L8/033/033/LC08_L1TP_033033_20180706_20180717_01_T1/LC08_L1TP_033033_20180706_20180717_01_T1_B4.TIF'\n",
    "mtl_filename = 'https://landsat-pds.s3.amazonaws.com/c1/L8/033/033/LC08_L1TP_033033_20180706_20180717_01_T1/LC08_L1TP_033033_20180706_20180717_01_T1_MTL.json'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def download_file(in_filename, out_filename):\n",
    "    if not os.path.exists(out_filename):\n",
    "        print(\"Downloading\", in_filename)\n",
    "        response = requests.get(in_filename)\n",
    "        with open(out_filename, 'wb') as f:\n",
    "            f.write(response.content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "download_file(nir_filename, 'nir.tif')\n",
    "download_file(red_filename, 'red.tif')\n",
    "download_file(mtl_filename, 'meta.json')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check image metadata\n",
    "\n",
    "Let's see if the image is tiled so we can select a chunk size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = rasterio.open('red.tif')\n",
    "print(img.is_tiled)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img.block_shapes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image has separate blocks for each band with block size 512 x 512. \n",
    "\n",
    "## Create XArray datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "red = xr.open_rasterio('red.tif', chunks={'band': 1, 'x': 1024, 'y': 1024})\n",
    "nir = xr.open_rasterio('nir.tif', chunks={'band': 1, 'x': 1024, 'y': 1024})\n",
    "nir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each dataset's data is a Dask array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red.variable.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optional: create a Dask client\n",
    "\n",
    "You can start a Dask client to monitor execution with the dashboard."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dask\n",
    "from dask.distributed import Client\n",
    "client = Client(processes=False)\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rescale bands using Landsat metadata\n",
    "\n",
    "The Landsat Level 1 images are delivered in a quantized format. This has to be [converted to top-of-atmosphere reflectance](https://landsat.usgs.gov/using-usgs-landsat-8-product) using the provided metadata.\n",
    "\n",
    "First we define convenience functions to load the rescaling factors and transform a dataset. The red band is band 4 and near infrared is band 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_scale_factors(filename, band_number):\n",
    "    with open(filename) as f:\n",
    "        metadata = json.load(f)\n",
    "    M_p = metadata['L1_METADATA_FILE'] \\\n",
    "                  ['RADIOMETRIC_RESCALING'] \\\n",
    "                  ['REFLECTANCE_MULT_BAND_{}'.format(band_number)]\n",
    "    A_p = metadata['L1_METADATA_FILE'] \\\n",
    "                  ['RADIOMETRIC_RESCALING'] \\\n",
    "                  ['REFLECTANCE_ADD_BAND_{}'.format(band_number)]\n",
    "    return M_p, A_p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_reflectance(ds, band_number, metafile='meta.json'):\n",
    "    M_p, A_p = load_scale_factors(metafile, band_number)\n",
    "    toa = M_p * ds + A_p\n",
    "    return toa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red_toa = calculate_reflectance(red, band_number=4)\n",
    "nir_toa = calculate_reflectance(nir, band_number=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the transformation is composed of arithmetic operations, execution is delayed and the operations are parallelized automatically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(red_toa.variable.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting image has floating point data with magnitudes appropriate to reflectance. This can be checked by computing the range of values in an image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red_max, red_min, red_mean = dask.compute(\n",
    "    red_toa.max(dim=['x', 'y']), \n",
    "    red_toa.min(dim=['x', 'y']),\n",
    "    red_toa.mean(dim=['x', 'y'])\n",
    ")\n",
    "print(red_max.item())\n",
    "print(red_min.item())\n",
    "print(red_mean.item())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate and display NDVI\n",
    "\n",
    "Now that we have the image as reflectance values, we are ready to compute NDVI.\n",
    "\n",
    "$$\n",
    "NDVI = \\frac{NIR - Red}{NIR + Red}\n",
    "$$\n",
    "\n",
    "This highlights areas of healthy vegetation with high NDVI values, which appear as green in the image below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi = (nir_toa - red_toa) / (nir_toa + red_toa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi2d = ndvi.squeeze()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "im = ndvi2d.compute().plot.imshow(cmap='BrBG', vmin=-0.5, vmax=1)\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}