landsat8-cog-ndvi.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Landsat 8 NDVI Analysis on the Cloud\n",
    "\n",
    "This notebook demonstrates a \"Cloud Native\" analysis of [Normalized Difference Vegetation Index (NDVI)](https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index) using [Landsat 8 satellite imagery](https://landsat.usgs.gov/landsat-8). \n",
    "\n",
    "**To run, simply hit Shift+Enter to run each code block**\n",
    "\n",
    "By \"Cloud Native\" we mean that images are not downloaded to your local machine. Instead, calculations are performed efficiently in parallel across many distributed machines on Google Cloud ([where the imagery is stored](https://cloud.google.com/storage/docs/public-datasets/landsat)).\n",
    "\n",
    "This workflow is possible because the Landsat 8 data is stored in [Cloud-Optimized Geotiff](http://www.cogeo.org) format, which can be accessed remotely via [xarray](http://xarray.pydata.org/en/stable/) and [rasterio](https://rasterio.readthedocs.io/en/latest/) Python libraries. Interactive, dynamically updating visualization is done with [Holoviews](http://holoviews.org). Distributed computing is enabled through a [Pangeo](http://pangeo.io) JupyterHub deployment with [Dask Kubernetes](https://github.com/dask/dask-kubernetes).\n",
    "\n",
    "Created on August 30, 2018 by:\n",
    "\n",
    "Scott Henderson, Daniel Rothenberg, Matthew Rocklin, Rich Signell, Joe Hamman, Ryan Abernathey, and Rob Fatland"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import required libraries\n",
    "import os\n",
    "import pandas as pd\n",
    "import rasterio\n",
    "import xarray as xr\n",
    "import requests\n",
    "import geoviews as gv\n",
    "import holoviews as hv\n",
    "import hvplot.xarray\n",
    "import hvplot.pandas\n",
    "import shapely\n",
    "\n",
    "import dask\n",
    "from dask_kubernetes import KubeCluster\n",
    "from dask.distributed import Client\n",
    "from dask.distributed import wait, progress\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print package versions\n",
    "print('Xarray version: ', xr.__version__)\n",
    "print('Rasterio version: ', rasterio.__version__)\n",
    "print('dask version: ', dask.__version__)\n",
    "print('hvplot version: ', hvplot.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set environment variables for cloud-optimized-geotiffs efficiency\n",
    "os.environ['GDAL_DISABLE_READDIR_ON_OPEN']='YES'\n",
    "os.environ['CPL_VSIL_CURL_ALLOWED_EXTENSIONS']='TIF'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use NASA Common Metadata Repository (CMR) to get Landsat 8 images\n",
    "\n",
    "[NASA CMR](https://earthdata.nasa.gov/about/science-system-description/eosdis-components/common-metadata-repository) is a new unified way to search for remote sensing assests across many archive centers. If you prefer a graphical user interface, NASA [Earthdata Search](https://search.earthdata.nasa.gov/search) is built on top of CMR. CMR returns download links through the USGS (https://earthexplorer.usgs.gov), but the same archive is mirrored as a (Google Public Dataset)[https://cloud.google.com/storage/docs/public-datasets/landsat], so we'll make a function that queries CMR and returns URLs to the imagery stored on Google Cloud."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def query_cmr_landsat(collection='Landsat_8_OLI_TIRS_C1',tier='T1', path=47, row=27):\n",
    "    \"\"\"Query NASA CMR for Collection1, Tier1 Landsat scenes from a specific path and row.\"\"\"\n",
    "    \n",
    "    data = [f'short_name={collection}',\n",
    "            f'page_size=2000',\n",
    "            f'attribute[]=string,CollectionCategory,{tier}',\n",
    "            f'attribute[]=int,WRSPath,{path}',\n",
    "            f'attribute[]=int,WRSRow,{row}',\n",
    "           ]\n",
    "\n",
    "    query = 'https://cmr.earthdata.nasa.gov/search/granules.json?' + '&'.join(data)\n",
    "\n",
    "    r = requests.get(query, timeout=100)\n",
    "    print(r.url)\n",
    "    \n",
    "    df = pd.DataFrame(r.json()['feed']['entry'])\n",
    "    \n",
    "    # Save results to a file\n",
    "    #print('Saved results to cmr-result.json')\n",
    "    #with open('cmr-result.json', 'w') as j:\n",
    "    #    j.write(r.text)\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_google_archive(pids, bands):\n",
    "    \"\"\"Turn list of product_ids into pandas dataframe for NDVI analysis.\"\"\"\n",
    "    \n",
    "    path =  pids[0].split('_')[2][1:3]\n",
    "    row =  pids[0].split('_')[2][-2:]\n",
    "    baseurl = f'https://storage.googleapis.com/gcp-public-data-landsat/LC08/01/0{path}/0{row}'\n",
    "    \n",
    "    dates = [pd.to_datetime(x.split('_')[3]) for x in pids]\n",
    "    df = pd.DataFrame(dict(product_id=pids, date=dates))\n",
    "    \n",
    "    for band in bands:\n",
    "        df[band] = [f'{baseurl}/{x}/{x}_{band}.TIF' for x in pids]\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Landsat scenes over a specific area are categorized by Path and Row\n",
    "df = query_cmr_landsat(path=47, row=27)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Holoviews maps are interactive, you can zoom in with the magnifying glass tool\n",
    "coordlist = df.polygons.iloc[0]\n",
    "lats = coordlist[0][0].split(' ')[::2]\n",
    "lons = coordlist[0][0].split(' ')[1::2]\n",
    "coords = [(float(lon),float(lat)) for lon, lat in zip(lons,lats)]\n",
    "poly = shapely.geometry.Polygon(coords)\n",
    "buffer = 1 #degrees\n",
    "left, bottom, right, top  = poly.bounds\n",
    "\n",
    "footprint = gv.Shape(poly, label=df.title.iloc[0]).options(alpha=0.5)\n",
    "tiles = gv.tile_sources.CartoEco.options(width=700, height=500).redim.range(Latitude=(bottom-1, top+1), Longitude=(left-1,right+1)) \n",
    "labels = gv.tile_sources.StamenLabels.options(level='annotation')\n",
    "tiles * footprint * labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all scenes for a given path and row, bands 4 and 5\n",
    "pids = df.title.tolist()\n",
    "# Don't use the most recent date since there can be a lag in data being on Google Storage\n",
    "df = make_google_archive(pids[:-1], ['B4', 'B5'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df.head()\n",
    "df.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Launch Dask Kubernetes Cluster\n",
    "\n",
    "This will allow us to distribute our analysis across many machines. In the default configuration for Pangeo Binder, each worker has 2 vCPUs and 7Gb of RAM. It may take several minutes to initialize these workers and make them available to Dask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select 10 'workers' under 'manual scaling' menu below and click 'Scale'\n",
    "# Click on the 'Dashboard link' to monitor calculation progress\n",
    "cluster = KubeCluster(n_workers=10)\n",
    "cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Attach Dask to the cluster\n",
    "client = Client(cluster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examine a single band Landsat image\n",
    "\n",
    "The *rasterio* library allows us to read Geotiffs on the web without downloading the entire image. *Xarray* has a built-in load_rasterio() function that allows us to open the file as a DataArray. Xarray also uses Dask for lazy reading, so we want to make sure the native block tiling of the image matches the dask \"chunk size\". These dask chunks are automatically distributed among all our workers when a computation is requested, so ideally they will fit in the worker memory. A chunk size of 2048x2048 with a float32 datatype implies a 16Mb array.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load with rasterio\n",
    "image_url = df.iloc[0]['B4']\n",
    "with rasterio.open(image_url) as src:\n",
    "    print(src.profile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that the blocksize of the image is 256 by 256, so we want xarray to use some multiple of that\n",
    "xchunk = 2048\n",
    "ychunk = 2048\n",
    "da = xr.open_rasterio(image_url, chunks={'band': 1, 'x': xchunk, 'y': ychunk})\n",
    "da"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### holoviews visualization\n",
    "If we request to compute something or plot these arrays, the necessary data chunks will be accessed on cloud storage. Watch the KubeCluster dashboard to see the worker activity when this command is run. Note that no data is stored on the disk here, it's all in memory\n",
    "\n",
    "Use the magnifying glass button on the right to interactively zoom in. The image resolution automatically updates based on zoom level. The cursor gives you UTM coordinates and the image value at that point!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band1 = da.sel(band=1).persist()\n",
    "img = band1.hvplot(rasterize=True, dynamic=True, width=700, height=500, cmap='magma')\n",
    "img"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load all Landsat bands into an xarray dataset\n",
    "\n",
    "Often we want to analyze a time series of satellite imagery, but we are constrained by computational resources. So we either download all the images, extract a small subset and then do our analysis. Or, we coarsen the resolution of all our images so that the entire set fits into our computer RAM. Because this notebook is running on Google Cloud with access to many resources in our Kube Cluster, we no longer have to worry about the computational constraints, and can conduct our analysis at full resoution!\n",
    "\n",
    "First we need to construct an xarray dataset object (which has data variables 'band4' and 'band5' in a n-dimensional array with x-coordinates representing UTM easting, y-coordinates representing UTM northing, and a time coordinate representing the image acquisition date).\n",
    "\n",
    "There are different ways to go about this, but we will load our images with a timestamp index since each image is taken on a different date. Typically, this is a chore if our images are not on the same grid to begin with, but xarray knows how to automatically align images based on their georeferenced coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that these landsat images are not necessarily the same shape or on the same grid:\n",
    "for image_url in df.B4[:5]:\n",
    "    with rasterio.open(image_url) as src:\n",
    "        print(src.shape, src.bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_multiband_dataset(row, bands=['B4','B5'], chunks={'band': 1, 'x': 2048, 'y': 2048}):\n",
    "    '''A function to load multiple landsat bands into an xarray dataset '''\n",
    "    \n",
    "    # Each image is a dataset containing both band4 and band5\n",
    "    datasets = []\n",
    "    for band in bands:\n",
    "        url = row[band]\n",
    "        da = xr.open_rasterio(url, chunks=chunks)\n",
    "        da = da.squeeze().drop(labels='band')\n",
    "        ds = da.to_dataset(name=band)\n",
    "        datasets.append(ds)\n",
    "\n",
    "    DS = xr.merge(datasets)\n",
    "    \n",
    "    return DS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Merge all acquisitions into a single large Dataset, this will take a minute\n",
    "datasets = []\n",
    "for i,row in df.iterrows():\n",
    "    try:\n",
    "        print('loading...', row.date)\n",
    "        ds = create_multiband_dataset(row)\n",
    "        datasets.append(ds)\n",
    "    except Exception as e:\n",
    "        print('ERROR loading, skipping acquistion!')\n",
    "        print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an xarray dataset\n",
    "DS = xr.concat(datasets, dim=pd.DatetimeIndex(df.date.tolist(), name='time'))\n",
    "print('Dataset size (Gb): ', DS.nbytes/1e9)\n",
    "DS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note that xarray has automatically expanded the dimensions to include the maximum extents of all the images, also the chunksize has been automatically adjusted.\n",
    "\n",
    "There is definitely some room for improvement here from a computational efficiency standpoint - in particular the dask chunks are no longer aligned with the image tiles. This is because each image starts at different coordinates and has different shapes, but xarray uses a single chunk size for the entire datasets. There will also be many zeros in this dataset, so future work could take advantage of sparse arrays. \n",
    "\n",
    "These points aside, our KubeCluster will automatically parallelize our computations for us, so we can not worry too much about optimal efficiency and just go ahead and run our analysis!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Again, only metadata is retrieved at this point, wich is why it's so quick!\n",
    "da = DS.sel(time='2013-04-21')['B4']\n",
    "print('Image size (Gb): ', da.nbytes/1e9)\n",
    "da"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distributed NDVI computations\n",
    "\n",
    "Set up our NDVI dataset. Note that NDVI is not actually computed until we call the Dask compute(), persist(), or call other functions such as plot() that require actually operate on the data!\n",
    "\n",
    "Because we now have a timeseries of images Holoviews automatically will add a time slider to our visualization! We add a 'SingleTap' features as well to keep track of interesting coordinates (by clicking on the image a white dot will appear and the coordinates are stored in the 'taps' list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NDVI = (DS['B5'] - DS['B4']) / (DS['B5'] + DS['B4'])\n",
    "NDVI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "taps = []\n",
    "def record_coords(x, y):\n",
    "    if None not in [x,y]:\n",
    "        taps.append([x, y])\n",
    "    return hv.Points(taps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: this will take a minute to load and is best viewed on a wide monitor\n",
    "# the time slider can get hidden on small screens\n",
    "img = NDVI.hvplot('x', 'y', groupby='time', dynamic=True, rasterize=True, width=700, height=500, cmap='magma')\n",
    "tap = hv.streams.SingleTap(transient=True, source=img)\n",
    "clicked_points = hv.DynamicMap(record_coords, streams=[tap])\n",
    "\n",
    "img * clicked_points.options(size=10, color='w')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Points clicked are stored in the 'taps list'\n",
    "if len(taps) == 0:\n",
    "    taps = [(562370, 5312519)]\n",
    "\n",
    "print('Selected points:')\n",
    "taps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract time series for region around point selected from map\n",
    "\n",
    "This uses a buffer around a selected point and does monthly resampling and will probably take a minute or so to pull the necessary data and run computations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xcen,ycen = taps[0]\n",
    "buf = 5000  # look at point +/- 5km\n",
    "ds = NDVI.sel(x=slice(xcen-buf,xcen+buf), y=slice(ycen-buf,ycen+buf))\n",
    "timeseries = ds.resample(time='1MS').mean(dims=['x','y']).persist()\n",
    "# Store as pandas series\n",
    "s = timeseries.to_series() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Holoviews is also great for interative 2D plots\n",
    "line = s.hvplot(width=700, height=300, legend=False)\n",
    "points = s.hvplot.scatter(width=700, height=300, legend=False)\n",
    "label = f'Mean NDVI: easting={xcen:g} , northing={ycen:g}'\n",
    "\n",
    "(line * points).relabel(label)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot subset of selected region at full resolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.sel(time=slice('2015-01-01', '2015-06-15')).plot.imshow('x', 'y', col='time', col_wrap=4, cmap='magma', vmin=0, vmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In conclusion\n",
    "\n",
    "* This notebook demonstrates the power of storing data publically in the Cloud as optimized geotiffs - scientists can conduct scalable analysis without downloading the data to a local machine. Only derived subsets and figures need to be downloaded!\n",
    "* We used a crude NDVI calculation, designed to demonstrate the syntax and tools - a proper analysis should take into account cloud masks and other corrections"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}