TheBasics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# You'll probably want to set our data rate higher for this notebook. \n",
    "# follow: http://stackoverflow.com/questions/43288550/iopub-data-rate-exceeded-when-viewing-image-in-jupyter-notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup \n",
    "Let's setup our environment. We'll pull in the the usual gis suspects and setup a leaflet map, read our API keys from a json file, and setup our Planet client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# See requirements.txt to set up your dev environment.\n",
    "import sys\n",
    "import os\n",
    "import json\n",
    "import scipy\n",
    "import urllib\n",
    "import datetime \n",
    "import urllib3\n",
    "import rasterio\n",
    "import subprocess\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from osgeo import gdal\n",
    "from planet import api\n",
    "from planet.api import filters\n",
    "from traitlets import link\n",
    "import rasterio.tools.mask as rio_mask\n",
    "from shapely.geometry import mapping, shape\n",
    "from IPython.display import display, Image, HTML\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "urllib3.disable_warnings()\n",
    "from ipyleaflet import (\n",
    "    Map,\n",
    "    Marker,\n",
    "    TileLayer, ImageOverlay,\n",
    "    Polyline, Polygon, Rectangle, Circle, CircleMarker,\n",
    "    GeoJSON,\n",
    "    DrawControl\n",
    ")\n",
    "\n",
    "%matplotlib inline\n",
    "# will pick up api_key via environment variable PL_API_KEY\n",
    "# but can be specified using `api_key` named argument\n",
    "api_keys = json.load(open(\"apikeys.json\",'r'))\n",
    "client = api.ClientV1(api_key=api_keys[\"PLANET_API_KEY\"])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make a slippy map to get GeoJSON\n",
    "\n",
    "* The planet API allows you to query using a [geojson](https://en.wikipedia.org/wiki/GeoJSON) which is a special flavor of json.\n",
    "* We are going to create a slippy map using leaflet and apply the Planet 2017 Q1 mosaic as the basemap. This requires our api key.\n",
    "* We are going to add a special draw handler that shoves a draw region into a object so we get the geojson.\n",
    "* If you don't want to do this, or need a fixed query try [geojson.io](http://geojson.io/#map=2/20.0/0.0)\n",
    "* To install and run:\n",
    "```\n",
    "$ pip install ipyleaflet\n",
    "$ jupyter nbextension enable --py --sys-prefix ipyleaflet\n",
    "$ jupyter nbextension enable --py --sys-prefix widgetsnbextension\n",
    "```\n",
    "* [More information](https://github.com/ellisonbg/ipyleaflet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Basemap Mosaic (v1 API)\n",
    "mosaicsSeries = 'global_quarterly_2017q1_mosaic'\n",
    "# Planet tile server base URL (Planet Explorer Mosaics Tiles)\n",
    "mosaicsTilesURL_base = 'https://tiles0.planet.com/experimental/mosaics/planet-tiles/' + mosaicsSeries + '/gmap/{z}/{x}/{y}.png'\n",
    "# Planet tile server url\n",
    "mosaicsTilesURL = mosaicsTilesURL_base + '?api_key=' + api_keys[\"PLANET_API_KEY\"]\n",
    "# Map Settings \n",
    "# Define colors\n",
    "colors = {'blue': \"#009da5\"}\n",
    "# Define initial map center lat/long\n",
    "center = [45.5231, -122.6765]\n",
    "# Define initial map zoom level\n",
    "zoom = 11\n",
    "# Set Map Tiles URL\n",
    "planetMapTiles = TileLayer(url= mosaicsTilesURL)\n",
    "# Create the map\n",
    "m = Map(\n",
    "    center=center, \n",
    "    zoom=zoom,\n",
    "    default_tiles = planetMapTiles # Uncomment to use Planet.com basemap\n",
    ")\n",
    "# Define the draw tool type options\n",
    "polygon = {'shapeOptions': {'color': colors['blue']}}\n",
    "rectangle = {'shapeOptions': {'color': colors['blue']}} \n",
    "\n",
    "# Create the draw controls\n",
    "# @see https://github.com/ellisonbg/ipyleaflet/blob/master/ipyleaflet/leaflet.py#L293\n",
    "dc = DrawControl(\n",
    "    polygon = polygon,\n",
    "    rectangle = rectangle\n",
    ")\n",
    "# Initialize an action counter variable\n",
    "actionCount = 0\n",
    "AOIs = {}\n",
    "\n",
    "# Register the draw controls handler\n",
    "def handle_draw(self, action, geo_json):\n",
    "    # Increment the action counter\n",
    "    global actionCount\n",
    "    actionCount += 1\n",
    "    # Remove the `style` property from the GeoJSON\n",
    "    geo_json['properties'] = {}\n",
    "    # Convert geo_json output to a string and prettify (indent & replace ' with \")\n",
    "    geojsonStr = json.dumps(geo_json, indent=2).replace(\"'\", '\"')\n",
    "    AOIs[actionCount] = json.loads(geojsonStr)\n",
    "    \n",
    "# Attach the draw handler to the draw controls `on_draw` event\n",
    "dc.on_draw(handle_draw)\n",
    "m.add_control(dc)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Querying the Planet API.\n",
    "* First we'll grab our geojson area of interest (AOI) and use it to construct a query.\n",
    "* We'll then build a search to search that area looking for PSScene3Band\n",
    "* We have lots of products: RapidEye, PlanetScope (PS) 3 and 4 band, LandSat, and Sentinel are all possible.\n",
    "* Once we have our query, we'll do the search. We will then iterate over the results, slurp up the data, and put them in a pandas data frame for easy sorting.\n",
    "* We'll print the first few so we're sure it works. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print AOIs[1]\n",
    "myAOI = AOIs[1][\"geometry\"]\n",
    "\n",
    "# build a query using the AOI and\n",
    "# a cloud_cover filter that excludes 'cloud free' scenes\n",
    "\n",
    "old = datetime.datetime(year=2013,month=1,day=1)\n",
    "\n",
    "query = filters.and_filter(\n",
    "    filters.geom_filter(myAOI),\n",
    "    filters.range_filter('cloud_cover', lt=50),\n",
    "    filters.date_range('acquired', gt=old)\n",
    ")\n",
    "\n",
    "# build a request for only PlanetScope imagery\n",
    "request = filters.build_search_request(\n",
    "    query, item_types=['PSScene3Band']\n",
    ")\n",
    "\n",
    "# if you don't have an API key configured, this will raise an exception\n",
    "result = client.quick_search(request)\n",
    "scenes = []\n",
    "planet_map = {}\n",
    "for item in result.items_iter(limit=500):\n",
    "    planet_map[item['id']]=item\n",
    "    props = item['properties']\n",
    "    props[\"id\"] = item['id']\n",
    "    props[\"geometry\"] = item[\"geometry\"]\n",
    "    props[\"thumbnail\"] = item[\"_links\"][\"thumbnail\"]\n",
    "    scenes.append(props)\n",
    "scenes = pd.DataFrame(data=scenes)\n",
    "display(scenes)\n",
    "print len(scenes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cleanup\n",
    "* The data we got back is good, but we need some more information\n",
    "* We got back big scenes, but we only care about our area of interest. The scene may not cover the whole area of interest.\n",
    "* We can use the [Shapely](http://toblerity.org/shapely/manual.html) library to quickly figure out how much each scene overlaps our AOI\n",
    "* We will convert our AOI and the geometry of each scene to calculate overlap using a shapely call.\n",
    "* The returned acquisition, publish, and update times are strings, we'll convert them to datatime objects so we wan search."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now let's clean up the datetime stuff\n",
    "# make a shapely shape from our aoi\n",
    "portland = shape(myAOI)\n",
    "footprints = []\n",
    "overlaps = []\n",
    "# go through the geometry from our api call, convert to a shape and calculate overlap area.\n",
    "# also save the shape for safe keeping\n",
    "for footprint in scenes[\"geometry\"].tolist():\n",
    "    s = shape(footprint)\n",
    "    footprints.append(s)\n",
    "    overlap = 100.0*(portland.intersection(s).area / portland.area)\n",
    "    overlaps.append(overlap)\n",
    "# take our lists and add them back to our dataframe\n",
    "scenes['overlap'] = pd.Series(overlaps, index=scenes.index)\n",
    "scenes['footprint'] = pd.Series(footprints, index=scenes.index)\n",
    "# now make sure pandas knows about our date/time columns.\n",
    "scenes[\"acquired\"] = pd.to_datetime(scenes[\"acquired\"])\n",
    "scenes[\"published\"] = pd.to_datetime(scenes[\"published\"])\n",
    "scenes[\"updated\"] = pd.to_datetime(scenes[\"updated\"])\n",
    "scenes.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filtering our search using pandas.\n",
    "* Using our dataframe we will filter the scenes to just what we want.\n",
    "* First we want scenes with less than 10% clouds.\n",
    "* Second we want standard quality images. Test images may not be high quality.\n",
    "* Third well only look for scenes since January.\n",
    "* Finally we will create a new data frame with our queries and print the results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's get it down to just good, recent, clear scenes\n",
    "clear = scenes['cloud_cover']<0.1\n",
    "good = scenes['quality_category']==\"standard\"\n",
    "recent = scenes[\"acquired\"] > datetime.date(year=2017,month=1,day=1)\n",
    "partial_coverage = scenes[\"overlap\"] > 30\n",
    "good_scenes = scenes[(good&clear&recent&partial_coverage)]\n",
    "display(good_scenes)\n",
    "print len(good_scenes)\n",
    "\n",
    "# Now let's get it down to just good, recent, clear scenes\n",
    "clear = scenes['cloud_cover']<0.5\n",
    "good = scenes['quality_category']==\"standard\"\n",
    "all_time = scenes[\"acquired\"] > datetime.date(year=2014,month=1,day=1)\n",
    "full_coverage = scenes[\"overlap\"] >= 60\n",
    "all_scenes = scenes[(good&clear&all_time&full_coverage)]\n",
    "display(all_scenes)\n",
    "print len(all_scenes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing scene foot prints overlap with our AOI\n",
    "* We know these scenes intersect with our AOI, but we aren't quite sure about the geometry.\n",
    "* We are going to plot our scene footprints and original AOI on our slippy map.\n",
    "* To do this we create GeoJson objects with properties. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first create a list of colors\n",
    "colors = [\"#ff0000\",\"#00ff00\",\"#0000ff\",\"#ffff00\",\"#ff00ff\",\"#00ffff\"]\n",
    "# grab our scenes from the geometry/footprint geojson\n",
    "footprints = good_scenes[\"geometry\"].tolist()\n",
    "# for each footprint/color combo\n",
    "for footprint,color in zip(footprints,colors):\n",
    "    # create the leaflet object\n",
    "    feat = {'geometry':footprint,\"properties\":{\n",
    "            'style':{'color': color,'fillColor': color,'fillOpacity': 0.2,'weight': 1}},\n",
    "            'type':u\"Feature\"}\n",
    "    # convert to geojson\n",
    "    gjson = GeoJSON(data=feat)\n",
    "    # add it our map\n",
    "    m.add_layer(gjson)\n",
    "# now we will draw our original AOI on top \n",
    "feat = {'geometry':myAOI,\"properties\":{\n",
    "            'style':{'color': \"#FFFFFF\",'fillColor': \"#FFFFFF\",'fillOpacity': 0.5,'weight': 1}},\n",
    "            'type':u\"Feature\"}\n",
    "gjson = GeoJSON(data=feat)\n",
    "m.add_layer(gjson)   \n",
    "m \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's see what we got. \n",
    "* The API returns a handy thumbnail link.\n",
    "* Let's tell jupyter to show it.\n",
    "* You may need to login to planet explorer to have auth. \n",
    "    * If this is the case just print the urls and paste them into your browser."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imgs = []\n",
    "# loop through our thumbnails and add display them\n",
    "for img in good_scenes[\"thumbnail\"].tolist():\n",
    "    imgs.append(Image(url=img))\n",
    "    print img\n",
    "display(*imgs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Product Activation and Downloading\n",
    "* There are two things we need to know, the satellite type (asset) and image type (product).\n",
    "* Full resolution uncompressed satellite images are *big* and there are lots of ways to view them.\n",
    "* For this reason Planet generally keeps images in their native format and only processes them on customer requests. There is some caching of processed scenes, but this is the exception not the rule.\n",
    "* All images must be activated prior to downloading and this can take some time based on demand.\n",
    "* Additionally we need to determine what sort of product we want to download. Generally speaking there are three kinds of scenes:\n",
    "    * Analytic - multi-band full resolution images that have not been processed. These are like raw files for DSLR camers.\n",
    "    * Visual - these are color corrected rectified tifs. If you are just starting out this is your best call.\n",
    "    * UDM - Usable data mask. This mask can be used to find bad pixels and columns and to mask out areas with clouds.\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def get_products(client, scene_id, asset_type='PSScene3Band'):    \n",
    "    \"\"\"\n",
    "    Ask the client to return the available products for a \n",
    "    given scene and asset type. Returns a list of product \n",
    "    strings\n",
    "    \"\"\"\n",
    "    out = client.get_assets_by_id(asset_type,scene_id)\n",
    "    temp = out.get()\n",
    "    return temp.keys()\n",
    "\n",
    "def activate_product(client, scene_id, asset_type=\"PSScene3Band\",product=\"analytic\"):\n",
    "    \"\"\"\n",
    "    Activate a product given a scene, an asset type, and a product.\n",
    "    \n",
    "    On success return the return value of the API call and an activation object\n",
    "    \"\"\"\n",
    "    temp = client.get_assets_by_id(asset_type,scene_id)  \n",
    "    products = temp.get()\n",
    "    if( product in products.keys() ):\n",
    "        return client.activate(products[product]),products[product]\n",
    "    else:\n",
    "        return None \n",
    "\n",
    "def download_and_save(client,product):\n",
    "    \"\"\"\n",
    "    Given a client and a product activation object download the asset. \n",
    "    This will save the tiff file in the local directory and return its \n",
    "    file name. \n",
    "    \"\"\"\n",
    "    out = client.download(product)\n",
    "    fp = out.get_body()\n",
    "    fp.write()\n",
    "    return fp.name\n",
    "\n",
    "def scenes_are_active(scene_list):\n",
    "    \"\"\"\n",
    "    Check if all of the resources in a given list of\n",
    "    scene activation objects is read for downloading.\n",
    "    \"\"\"\n",
    "    retVal = True\n",
    "    for scene in scene_list:\n",
    "        if scene[\"status\"] != \"active\":\n",
    "            print \"{} is not ready.\".format(scene)\n",
    "            return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scenes ACTIVATE!\n",
    "* Given our good scenes list we will convert the data frame \"id\" column into a list and activate every item in that list. \n",
    "* For this example we are going to default to using a 3Band visual product but I have included some four band methods to help you out.\n",
    "* Activation usually takes about 5-15 minutes so get some coffee."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_get = good_scenes[\"id\"].tolist()\n",
    "activated = []\n",
    "# for each scene to get\n",
    "for scene in to_get:\n",
    "    # get the product \n",
    "    product_types = get_products(client,scene)\n",
    "    for p in product_types:\n",
    "        # if there is a visual product\n",
    "        if p == \"visual\": # p == \"basic_analytic_dn\"\n",
    "            print \"Activating {0} for scene {1}\".format(p,scene)\n",
    "            # activate the product\n",
    "            _,product = activate_product(client,scene,product=p)\n",
    "            activated.append(product)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Download Scenes\n",
    "* In this section we will see if our scenes have been activated.\n",
    "* If they are activated the client object will have its status flag set to active.\n",
    "* Once that is done we will then save the scenes to the local directory.\n",
    "* A smart engineer would set a path variable to store these files and check if the asset has already been downloaded prior to downloading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tiff_files = []\n",
    "asset_type = \"_3B_Visual\"\n",
    "# check if our scenes have been activated\n",
    "if True: #scenes_are_active(activated):\n",
    "    for to_download,name in zip(activated,to_get):\n",
    "        # create the product name\n",
    "        name = name + asset_type + \".tif\"\n",
    "        # if the product exists locally\n",
    "        if( os.path.isfile(name) ):\n",
    "            # do nothing \n",
    "            print \"We have scene {0} already, skipping...\".format(name)\n",
    "            tiff_files.append(name)\n",
    "        elif to_download[\"status\"] == \"active\":\n",
    "            # otherwise download the product\n",
    "            print \"Downloading {0}....\".format(name)\n",
    "            fname = download_and_save(client,to_download)\n",
    "            tiff_files.append(fname)\n",
    "            print \"Download done.\"\n",
    "        else:\n",
    "            print \"Could not download, still activating\"\n",
    "else:\n",
    "    print \"Scenes aren't ready yet\"\n",
    "\n",
    "print tiff_files \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loading Images\n",
    "* There are a varitety of ways to load tif data including Rasterio, GDAL, OpenCV, SKImage. \n",
    "* Today we are going to use rasterio and load each channel into a numpy array.\n",
    "* Since the visual 3Band products are rotated we can also open a mask layer for processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def load_image4(filename):\n",
    "    \"\"\"Return a 4D (r, g, b, nir) numpy array with the data in the specified TIFF filename.\"\"\"\n",
    "    path = os.path.abspath(os.path.join('./', filename))\n",
    "    if os.path.exists(path):\n",
    "        with rasterio.open(path) as src:\n",
    "            b, g, r, nir = src.read()\n",
    "            return np.dstack([r, g, b, nir])\n",
    "        \n",
    "def load_image3(filename):\n",
    "    \"\"\"Return a 3D (r, g, b) numpy array with the data in the specified TIFF filename.\"\"\"\n",
    "    path = os.path.abspath(os.path.join('./', filename))\n",
    "    if os.path.exists(path):\n",
    "        with rasterio.open(path) as src:\n",
    "            b,g,r,mask = src.read()\n",
    "            return np.dstack([b, g, r])\n",
    "        \n",
    "def get_mask(filename):\n",
    "    \"\"\"Return a 1D mask numpy array with the data in the specified TIFF filename.\"\"\"\n",
    "    path = os.path.abspath(os.path.join('./', filename))\n",
    "    if os.path.exists(path):\n",
    "        with rasterio.open(path) as src:\n",
    "            b,g,r,mask = src.read()\n",
    "            return np.dstack([mask])\n",
    "\n",
    "def rgbir_to_rgb(img_4band):\n",
    "    \"\"\"Convert an RGBIR image to RGB\"\"\"\n",
    "    return img_4band[:,:,:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read Images and Use Matplotlib to show them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_files = []\n",
    "masks = []\n",
    "# load the images and masks\n",
    "for fname in tiff_files[0:2]:\n",
    "    img_files.append(load_image3(fname))\n",
    "    masks.append(get_mask(fname))\n",
    "i = 0\n",
    "# use matplotlib to display the map\n",
    "for img,name in zip(img_files,tiff_files):\n",
    "    plt.figure(i,figsize=(18,36))\n",
    "    plt.imshow(img)\n",
    "    plt.title(name)\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quick Histogram\n",
    "* Next up we'll plot the histogram of the image.\n",
    "* A histogram is just a plot of the number of pixels with a specific intensity for a given color. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy.ma as ma\n",
    "def plot_hist4(img_4band,title=\"\"):\n",
    "    # Plot a four band histogram\n",
    "    r, g, b, nir = img_4band[:, :, 0], img_4band[:, :, 1], img_4band[:, :, 2], img_4band[:, :, 3]\n",
    "    for slice_, name, color in ((r,'r', 'red'),(g,'g', 'green'),(b,'b', 'blue'), (nir, 'nir', 'magenta')):\n",
    "        plt.hist(slice_.ravel(), bins=100, \n",
    "                 range=[0,img_4band.max()], \n",
    "                 label=name, color=color, histtype='step')\n",
    "    plt.title(title)\n",
    "    plt.legend()\n",
    "    \n",
    "def plot_hist3(img_3band,mask,title=\"\"):\n",
    "    # plot a three band histogramwaiter = []\n",
    "    r, g, b = img_3band[:, :, 0], img_3band[:, :, 1], img_3band[:, :, 2]\n",
    "    r = ma.masked_array(r,mask=mask)\n",
    "    g = ma.masked_array(g,mask=mask)\n",
    "    b = ma.masked_array(b,mask=mask)\n",
    "    for slice_, name, color in ((r,'r', 'red'),(g,'g', 'green'),(b,'b', 'blue')):\n",
    "        plt.hist(slice_.ravel(), bins=25, \n",
    "                 range=[0,img_3band.max()], \n",
    "                 label=name, color=color, histtype='step')\n",
    "    plt.title(title)\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0\n",
    "for img,name,mask in zip(img_files,tiff_files,masks):\n",
    "    plt.figure(i,figsize=(9,18))\n",
    "    plot_hist3(img,mask=mask,title=name)\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Decomposing Channels\n",
    "* We can also decompose the channels of the image. \n",
    "* Sometimes it is useful to work just in a single channel.\n",
    "* Other times channels can be used to do useful things, like filter out clouds.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_bands4(img,title=\"\",i=0):\n",
    "    fig = plt.figure(i)\n",
    "    fig.set_size_inches(24, 3)\n",
    "    r, g, b, nir = img[:, :, 0], img[:, :, 1], img[:, :, 2], img[:, :, 3]\n",
    "    fig.suptitle(title)\n",
    "    for i, (x, c) in enumerate(((r, 'r'), (g, 'g'), (b, 'b'), (nir, 'near-ir'))):\n",
    "        a = fig.add_subplot(1, 4, i+1)\n",
    "        a.set_title(c)\n",
    "        plt.imshow(x)\n",
    "        \n",
    "def plot_bands3(img,title=\"\",i=0):\n",
    "    fig = plt.figure(i)\n",
    "    fig.set_size_inches(24, 5)\n",
    "    r, g, b = img[:, :, 0], img[:, :, 1], img[:, :, 2]\n",
    "    fig.suptitle(title)\n",
    "    for i, (x, c) in enumerate(((r, 'r'), (g, 'g'), (b, 'b'))):\n",
    "        a = fig.add_subplot(1, 4, i+1)\n",
    "        a.set_title(c)\n",
    "        plt.imshow(x)        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_bands3(img_files[0],title=tiff_files[0],i=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# But all of these scenes are big, and we want downtown Portland \n",
    "* We can clip all of the scenes to the AOI we selected at the start of the notebook\n",
    "* First we'll dump the geojson to a file.\n",
    "* Since geospatial data is \"big\" we often work with files and get stuff out of memory ASAP.\n",
    "* For each of our scenes we'll create a 'clip' file.\n",
    "* We will use a tool called GDAL to clip the scene to our AOI\n",
    "* GDAL stands for [Geospatial Data Abstraction Library](http://www.gdal.org/)\n",
    "* GDAL is a C++ library that is often run from the command line, but it does have SWIG bindings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aoi_file =\"portland.geojson\" \n",
    "# write our input AOI to a geojson file.\n",
    "with open(aoi_file,\"w\") as f:\n",
    "    f.write(json.dumps(myAOI))\n",
    "    \n",
    "\n",
    "# create our full input and output names\n",
    "clip_names = [os.path.abspath(tiff[:-4]+\"_clip\"+\".tif\") for tiff in tiff_files]\n",
    "full_tif_files = [os.path.abspath(\"./\"+tiff) for tiff in tiff_files]\n",
    "\n",
    "for in_file,out_file in zip(tiff_files,clip_names):\n",
    "    commands = [\"gdalwarp\", # t\n",
    "           \"-t_srs\",\"EPSG:3857\",\n",
    "           \"-cutline\",aoi_file,\n",
    "           \"-crop_to_cutline\",\n",
    "           \"-tap\",\n",
    "            \"-tr\", \"3\", \"3\"\n",
    "           \"-overwrite\"]\n",
    "    subprocess.call([\"rm\",out_file])\n",
    "    commands.append(in_file)\n",
    "    commands.append(out_file)\n",
    "    print \" \".join(commands)\n",
    "    subprocess.call(commands)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Awesome, Let's take a look at what we got."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clip_img_files = [load_image3(fname) for fname in clip_names]\n",
    "    \n",
    "i = 0\n",
    "for img,name in zip(clip_img_files,clip_names):\n",
    "    plt.figure(i,figsize=(6,12))\n",
    "    plt.imshow(img)\n",
    "    plt.title(name)\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hrm... that's not right.\n",
    "* You'll notice that a lot of these scenes don't fill our AOI.\n",
    "* A lot of theses images were taken roughly at the same time.\n",
    "* We should try to merge these scenes together to make one big scene.\n",
    "* This process is called mosaicking, and GDAL can help. \n",
    "* We will call GDAL from the command line using subprocess to do this for us.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subprocess.call([\"rm\",\"merged.tif\"])\n",
    "commands = [\"gdalwarp\", # t\n",
    "           \"-t_srs\",\"EPSG:3857\",\n",
    "           \"-cutline\",aoi_file,\n",
    "           \"-crop_to_cutline\",\n",
    "           \"-tap\",\n",
    "            \"-tr\", \"3\", \"3\"\n",
    "           \"-overwrite\"]\n",
    "output_mosaic = \"merged.tif\"\n",
    "for tiff in tiff_files[0:2]:\n",
    "    commands.append(tiff)\n",
    "commands.append(output_mosaic)\n",
    "print \" \".join(commands)\n",
    "subprocess.call(commands)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's take a look.... looks much better"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged = load_image3(\"./merged.tif\")\n",
    "plt.figure(i,figsize=(6,12))\n",
    "plt.imshow(merged)\n",
    "plt.title(\"merged\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now let's pull it all together to do something interesting.\n",
    "* First we'll download and activate all of our targe scenes.\n",
    "* Then we'll clip them using GDAL to the small AOI we selected above.\n",
    "* Finally we'll export them and use that data to make a mosaic. \n",
    "* We'll use [ImageMagick](https://www.imagemagick.org/script/index.php) to convert our tifs to gifs, and our multiple gifs to an animated gif. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Activate\n",
    "to_get = all_scenes[\"id\"].tolist()\n",
    "activated = []\n",
    "for scene in to_get:\n",
    "    product_types = get_products(client,scene)\n",
    "    for p in product_types:\n",
    "        if p == \"visual\": # p == \"basic_analytic_dn\"\n",
    "            print \"Activating {0} for scene {1}\".format(p,scene)\n",
    "            _,product = activate_product(client,scene,product=p)\n",
    "            activated.append(product)\n",
    "            \n",
    "# Download  \n",
    "tiff_files = []\n",
    "asset_type = \"_3B_Visual\"\n",
    "if True: #scenes_are_active(activated):\n",
    "    for to_download,name in zip(activated,to_get):\n",
    "        name = name + asset_type + \".tif\"\n",
    "        if( os.path.isfile(name) ):\n",
    "            print \"We have scene {0} already, skipping...\".format(name)\n",
    "            tiff_files.append(name)\n",
    "        elif to_download[\"status\"] == \"active\":\n",
    "            print \"Downloading {0}....\".format(name)\n",
    "            fname = download_and_save(client,to_download)\n",
    "            tiff_files.append(fname)\n",
    "            print \"Download done.\"\n",
    "        else:\n",
    "            print \"Could not download, still activating\"\n",
    "else:\n",
    "    print \"Scenes aren't ready yet\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finally let's process the scenes we just downloaded and make a gif."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tiff_files = sorted(tiff_files)\n",
    "# Create a list of tif file names. \n",
    "for tiff in tiff_files:\n",
    "    clip_names.append(os.path.abspath(tiff[:-4]+\"_clip\"+\".tif\"))\n",
    "\n",
    "full_tif_files = []\n",
    "for tiff in tiff_files:\n",
    "    full_tif_files.append(os.path.abspath(\"./\"+tiff))\n",
    "\n",
    "    # Run GDAL to crop our file down.\n",
    "for in_file,out_file in zip(tiff_files,clip_names):\n",
    "    commands = [\"gdalwarp\", # t\n",
    "           \"-t_srs\",\"EPSG:3857\",\n",
    "           \"-cutline\",aoi_file,\n",
    "           \"-crop_to_cutline\",\n",
    "           \"-tap\",\n",
    "            \"-tr\", \"3\", \"3\"\n",
    "           \"-overwrite\"]\n",
    "    subprocess.call([\"rm\",out_file])\n",
    "    commands.append(in_file)\n",
    "    commands.append(out_file)\n",
    "    print \" \".join(commands)\n",
    "    subprocess.call(commands)    \n",
    "\n",
    "temp_names = []\n",
    "i = 0 \n",
    "# use image magic convert to \n",
    "for in_file in clip_names:\n",
    "    temp_name = \"img{0}.gif\".format(i) \n",
    "    command = [\"convert\", in_file, \"-sample\", \"30x30%\",temp_name]\n",
    "    temp_names.append(temp_name)\n",
    "    i += 1 \n",
    "    subprocess.call(command)\n",
    "\n",
    "magic = \"portland.gif\"\n",
    "last_call = [\"convert\",\"-delay\", \"40\",\"-loop\",\"0\", \"img*.gif\",magic]\n",
    "subprocess.call(last_call)\n",
    "print \"done!\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<img src=\"./XXX.gif\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}