03_multidim_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-dimensional Analysis with xarray\n",
    "\n",
    "### Questions\n",
    "- What Python tools can I use work with multidimensional data like NetCDF files? \n",
    "\n",
    "### Objectives\n",
    "- Learn how to use xarray to conscisely work with multidimensional data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Introduction to multidimensional arrays\n",
    "\n",
    "Unlabelled, N-dimensional arrays of numbers, such as NumPy's ndarray, are the most widely used data structure in scientific computing. Geoscientists have a particular need for structuring their data as arrays. For example, we commonly work with sets of climate variables (e.g. temperature and precipitation) that vary in space and time and are represented on a regularly-spaced grid. Often we need to subset a large global grid to look at data for a particular region, or select a specific time slice. Then we might want to apply statistical functions to these subsetted groups to generate summary information.\n",
    "\n",
    "<br>\n",
    "<img src=\"http://xarray.pydata.org/en/stable/_images/dataset-diagram.png\" width = \"800\" border = \"10\">\n",
    "<br>\n",
    "\n",
    "> ## Isn't this the same as raster processing?\n",
    "> The tools in this tutorial have some similarity to raster image processing tools.\n",
    "> Both require computational engines that can manipulate large stacks of data formatted as arrays.\n",
    "> Here we focus on tools that are optimized to handle data that have many variables spanning dimensions\n",
    "> of time and space. See the raster tutorials for tools that are optimized for image processing of remote sensing datasets.\n",
    "\n",
    "\n",
    "### Conventional Approach: Working with Unlabelled Arrays\n",
    "\n",
    "Multidimensional array data are often stored in user-defined binary formats, and distributed with custom Fortran\n",
    "or C++ libraries used to read and process the data. Users are responsible for setting up their own file structures and custom codes to handle these files. Subsetting the data involves reading everything into an in-memory array, and then using a series of nested loops with conditional statements to look for a specific range of index values associated with the temporal or spatial slice needed. Also, clever use of matrix algebra is often used to summarize data across spatial and temporal dimensions.\n",
    "\n",
    "### Challenges:\n",
    "\n",
    "The biggest challenge in working with N-dimensional arrays in this fashion is the fact that the data are almost disassociated from their metadata. Users are left with the task of tracking the meaning behind array indices using domain-specific software, often leading to inefficiencies and errors. Common pitfalls often occur in in the form of questions like \"is the time axis of my array in the first or third index position?\", or \"does my array of timestamps still align with my data after resampling?\"."
   ]
  },
  {
   "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": [
    "import xarray as xr "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Recall that we are interested in looking at landuse over the State of Pará in Brazil, where extensive logging and illegal deforestation is happening. The Landsat tile we will be looking at is Path 227, Row 065. The date for the file we will be accessing is 8 June, 2020 and we will extract the NIR, red band and metadata file from the AWS s3 bucket**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open path to file on s3 bucket with rasterio\n",
    "print('Landsat on AWS:')\n",
    "filepath = 'http://landsat-pds.s3.amazonaws.com/c1/L8/227/065/LC08_L1TP_227065_20200608_20200626_01_T1/LC08_L1TP_227065_20200608_20200626_01_T1_B4.TIF'\n",
    "with rasterio.open(filepath) as src:\n",
    "    print(src.profile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "date = '2020-06-08'\n",
    "url = 'http://landsat-pds.s3.amazonaws.com/c1/L8/227/065/LC08_L1TP_227065_20200608_20200626_01_T1/'\n",
    "redband = 'LC08_L1TP_227065_20200608_20200626_01_T1_B{}.TIF'.format(4)\n",
    "nirband = 'LC08_L1TP_227065_20200608_20200626_01_T1_B{}.TIF'.format(5)\n",
    "mtlfile = 'LC08_L1TP_227065_20200608_20200626_01_T1_{}.json'.format('MTL')\n",
    "\n",
    "# Overviews are reduced resolution versions of your dataset that can speed up rendering when you don’t need \n",
    "# full resolution. By precomputing the upsampled pixels, rendering can be significantly faster when zoomed out.\n",
    "# More info here: https://rasterio.readthedocs.io/en/latest/topics/overviews.html\n",
    "\n",
    "with rasterio.open(url+redband) as src:\n",
    "    profile = src.profile\n",
    "    oviews = src.overviews(1) # list of overviews from biggest to smallest\n",
    "    oview = oviews[1]  # Use second-highest resolution overview\n",
    "    print('Decimation factor= {}'.format(oview))\n",
    "    red = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's look at what the red band image looks like\n",
    "\n",
    "plt.imshow(red)\n",
    "plt.colorbar()\n",
    "plt.title('{}\\nRed {}'.format(redband, red.shape))\n",
    "plt.xlabel('Column #')\n",
    "plt.ylabel('Row #')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's download the files \n",
    "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": [
    "nir_filename = url + nirband\n",
    "red_filename = url + redband\n",
    "mtl_filename = url + mtlfile\n",
    "\n",
    "download_file(nir_filename, 'nir.tif')\n",
    "download_file(red_filename, 'red.tif')\n",
    "download_file(mtl_filename, 'meta.json')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the shape size for the red band image\n",
    "\n",
    "red = rasterio.open(red_filename)\n",
    "print(red.is_tiled)\n",
    "red.block_shapes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Technically you do not need to download the files to read the data\n",
    "\n",
    "red = xr.open_rasterio(red_filename, chunks={'band': 1, 'x': 1024, 'y': 1024})\n",
    "nir = xr.open_rasterio(nir_filename, chunks={'band': 1, 'x': 1024, 'y': 1024})\n",
    "red"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  <font color=blue>Exercise 1: Inspect the dimensions in the dataset. How would you only view the dimensions along the x dimension? </font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer for Exercise 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate and plot NDVI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Landsat Level 1 images are delivered in a quantized format. This has to be converted to top-of-atmosphere reflectance 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.\n",
    "\n",
    "We will also introduce here the concept of Dask, a flexible library for parallel computing in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "N𝐷𝑉𝐼=𝑁𝐼𝑅−𝑅𝑒𝑑 / 𝑁𝐼𝑅+𝑅𝑒𝑑"
   ]
  },
  {
   "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": "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(red_toa.variable.data)"
   ]
  },
  {
   "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": "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}