gridded_data_4_rasterio.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Waterhackweek cyberseminar\n",
    "# Workflows for gridded climate datasets\n",
    "## Diana Gergel and Bart Nijssen\n",
    "### University of Washington\n",
    "### February 14, 2019"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Abstract:** Climate change, forecasting, satellite datasets, large model ensembles ... Large gridded datasets are everywhere in hydrology and earth science. While accessing and analyzing these datasets required some serious programming skills not so long ago, a number of toolkits are now available that let you easily access, ingest, analyze and display gridded climate datasets. In this webinar we’ll discuss one of the most common file formats used in our field for large data sets, the Network Common Data Format (NetCDF), and step through a Jupyter notebook to showcase python packages, such as xarray and cartopy, that can be used to examine them. No prior experience required, although we will build on some of the skills you have acquired in earlier webinars in the series."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GitHub repository\n",
    "\n",
    "The GitHub repository for this presentation can be found at: https://github.com/bartnijssen/gridded_data\n",
    "\n",
    "The README.md file in this repository explains how to configure the environment you need to run the sample notebooks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline\n",
    "\n",
    "* `gridded_data_1_background.ipynb`\n",
    "\n",
    "\n",
    "* `gridded_data_2_xarray_and_dask.ipynb`\n",
    "\n",
    "\n",
    "* `gridded_data_3_cartopy_and_seaborn.ipynb`\n",
    "\n",
    "\n",
    "* `gridded_data_4_rasterio.ipynb` (this notebook)\n",
    "    + Opening and displaying a GeoTIFF using rasterio\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other data formats\n",
    "\n",
    "Clearly there are many standardized data formats other than NetCDF, even though that is the file format we tend to use most commonly in our line of research. Fortunately the community contributing to python is large and you will be able to find tools to access most common (and many uncommon) data formats. If you can't find one, maybe you can provide such a package to the rest of the python community or perhaps you should seriously consider using a more widely shared data format to store your data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## rasterio\n",
    "\n",
    "In this last notebook we will very briefly show one other python package commonly used to access data files that are commonly used in the GIS and geospatial communities. These are data files that also include some metadata, in particular projection information that maps each data point in the file to a geospatial location.\n",
    "\n",
    "As a caveat - I hardly ever work with these files and for the following commands I borrow heavily from the rasterio documentation at https://rasterio.readthedocs.io/ and from an earlier geohackweek tutorial on working with raster datasets at https://geohackweek.github.io/raster/04-workingwithrasters/. Please look there if you want more information.\n",
    "\n",
    "Here, we will limit ourselves at opening, querying and displaying some geotiff files and showcasing a few features of rasterio.\n",
    "\n",
    "rasterio is built on [GDAL](https://gdal.org/), the Geospatial Data Abstraction Library, which is a translator library for raster and vector geospatial data formats. There are command-line utilities that are built on this library as well as python bindings, but they can be difficult to use. rasterio's goal is to bring all of GDAL's raster functionality cleanly into python, while performing as fast as GDAL’s Python bindings. It reads raster files and returns them as Numpy arrays, which can then be manipulated by other python code. Note that rasterio also has a set of command-line utilities. \n",
    "\n",
    "Basically, rasterio will make it easy to read and write geo-registered files. Computations that are done on the data you read can be done with other python libraries such as numpy or even xarray as we will see at the end of this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Opening a dataset\n",
    "\n",
    "In this example, we will use rasterio to read some Geotiff files. Geotiff is a common format in GIS applications, because it allows georeferencing information to be embedded within a TIFF file. For example, orthorectified aerial photos and satellite images are often distributed in GeoTIFF format. \n",
    "\n",
    "Here we will take advantage of rasterio to directly read data from a URL. \n",
    "\n",
    "First we import rasterio:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and then define the URL we want to read. In this case, we will read from a landsat image stored on AWS (note that this is the image used in the [geohackweek tutorial](https://geohackweek.github.io/raster/04-workingwithrasters/)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_url = 'http://landsat-pds.s3.amazonaws.com/c1/L8/042/034/LC08_L1TP_042034_20170616_20170629_01_T1/LC08_L1TP_042034_20170616_20170629_01_T1_B4.TIF'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a landsat-8 image from path 42, row 34, taken in June 2017. 'T1' stands for 'Tier 1' (for analytic use) and 'B4' stands for 'Band 4', which indicates a specific wavelength band."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use raterio to read the image metadata, without downloading the entire image (remote sensing images can be large):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(image_url) as src:\n",
    "    print(src.profile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a bit hard on the eyes, so let's format it a bit better"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(image_url) as src:\n",
    "    for key, value in src.profile.items():\n",
    "        print(\"{} : {}\".format(key, value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So `src.profile` is a collection of metadata for the file. We see that it is a Geotiff (Gtiff), the image values are unsigned integer format, nodata values are not assigned, the image has a dimensions of 7711x7531, is a single band, is in UTM coordinates, has a simple affine transformation, is chunked into smaller 512x512 arrays, tiled and compressed on the AWS hard drive where it is stored."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an aside: The `with` construction is commonly used in python (see [here](https://docs.python.org/3/reference/compound_stmts.html#the-with-statement)) for more detail. The advantage of using it with file access is that the file is automatically closed at the end of the block, even if the block fails. So when you do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(image_url) as src:\n",
    "    print('In block\\t: {}'.format(src))\n",
    "print('After block\\t: {}'.format(src))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You see that `src` is an `open DatasetReader` within the `with` block and a `closed DatasetReader` after the block, even though we never closed the reader explicitly. That is often very useful, but not so much in this case, so we will just open `image_url` without the `with` block:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src = rasterio.open(image_url)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then quickly see what all the elements are in `src`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir(src)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To read the actual data from `src`, we would use `src.read(n)`, where `n` is the band we want to read. Since this dataset has only one band (Band 4), we would use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = src.read(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can inspect the data type of `data`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and see that this is now a `numpy ndarray` or n-dimensional array. You can now do computations on this like on any other numpy array. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get a quick overview of the content of the file itself, we can use rasterio's plot functionality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio.plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rasterio.plot.show(src)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the x and y-labels in this case are actually georeferenced coordinates, rather than just image and line numbers. If you would plot the data contained in the file using `matplotlib`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "you get the same image, but the georeferencing is lost. We'll leave it as an exercise for you how to plot this data using `cartopy`. Take a look at the [geohackweek tutorial](https://geohackweek.github.io/raster/04-workingwithrasters/) to see how you can subset the image, generate thumbnails, or combine multiple bands to calculate the normalized vegetation difference index (NDVI)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## xarray and rasterio\n",
    "\n",
    "Lastly, in python there are often many ways to do the same thing. For example, having seen some of the `xarray` functionality, you may wonder how you can convert the data you read with rasterio into `xarray`. An experimental interface exists to read files that rasterio can read directly into `xarray` as a `DataArray`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da = xr.open_rasterio(image_url, parse_coordinates=True, chunks=None, cache=None, lock=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da[0].plot.imshow()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "whw_gridded_data",
   "language": "python",
   "name": "whw_gridded_data"
  },
  "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
}