1-multiple-cogs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Best-practices for Cloud-Optimized Geotiffs\n",
    "\n",
    "**Part 2. Multiple COGs**\n",
    "\n",
    "This notebook goes over ways to construct a multidimensional xarray DataArray from many 2D COGS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dask\n",
    "import s3fs\n",
    "import intake\n",
    "import os\n",
    "import xarray as xr\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use the same GDAL environment settings as we did for the single COG case\n",
    "env = dict(GDAL_DISABLE_READDIR_ON_OPEN='EMPTY_DIR', \n",
    "           AWS_NO_SIGN_REQUEST='YES',\n",
    "           GDAL_MAX_RAW_BLOCK_CACHE_SIZE='200000000',\n",
    "           GDAL_SWATH_SIZE='200000000',\n",
    "           VSI_CURL_CACHE_SIZE='200000000')\n",
    "os.environ.update(env)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up a connection with credentials and other settings\n",
    "s3 = s3fs.S3FileSystem(anon=True)\n",
    "objects = s3.ls('sentinel-s1-rtc-indigo/tiles/RTC/1/IW/10/T/ET/2020/')\n",
    "images = ['s3://' + obj + '/Gamma0_VV.tif' for obj in objects]\n",
    "print(len(images))\n",
    "images[:6] #january 2020 scenes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GDAL VRT\n",
    "\n",
    "A GDAL VRT file is an XML format that can group together many separate files into separate bands. It's common to create such a file with a the GDAL command line tool `gdalbuildvrt`, illustrated below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#step 1) write a file list that points to the data. GDAL requires special prefixes for this /vsis3/ or /vsicurl/\n",
    "with open('files.txt', 'w') as f:\n",
    "    lines = [x.replace('s3://', '/vsis3/') + '\\n' for x in images[:6]]\n",
    "    f.writelines(lines)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# step 2) create a VRT file\n",
    "!gdalbuildvrt stack.vrt -separate -input_file_list files.txt "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# step 4) open with xarray\n",
    "chunks=dict(band=1, x=2745, y=2745)\n",
    "da = xr.open_rasterio('stack.vrt', chunks=chunks)\n",
    "da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# step 5) optionally modify coordinates (e.g. time dimension extracted from file name)\n",
    "da = da.rename({'band':'time'})\n",
    "da['time'] = [pd.to_datetime(x[60:68]) for x in images[:6]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Recap\n",
    "\n",
    "1. `xr.open_rasterio(stack.vrt)` stores band coordinates as sequential integers (we lose file name and metadata from each individual COG, so it's common to alter the coordinates after opening the dataset)\n",
    "2. data is tied to a reference to a local file ('stack.vrt'), which can cause problems with distributed computing if you don't have access to the local filesystem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## intake-xarray\n",
    "\n",
    "[intake-xarray](https://github.com/intake/intake-xarray) is a plugin for the intake library. It uses fsspec/s3fs under the hood to facilitate loading data into python objects. the function `intake.open_rasterio()` accepts a list of paths. it returns an intake object with a `to_dask()` function that returns an xarray DataArray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# ~ 1s for 6 files\n",
    "\n",
    "# this loads the image ID into xarray's band coordinates. \n",
    "\n",
    "pattern = 's3://sentinel-s1-rtc-indigo/tiles/RTC/1/IW/10/T/ET/2020/{band}/Gamma0_VV.tif'\n",
    "chunks=dict(band=1, x=2745, y=2745)\n",
    "sources = intake.open_rasterio(images[:6], chunks=chunks, path_as_pattern=pattern, concat_dim='band')\n",
    "da = sources.to_dask() \n",
    "da"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### recap:\n",
    "\n",
    "* This is a convient way to avoid constructing a VRT and load a bunch of COGs. It works well as long as the COG urls follow a distinct pattern. Metadata is also lost (we have attributes from the first COG, not others)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Custom\n",
    "\n",
    "You can also just use xarray and dask to construct a larger datacube from many COGS. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# 4 - 8 s\n",
    "# Load all the images\n",
    "\n",
    "chunks=dict(band=1, x=2745, y=2745)\n",
    "dataArrays = [xr.open_rasterio(url, chunks=chunks) for url in images]\n",
    "\n",
    "# note use of join='override' b/c we know these COGS have the same coordinates\n",
    "da = xr.concat(dataArrays, dim='band', join='override', combine_attrs='drop')\n",
    "da = da.rename({'band':'time'})\n",
    "da['time'] = [pd.to_datetime(x[60:68]) for x in images]\n",
    "da"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### recap:\n",
    "\n",
    "* The cell above is essentially a for-loop that iterates over each COG in sequence. 50ms-200ms * 80 ~ 4-16 seconds. The next notebook will look at using Dask to speed things up by opening the files in parallel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize\n",
    "\n",
    "Here is an example of interactive visualization again using hvplot. Since we're using full resolution arrays it's key to set the `rasterize=True` keyword argument. That uses the datashader library to pre-render images before sending them to the browser.\n",
    "\n",
    "This is extremely powerful because, resolution updates as you zoom in, and you can scrub through the data cube with an interactive slider widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hvplot.xarray\n",
    "da.hvplot.image(rasterize=True, aspect='equal', cmap='gray', clim=(0,0.4))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:notebook] *",
   "language": "python",
   "name": "conda-env-notebook-py"
  },
  "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}