add mosaicing exercise & key for SciPy class

Author: Sara Safavi

jupyter-notebooks/in-class-exercises/mosaicing-and-masking/data/mt-dana-small.geojson

@@ -0,0 +1 @@
+{"type":"FeatureCollection","features":[{"type":"Feature","properties":{},"geometry":{"type":"Polygon","coordinates":[[[-119.16183471679688,37.82903964181452],[-119.14947509765626,37.83663205340172],[-119.13745880126953,37.846392577323286],[-119.13574218750001,37.856422880849514],[-119.13883209228514,37.86645181975611],[-119.12406921386719,37.86916210952103],[-119.12200927734376,37.875937397778614],[-119.1212688230194,37.90572368618133],[-119.13740499245301,37.930641295117404],[-119.16595458984376,37.92659678938742],[-119.18243408203126,37.9447389942697],[-119.2088161252655,37.95257263611974],[-119.25516469704283,37.92522514171301],[-119.2630611203827,37.88215253011582],[-119.25104482399598,37.84474832157969],[-119.18203695046083,37.82576791597315],[-119.16183471679688,37.82903964181452]]]}}]}

jupyter-notebooks/in-class-exercises/mosaicing-and-masking/images/explorer-data-order.png

@@ -0,0 +1,345 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Creating a composite image from multiple PlanetScope scenes\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this exercise, you'll learn how to create a composite image (or mosaic) from multiple PlanetScope satellite images that cover an area of interest (AOI). We'll use `rasterio`, along with its vector-data counterpart `fiona`, to do this. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step 1. Aquiring Imagery"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In order to visually search for imagery in our AOI, we'll use [Planet Explorer](https://www.planet.com/explorer/).\n",
+    "\n",
+    "For this exercise, we're going to visit Yosemite National Park. In the screenshot below you'll see an AOI drawn around [Mount Dana](https://en.wikipedia.org/wiki/Mount_Dana) on the eastern border of Yosemite. You can use [data/mt-dana-small.geojson](data/mt-dana-small.geojson) to search for this same AOI yourself.\n",
+    "\n",
+    "Here we want an image that depicts the mountain on a clear summer day, so for this data search in Planet Explorer we'll set the filters to show only scenes with less than 5% cloud cover, and narrow down the date range to images captured between July 1-July 31, 2017. Since we're only interested in PlanetScope data, and we're creating a visual - not analytic - product, we'll set the Source to `3-band PlanetScope Scene`. Finally, since we want to create a mosaic that includes our entire AOI, we'll set the Area coverage to full coverage."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "![Mount Dana in Planet Explorer](images/explorer-mt-dana.gif)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As you can see in the animated gif above, this search yields multiple days within July 2017 that match our filters. After previewing a few days, I decided I like the look of July 21, 2017.\n",
+    "\n",
+    "After selecting a single day, you can roll over the individual images to preview their coverage. In the gif above, you'll notice that it takes three individual images to completely cover our AOI. In this instance, as I roll over each item in Planet Explorer I can see that the scenes' rectangular footprints extend far beyond Mount Dana. All three scenes overlap slightly, and one scene touches only a small section at the bottom of the AOI. Still, taken together the images provide 100% coverage, so we'll go ahead and place an order for the Visual imagery products for these three scenes."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "![Download imagery Planet Explorer](images/explorer-data-order.png)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Once the order is ready, download the images, extract them from the .zip, and move them into the `data/` directory adjacent to this Notebook."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step 2. Inspecting Imagery"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Load our 3 images using rasterio\n",
+    "\n",
+    "import rasterio\n",
+    "\n",
+    "img1 = rasterio.open('data/20170721_175836_103c/20170721_175836_103c_3B_Visual.tif')\n",
+    "img2 = rasterio.open('data/20170721_175837_103c/20170721_175837_103c_3B_Visual.tif')\n",
+    "img3 = rasterio.open('data/20170721_175838_103c/20170721_175838_103c_3B_Visual.tif')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "At this point we can use `rasterio` to inspect the metadata of these three images. Specifically, in order to create a composite from these images, we want to verify that all three images have the same data type, the same coordinate reference systems and the same band count:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(img1.meta['dtype'], img1.meta['crs'], img1.meta['count'])\n",
+    "print(img2.meta['dtype'], img2.meta['crs'], img2.meta['count'])\n",
+    "print(img3.meta['dtype'], img3.meta['crs'], img3.meta['count'])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Success - they do! But wait, I thought we were using a \"Visual\" image, and expecting only 3 bands of information (RGB)?\n",
+    "\n",
+    "Let's take a closer look at what these bands contain:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Read in color interpretations of each band in img1 - here we'll assume img2 and img3 have the same values\n",
+    "colors = [img1.colorinterp[band] for band in range(img1.count)]\n",
+    "\n",
+    "# take a look at img1's band types:\n",
+    "for color in colors:\n",
+    "    print(color.name)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The fourth channel is actually a binary alpha mask: this is common in satellite color models, and can be confirmed in Planet's [documentation on the PSSCene3Band product](https://developers.planet.com/docs/api/psscene3band/).\n",
+    "\n",
+    "Now that we've verified all three satellite images have the same critical metadata, we can safely use `rasterio.merge` to stitch them together."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step 3. Creating the Mosaic"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from rasterio.merge import merge\n",
+    "\n",
+    "# merge returns the mosaic & coordinate transformation information\n",
+    "(mosaic, transform) = merge([img1, img2, img3])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Once that process is complete, take a moment to congratulate yourself. At this stage you've successfully acquired adjacent imagery, inspected metadata, and performed a compositing process in order to generate a new mosaic. Well done!\n",
+    "\n",
+    "Before we go further, let's use `rasterio.plot` (a matplotlib interface) to preview the results of our mosaic. This will just give us a quick-and-dirty visual representation of the results, but it can be useful to verify the compositing did what we expected."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from rasterio.plot import show\n",
+    "\n",
+    "show(mosaic)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "At this point we're ready to write our mosaic out to a new GeoTIFF file. To do this, we'll want to grab the geospatial metadata from one of our original images (again, here we'll use img1 to represent the metadata of all 3 input images)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Grab a copy of our source metadata, using img1\n",
+    "meta = img1.meta\n",
+    "\n",
+    "# Update the original metadata to reflect the specifics of our new mosaic\n",
+    "meta.update({\"transform\": transform,\n",
+    "    \"height\":mosaic.shape[1],\n",
+    "    \"width\":mosaic.shape[2]})\n",
+    "\n",
+    "with rasterio.open('data/mosaic.tif', 'w', **meta) as dst:\n",
+    "    dst.write(mosaic)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step 4. Clip the Mosaic to AOI Boundaries"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now that we've successfully created a composite mosaic of three input images, the final step is to clip that mosaic to our area of interest. To do that, we'll create a mask for our mosaic based on the AOI boundaries, and crop the mosaic to the extents of that mask.\n",
+    "\n",
+    "You'll recall that we used Explorer to search for Mount Dana, in Yosemite National Park. The GeoJSON file we used for that search can also be used here, to provide a mask outline for our mosaic.\n",
+    "\n",
+    "For this step we're going to do a couple things: first, we'll use rasterio's sister-library `fiona` to read in the GeoJSON file. Just as `rasterio` is used to manipulate raster data, `fiona` works similarly on vector data. Where `rasterio` represents raster imagery as numpy arrays, `fiona` represents vector data as GeoJSON-like Python dicts. You can learn [more about fiona here](http://toblerity.org/fiona/manual.html).\n",
+    "\n",
+    "After reading in the GeoJSON you'll want to extract the _geometry_ of the AOI (_hint:_ `geometry` will be the dict key). \n",
+    "\n",
+    "#### A note about Coordinate Reference Systems\n",
+    "\n",
+    "If you attempt to apply the AOI to the mosaic imagery now, rasterio will throw errors, telling you that the two datasets do not overlap. This is because the Coordinate Reference System (CRS) used by each dataset do not match. You can verify this by reading the `crs` attribute of the Collection object generated by `fiona.open()`.\n",
+    "\n",
+    "_Hint: the CRS of mt-dana-small.geojson is:_  `'epsg:4326'`\n",
+    "\n",
+    "You'll recall that earlier we validated the metadata of the original input imagery, and learned it had a CRS of `'epsg:32611'`. \n",
+    "\n",
+    "Before the clip can be applied, you will need to to transform the geometry of the AOI to match the CRS of the imagery. Luckily, fiona is smart enough to apply the necessary mathematical transformation to a set of coordinates in order to convert them to new values: apply `fiona.transform.transform_geom` to the AOI geometry to do this, specifying the GeoJSON's CRS as the source CRS, and the imagery's CRS as the destination CRS."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# use rasterio's sister-library for working with vector data\n",
+    "import fiona\n",
+    "\n",
+    "# use fiona to open our original AOI GeoJSON\n",
+    "with fiona.open('data/mt-dana-small.geojson') as mt:\n",
+    "    aoi = [feature[\"geometry\"] for feature in mt]\n",
+    "    \n",
+    "# transform AOI to match mosaic CRS\n",
+    "from fiona.transform import transform_geom\n",
+    "transformed_coords = transform_geom('EPSG:4326', 'EPSG:32611', aoi[0])\n",
+    "\n",
+    "aoi = [transformed_coords]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "At this stage you have read in the AOI geometry and transformed its coordinates to match the mosaic. We're now ready to use `rasterio.mask.mask` to create a mask over our mosaic, using the AOI geometry as the mask line. \n",
+    "\n",
+    "Passing `crop=True` to the mask function will automatically crop the bits of our mosaic that fall outside the mask boundary: you can think of it as applying the AOI as a cookie cutter to the mosaic."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# import rasterio's mask tool\n",
+    "from rasterio.mask import mask\n",
+    "\n",
+    "# apply mask with crop=True to cut to boundary\n",
+    "with rasterio.open('data/mosaic.tif') as mosaic:\n",
+    "    clipped, transform = mask(mosaic, aoi, crop=True)\n",
+    "    \n",
+    "# See the results!\n",
+    "show(clipped)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Congratulations! You've created a clipped mosaic, showing only the imagery that falls within our area of interest.\n",
+    "\n",
+    "From here, the only thing left to do is save our results to a final output GeoTIFF."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# save the output to a final GeoTIFF\n",
+    "\n",
+    "# use the metadata from our original mosaic\n",
+    "meta = mosaic.meta.copy()\n",
+    "\n",
+    "# update metadata with new, clipped mosaic's boundaries\n",
+    "meta.update({\"transform\": transform,\n",
+    "    \"height\":clipped.shape[1],\n",
+    "    \"width\":clipped.shape[2]})\n",
+    "\n",
+    "# write the output to a GeoTIFF\n",
+    "with rasterio.open('data/clipped_mosaic.tif', 'w', **meta) as dst:\n",
+    "    dst.write(masked)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Final results\n",
+    "\n",
+    "If you like, you might take a closer look at the final clipped mosaic using QGIS:"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "![Final clipped mosaic viewed in QGIS](images/final_in_qgis.png)"
+   ]
+  }
+ ],
+ "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.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}