1-raster.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with raster datasets\n",
    "\n",
    "questions:\n",
    "  - \"How can I extract pixel values from rasters and perform computations?\"\n",
    "  - \"How might I write pixel values out to a new raster file?\"\n",
    "\n",
    "objectives:\n",
    "  - \"Understand the basic components of a raster dataset and how to access them from a python program.\"\n",
    "  - \"Perform numerical operations on pixel values.\"\n",
    "  - \"Read from and write to raster datasets.\"\n",
    "\n",
    "keypoints:\n",
    "  - \"Rasterio is built around the GDAL library (recall section 3), to facilitate raster operations in Python.\"\n",
    "  - \"Pixel values of rasters can be extracted to a numpy array.\"\n",
    "  - \"Computation is done in local memory on numpy arrays, then saved to the raster format of choice.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Background\n",
    "\n",
    "[GDAL](https://en.wikipedia.org/wiki/GDAL) is a powerful and mature library for reading, writing and warping raster datasets, written in C++ with bindings to other languages. There are a variety of geospatial libraries available on the python package index, and almost all of them depend\n",
    "on GDAL.  One such python library developed and supported by Mapbox, ``rasterio``, builds on top of GDAL's many features, but provides a more pythonic interface and supports many of the features and formats that GDAL supports. Both GDAL and rasterio are constantly being updated and improved: As of writing this tutorial (July 2018), GDAL is at version 2.3.1 and rasterio is at version 1.0.2.\n",
    "\n",
    "**When should you use GDAL directly?**\n",
    "* If you are comfortable with the terminal, [GDAL's command line utilities](https://www.gdal.org/gdal_utilities.html) are very useful for scripting.\n",
    "* Note that GDAL also has auto-generated [Python bindings](https://pypi.org/project/GDAL/), but we recommend using rasterio instead!\n",
    "\n",
    "**When should you use ``rasterio`` instead of GDAL?**\n",
    "* Maybe always?! Rasterio also has a set of command line tools \n",
    "* If you are working in a Python environment (ipython, scripts, jupyter lab)\n",
    "\n",
    "**When might these not be the best tools?**\n",
    "* Both libraries are critical for input/output operations, but you'll draw on other libraries for computation (e.g. numpy)\n",
    "    * That said, GDAL does have some standard processing scripts (for example [pan-sharpening](https://www.gdal.org/gdal_pansharpen.html)) and rasterio provides a [plugin interface](https://github.com/mapbox/rasterio/wiki/Rio-plugin-registry) for workflows\n",
    "* For polished map creation and interactive visualization, a desktop GIS software like [QGIS](https://qgis.org/en/site/) may be a better, more fully-featured choice.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Loading and viewing Landsat 8 imagery\n",
    "\n",
    "The Landsat program is the longest-running civilian satellite imagery program, with the first satellite launched in 1972 by the US Geological Survey. Landsat 8 is the latest satellite in this program, and was launched in 2013. Landsat observations are processed into “scenes”, each of which is approximately 183 km x 170 km, with a spatial resolution of 30 meters and a temporal resolution of 16 days. The duration of the landsat program makes it an attractive source of medium-scale imagery for land surface change analyses. Landsat is multiband imagery, you can read more about it [here](https://landsat.usgs.gov).\n",
    "\n",
    "One reason we'll use Landsat 8 for this demo is that the entire Landsat 8 archive is hosted by various commercial Cloud providers with free public access ([AWS](https://registry.opendata.aws/landsat-8/) and [Google Cloud](https://cloud.google.com/storage/docs/public-datasets/landsat))! \n",
    "\n",
    "We start by importing all the python libraries we need in this tutorial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "import rasterio.plot\n",
    "import pyproj\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print('rasterio version: ', rasterio.__version__)\n",
    "print('pyproj version: ', pyproj.__version__)\n",
    "print('matplotlib version: ', matplotlib.__version__)\n",
    "print('numpy version: ', np.__version__)\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Satellite archives on the Cloud\n",
    "\n",
    "A single Landsat8 scene is about 1 Gb in size since it contains a large array of data for each imagery band. Images are often grouped into a single .tar.gz file, which requires scientists to download the entire 1 Gb, even if you are just interested in a single band or subset of the image. To circumvent downloading, a new format has been proposed for Cloud storage called [Cloud-optimized Geotiffs (COGs)](http://www.cogeo.org), to allow easy access to individual bands or even subsets of single images. \n",
    "\n",
    "In a strict sense, the Landsat images on AWS and Google Cloud are not COGs (you can validate the format with online tools such as this - http://cog-validate.radiant.earth/html). Nevertheless, the files are stored as tiled Geotiffs, so many COG features still work.\n",
    "\n",
    "\n",
    "### AWS versus Google Cloud\n",
    "You can use [NASA's Earthdata Search](https://search.earthdata.nasa.gov) website to discover data. Landsat images are organized by 'path' and 'row'. We've chosen a scenes from path 42, row 34, that doesn't have many clouds present (LC08_L1TP_042034_20170616_20170629_01_T1). Note that 'T1' stands for 'Tier 1' (for analytic use), and 'RT' stands for 'Real-time', for which quality control is not as rigorous. Read more about the various Landsat formats and collections [here](https://landsat.usgs.gov/landsat-collections).\n",
    "\n",
    "At first glance it seems that you can find this data on both AWS or Google Cloud, for example look at the band 4 image for the date we selected:\n",
    "\n",
    "[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](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)\n",
    "\n",
    "\n",
    "Or, the same file is available on Google Cloud storage:\n",
    "\n",
    "[https://storage.googleapis.com/gcp-public-data-landsat/LC08/01/042/034/LC08_L1TP_042034_20170616_20170629_01_T1/LC08_L1TP_042034_20170616_20170629_01_T1_B4.TIF](https://storage.googleapis.com/gcp-public-data-landsat/LC08/01/042/034/LC08_L1TP_042034_20170616_20170629_01_T1/LC08_L1TP_042034_20170616_20170629_01_T1_B4.TIF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Landsat on Google:')\n",
    "filepath = 'https://storage.googleapis.com/gcp-public-data-landsat/LC08/01/042/034/LC08_L1TP_042034_20170616_20170629_01_T1/LC08_L1TP_042034_20170616_20170629_01_T1_B4.TIF'\n",
    "with rasterio.open(filepath) as src:\n",
    "    print(src.profile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Landsat on AWS:')\n",
    "filepath = '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'\n",
    "with rasterio.open(filepath) as src:\n",
    "    print(src.profile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What just happend?\n",
    "\n",
    "If you’re familiary with programming in python, you’ve probably seen context managers before. This context manager, **rasterio.open()** functions like the python standard library function open for opening files. The block of code within the with ... as statement is executed once the file is opened, and the file is closed when the context manager exits. This means that we don’t have to manually close the raster file, as the context manager handles that for us.\n",
    "\n",
    "Instead of a local file path, rasterio knows how to read URLs too, so we just passed the link to the file on AWS\n",
    "\n",
    "**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](http://www.spatialreference.org/ref/epsg/wgs-84-utm-zone-11n), has a simple affine transformation, is chunked into smaller 512x512 arrays, tiled and compressed on the AWS hard drive where it is stored.\n",
    "\n",
    "**Note that we have not actually downloaded the image! We just read the metadata.** \n",
    "\n",
    "Did you catch the subtle differences between the data on AWS versus the data on Google? The image block sizes, or \"tiling\" is different, as is the compression scheme. It also turns out that AWS stores pre-made image overviews, but Google does not. It is important to realize that not all archives are the same! In fact the Google Cloud storage archive is a complete mirror of the full USGS archive, whereas as of writing this, AWS only has collection 1 scenes after 2017. In general, performance will be best if you are accessing data locally (so if you are running this notebook on AWS, use the landsat public data on AWS)! We are going to use the AWS archive for the rest of this tutorial because of the convenience of overviews:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot a low-resolution overview"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The grid of raster values can be accessed as a numpy array and plotted:\n",
    "with rasterio.open(filepath) as src:\n",
    "    oviews = src.overviews(1) # list of overviews from biggest to smallest\n",
    "    oview = oviews[-1] # let's look at the smallest thumbnail\n",
    "    print('Decimation factor= {}'.format(oview))\n",
    "    # NOTE this is using a 'decimated read' (http://rasterio.readthedocs.io/en/latest/topics/resampling.html)\n",
    "    thumbnail = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))\n",
    "    \n",
    "print('array type: ',type(thumbnail))\n",
    "print(thumbnail)\n",
    "    \n",
    "plt.imshow(thumbnail)\n",
    "plt.colorbar()\n",
    "plt.title('Overview - Band 4 {}'.format(thumbnail.shape))\n",
    "plt.xlabel('Column #')\n",
    "plt.ylabel('Row #')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Earlier we saw in the metadata that a no-data value wasn't set, but pixels outside the imaged area are clearly set to \"0\". Images commonly look like this because of satellite orbits and the fact that the Earth is rotating as imagery is acquired! The colormap is often improved if we change the out of bounds area to NaN. To do this we have to convert the datatype from uint16 to float32 (so be aware the array with NaNs will take 2x the storage space). The serrated edge is due to the coarse sampling of the full resolution image that we are doing here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(filepath) as src:\n",
    "    oviews = src.overviews(1) \n",
    "    oview = oviews[-1]\n",
    "    print('Decimation factor= {}'.format(oview))\n",
    "    thumbnail = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))\n",
    "    \n",
    "    thumbnail = thumbnail.astype('f4')\n",
    "    thumbnail[thumbnail==0] = np.nan\n",
    "    \n",
    "plt.imshow(thumbnail)\n",
    "plt.colorbar()\n",
    "plt.title('Overview - Band 4 {}'.format(thumbnail.shape))\n",
    "plt.xlabel('Column #')\n",
    "plt.ylabel('Row #')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pull an image subset at full resolution\n",
    "\n",
    "A nice feature of COGs is that you can request a subset of the image and only that subset will be downloaded to your computer's memory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#https://rasterio.readthedocs.io/en/latest/topics/windowed-rw.html\n",
    "\n",
    "#rasterio.windows.Window(col_off, row_off, width, height)\n",
    "window = rasterio.windows.Window(1024, 1024, 1280, 2560)\n",
    "\n",
    "with rasterio.open(filepath) as src:\n",
    "    subset = src.read(1, window=window)\n",
    "\n",
    "plt.figure(figsize=(6,8.5))\n",
    "plt.imshow(subset)\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.title(f'Band 4 Subset\\n{window}')\n",
    "plt.xlabel('Column #')\n",
    "plt.ylabel('Row #')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Example computation: NDVI\n",
    "\n",
    "The [Normalized Difference Vegetation\n",
    "Index](https://en.wikipedia.org/wiki/Normalized_Difference_Vegetation_Index) is a simple indicator that can be used to assess whether the target includes healthy vegetation.   This calculation uses two bands of a multispectral image dataset, the Red and Near-Infrared (NIR) bands:\n",
    "\n",
    "$$NDVI = \\frac{(NIR - Red)}{(NIR + Red)}$$\n",
    "\n",
    "For this tutorial, we'll use the NIR and Red bands from a Landsat-8 scene above part of the central valley and the Sierra Nevada in California.  We'll be using [Level 1TP datasets](https://landsat.usgs.gov/landsat-processing-details), orthorectified, map-projected images containing radiometrically calibrated data. \n",
    "\n",
    "### Bands\n",
    "\n",
    "* Red: Band 4 \n",
    "* Near-Infrared: Band 5\n",
    "\n",
    "Because of the longevity of the landsat mission and because different sensors on the satellite record data at different resolutions, these bands are individually stored as single-band raster files.  Some other rasters may store multiple bands in the same file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the same example image:\n",
    "date = '2017-06-16'\n",
    "url = 'http://landsat-pds.s3.amazonaws.com/c1/L8/042/034/LC08_L1TP_042034_20170616_20170629_01_T1/'\n",
    "redband = 'LC08_L1TP_042034_20170616_20170629_01_T1_B{}.TIF'.format(4)\n",
    "nirband = 'LC08_L1TP_042034_20170616_20170629_01_T1_B{}.TIF'.format(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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)))\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": [
    "with rasterio.open(url+nirband) as src:\n",
    "    oviews = src.overviews(1) # list of overviews from biggest to smallest\n",
    "    oview = oviews[1]  # Use second-highest resolution overview\n",
    "    nir = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))\n",
    "    \n",
    "plt.imshow(nir)\n",
    "plt.colorbar()\n",
    "plt.title('{}\\nNIR {}'.format(nirband, nir.shape))\n",
    "plt.xlabel('Column #')\n",
    "plt.ylabel('Row #')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_ndvi(nir,red):\n",
    "    '''Calculate NDVI from integer arrays'''\n",
    "    nir = nir.astype('f4')\n",
    "    red = red.astype('f4')\n",
    "    ndvi = (nir - red) / (nir + red)\n",
    "    return ndvi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi = calc_ndvi(nir,red)\n",
    "plt.imshow(ndvi, cmap='RdYlGn')\n",
    "plt.colorbar()\n",
    "plt.title('NDVI {}'.format(date))\n",
    "plt.xlabel('Column #')\n",
    "plt.ylabel('Row #')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Save the NDVI raster to local disk\n",
    "So far, we have read in a cloud-optimized geotiff from the Cloud into our computer memory (RAM), and done a simple computation. What if we want to save this result locally for future use?\n",
    "\n",
    "Since we have used a subsampled overview, we have to modify the orginal metadata before saving! In particular, the Affine matrix describing the coordinates is different (new resolution and extents, and we've changed the datatype). We'll stick with Geotiff, but note that it is no longer cloud-optimized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "localname = 'LC08_L1TP_042034_20170616_20170629_01_T1_NDVI_OVIEW.tif'\n",
    "\n",
    "with rasterio.open(url+nirband) as src:\n",
    "    profile = src.profile.copy()\n",
    "    \n",
    "    aff = src.transform\n",
    "    newaff = rasterio.Affine(aff.a * oview, aff.b, aff.c,\n",
    "                             aff.d, aff.e * oview, aff.f)\n",
    "    profile.update({\n",
    "            'dtype': 'float32',\n",
    "            'height': ndvi.shape[0],\n",
    "            'width': ndvi.shape[1],\n",
    "            'transform': newaff})  \n",
    "    \n",
    "    with rasterio.open(localname, 'w', **profile) as dst:\n",
    "        dst.write_band(1, ndvi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reopen the file and plot\n",
    "with rasterio.open(localname) as src:\n",
    "    print(src.profile)\n",
    "    ndvi = src.read(1) # read the entire array\n",
    "\n",
    "plt.imshow(ndvi, cmap='RdYlGn')\n",
    "plt.colorbar()\n",
    "plt.title('NDVI {}'.format(date))\n",
    "plt.xlabel('Column #')\n",
    "plt.ylabel('Row #')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rasterio also has a 'convenience method' for plotting with georeferenced coordinates\n",
    "# in this case, coordinates are Easting [m] and Northing [m], and colorbar is default instead of RdYlGn\n",
    "with rasterio.open(localname) as src:\n",
    "    fig, ax = plt.subplots()\n",
    "    rasterio.plot.show(src, ax=ax, title='NDVI') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatial indexing and extracting values\n",
    "Raster images really have two sets of coordinates. First, 'image coordinates' correspond to the row and column for a specific pixel. Second, the 'spatial coordinates' correspond to the location of each pixel on the surface of the Earth. Rasterio makes it convenient to use both coordinate systems.\n",
    "\n",
    "Lets say you want the value of NDVI at a specific point in this scene. For example Fresno, CA (-119.770163586, 36.741997032). But the image is in UTM coordinates, so you have to first convert these points to UTM (Many websites will do this for you https://www.geoplaner.com), or you can use the pyproj library. Be warned that this only works if you read the dataset at it's full resolution (if you do decimated or windowed reads of the file, you need to adjust the affine transform)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(localname) as src:\n",
    "    # Use pyproj to convert point coordinates\n",
    "    utm = pyproj.Proj(src.crs) # Pass CRS of image from rasterio\n",
    "    lonlat = pyproj.Proj(init='epsg:4326')\n",
    "\n",
    "    lon,lat = (-119.770163586, 36.741997032)\n",
    "    east,north = pyproj.transform(lonlat, utm, lon, lat)\n",
    "\n",
    "    print('Fresno NDVI\\n-------')\n",
    "    print(f'lon,lat=\\t\\t({lon:.2f},{lat:.2f})')\n",
    "    print(f'easting,northing=\\t({east:g},{north:g})')\n",
    "\n",
    "    # What is the corresponding row and column in our image?\n",
    "    row, col = src.index(east, north) # spatial --> image coordinates\n",
    "    print(f'row,col=\\t\\t({row},{col})')\n",
    "\n",
    "    # What is the NDVI?\n",
    "    value = ndvi[row, col]\n",
    "    print(f'ndvi=\\t\\t\\t{value:.2f}')\n",
    "\n",
    "\n",
    "    # Or if you see an interesting feature and want to know the spatial coordinates:\n",
    "    row, col = 200, 450\n",
    "    east, north = src.xy(row,col) # image --> spatial coordinates\n",
    "    lon,lat = pyproj.transform(utm, lonlat, east, north)\n",
    "    value = ndvi[row, col]\n",
    "    print(f'''\n",
    "Interesting Feature\n",
    "-------\n",
    "row,col=          ({row},{col})\n",
    "easting,northing= ({east:g},{north:g})\n",
    "lon,lat=          ({lon:.2f},{lat:.2f})\n",
    "ndvi=              {value:.2f}\n",
    "''')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Calculate change in NDVI over time\n",
    "\n",
    "Let's take a look at the difference in NDVI between a scene in June 2013 and June 2017. If you went to the [AWS Landsat Archive page](https://registry.opendata.aws/landsat-8), you probably noticed that it isn't obvious how to search and discover images (most of the time you probably won't know the row, path, or full URL of images over your area of interest!) There are many options for searching, graphical web applications like [NASA's Earthdata Search](https://search.earthdata.nasa.gov), or convenient Python tools like DevSeed's [sat-search](https://github.com/sat-utils/sat-search), among many others! We won't go into these tools here, but we encourage you to experiment with your own image scenes. Here is a file from June 2018 that was easy to find with Landsat for AWS using the Path and Row from the earlier URL:\n",
    "\n",
    "https://landsatonaws.com/L8/042/034/LC08_L1TP_042034_20180619_20180703_01_T1\n",
    "\n",
    "With the band 4 imagery URL:\n",
    "\n",
    "http://landsat-pds.s3.amazonaws.com/c1/L8/042/034/LC08_L1TP_042034_20180619_20180703_01_T1/LC08_L1TP_042034_20180619_20180703_01_T1_B4.TIF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the same example image:\n",
    "date2 = '2018-06-19'\n",
    "url2 = 'http://landsat-pds.s3.amazonaws.com/c1/L8/042/034/LC08_L1TP_042034_20180619_20180703_01_T1/'\n",
    "redband2 = 'LC08_L1TP_042034_20180619_20180703_01_T1_B{}.TIF'.format(4)\n",
    "nirband2 = 'LC08_L1TP_042034_20180619_20180703_01_T1_B{}.TIF'.format(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filepath = url2+redband2\n",
    "with rasterio.open(filepath) as src:\n",
    "    print('Opening:', filepath)\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",
    "    red2 = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))\n",
    "\n",
    "filepath = url2+nirband2\n",
    "with rasterio.open(filepath) as src:\n",
    "    print('Opening:', filepath)\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",
    "    nir2 = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))\n",
    "\n",
    "ndvi2 = calc_ndvi(nir2, red2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,3, figsize=(14,6), sharex=True, sharey=True)\n",
    "\n",
    "plt.sca(axes[0])\n",
    "plt.imshow(ndvi, cmap='RdYlGn', vmin=-1, vmax=1)\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.title('NDVI {}'.format(date))\n",
    "plt.xlabel('Column #')\n",
    "plt.ylabel('Row #')\n",
    "\n",
    "plt.sca(axes[1])\n",
    "plt.imshow(ndvi2, cmap='RdYlGn', vmin=-1, vmax=1)\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.title('NDVI {}'.format(date2))\n",
    "\n",
    "plt.sca(axes[2])\n",
    "plt.imshow(ndvi2 - ndvi, cmap='bwr', vmin=-1, vmax=1)\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.title('Diff ({} - {})'.format(date2, date))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What just happened?\n",
    "\n",
    "We just loaded 4 decimated Landsat 8 band images into memory and computed the difference in NDVI between two dates. That was relatively easy because these two rasters happen to use the same coordinate system and grid. This is know as 'Analysis Ready Data'. In the geosciences, we commonly have data in WGS84 Lat, Lon. Or what if you want to compare to NDVI estimated from Sentinel-2 satellite acquisitions, which are on a different grid?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Advanced uses of rasterio \n",
    "\n",
    "Rasterio can also be used for masking, reprojecting, and regridding distinct datasets. Here is one simple example to reproject our local NDVI onto a WGS84 Lat/Lon Grid. The example creates a 'VRT' file, which is merely a ASCII text file describing the transformation and other image metadata. The array values do no need to be duplicated, the VRT file just contains a reference to the local file (LC08_L1TP_042034_20170616_20170629_01_T1_NDVI_OVIEW.tif). This is extremely useful to avoid duplicating lots of data and filling up your computer!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio.warp\n",
    "import rasterio.shutil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "localname = 'LC08_L1TP_042034_20170616_20170629_01_T1_NDVI_OVIEW.tif'\n",
    "vrtname = 'LC08_L1TP_042034_20170616_20170629_01_T1_NDVI_OVIEW_WGS84.vrt'\n",
    "\n",
    "with rasterio.open(localname) as src:\n",
    "    with rasterio.vrt.WarpedVRT(src, crs='epsg:4326', resampling=rasterio.enums.Resampling.bilinear) as vrt:\n",
    "        rasterio.shutil.copy(vrt, vrtname, driver='VRT')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE our coordinates have changed to lat, lon. we should probably crop the edge artifacts do to reprojection too!\n",
    "with rasterio.open(vrtname) as src:\n",
    "    rasterio.plot.show(src, title='NDVI', cmap='RdYlGn', vmin=-1, vmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The issue with VRTs is now you have two files and let's say you want to share this image with a colleague. Also if you move things around on your computer, the paths might get mixed up, so in some cases it's nice to save a complete reprojected file. The second code block does this, saving our NDVI as a Geotiff in WGS84."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "localname = 'LC08_L1TP_042034_20170616_20170629_01_T1_NDVI_OVIEW.tif'\n",
    "tifname = 'LC08_L1TP_042034_20170616_20170629_01_T1_NDVI_OVIEW_WGS84.tif'\n",
    "\n",
    "dst_crs = 'EPSG:4326'\n",
    "\n",
    "with rasterio.open(localname) as src:\n",
    "    profile = src.profile.copy()\n",
    "    \n",
    "    transform, width, height = rasterio.warp.calculate_default_transform(\n",
    "        src.crs, dst_crs, src.width, src.height, *src.bounds)\n",
    "    \n",
    "    profile.update({\n",
    "        'crs': dst_crs,\n",
    "        'transform': transform,\n",
    "        'width': width,\n",
    "        'height': height\n",
    "    })\n",
    "    \n",
    "    with rasterio.open(tifname, 'w', **profile) as dst:\n",
    "        rasterio.warp.reproject(\n",
    "            source=rasterio.band(src, 1),\n",
    "            destination=rasterio.band(dst, 1),\n",
    "            src_transform=src.transform,\n",
    "            src_crs=src.crs,\n",
    "            dst_transform=transform,\n",
    "            dst_crs=dst_crs,\n",
    "            resampling=rasterio.warp.Resampling.bilinear)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE our coordinates have changed to lat, lon. we should probably crop the edge artifacts do to reprojection too!\n",
    "with rasterio.open(tifname) as src:\n",
    "    rasterio.plot.show(src, title='NDVI', cmap='RdYlGn', vmin=-1, vmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Additional resources\n",
    "\n",
    "This just scratches the surface of what you can do with rasterio. Check out the excellent official documentation as you continue exploring:\n",
    "\n",
    "https://rasterio.readthedocs.io/en/latest/"
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}