MovieTime.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# See requirements.txt to set up your dev environment.\n",
    "import os\n",
    "import cv2\n",
    "import sys\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, ogr, osr\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",
    "#from scipy import ndimage\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\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's pull it all together to do something cool.\n",
    "* Let's reuse a lot of our code to make a movie of our travel around Portland.\n",
    "* We'll first select a bunch of recent scenes, activate, and download them.\n",
    "* After that we'll create a mosaic, a path, and trace the path through the moasic. \n",
    "* We'll use the path to crop subregions, save them as images, and create a video.\n",
    "* First step is to trace our AOI and a path through it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "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": [
    "# Query the API\n",
    "* Now we'll save the geometry for our AOI and the path.\n",
    "* We'll also filter and cleanup our data just like before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "areaAOI = AOIs[1][\"geometry\"]\n",
    "pathAOI = AOIs[2][\"geometry\"]\n",
    "\n",
    "aoi_file =\"portland.geojson\" \n",
    "with open(aoi_file,\"w\") as f:\n",
    "    f.write(json.dumps(areaAOI))\n",
    "# build a query using the AOI and\n",
    "# a cloud_cover filter that excludes 'cloud free' scenes\n",
    "\n",
    "old = datetime.datetime(year=2017,month=1,day=1)\n",
    "\n",
    "query = filters.and_filter(\n",
    "    filters.geom_filter(areaAOI),\n",
    "    filters.range_filter('cloud_cover', lt=5),\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": [
    "# Just like before we clean up our data and distill it down to just the scenes we want."
   ]
  },
  {
   "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(areaAOI)\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()\n",
    "\n",
    "# Now let's get it down to just good, recent, clear scenes\n",
    "clear = scenes['cloud_cover']<0.4\n",
    "good = scenes['quality_category']==\"standard\"\n",
    "recent = scenes[\"acquired\"] > datetime.date(year=2017,month=1,day=1)\n",
    "partial_coverage = scenes[\"overlap\"] > 10\n",
    "good_scenes = scenes[(good&clear&recent&partial_coverage)]\n",
    "print good_scenes\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# To make sure we are good we'll visually inspect the scenes in our slippy map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first create a list of colors\n",
    "colors = [\"#ff0000\",\"#00ff00\",\"#0000ff\",\"#ffff00\",\"#ff00ff\",\"#00ffff\",\"#ff0000\",\"#00ff00\",\"#0000ff\",\"#ffff00\",\"#ff00ff\",\"#00ffff\"]\n",
    "# grab our scenes from the geometry/footprint geojson\n",
    "# Chane this number as needed\n",
    "footprints = good_scenes[0:10][\"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':areaAOI,\"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 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This is from the previous notebook. We are just activating and downloading scenes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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",
    "    return True\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\n",
    "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": [
    "# Perform the actual activation ... go get coffee"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_get = good_scenes[\"id\"][0:10].tolist()\n",
    "to_get = sorted(to_get)\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 productfor p in labels:\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": [
    "# Downloand the scenes"
   ]
  },
  {
   "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 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 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now, just like before, we will mosaic those scenes.\n",
    "* It is easier to call out using subprocess and use the command line util.\n",
    "* Just iterate through the files and drop them into a single file portland_mosaic.tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subprocess.call([\"rm\",\"portland_mosaic.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 = \"portland_mosaic.tif\"\n",
    "for tiff in tiff_files:\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 at what we got"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "merged = load_image3(output_mosaic)\n",
    "plt.figure(0,figsize=(18,18))\n",
    "plt.imshow(merged)\n",
    "plt.title(\"merged\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now we are going to write a quick crop function.\n",
    "* this function takes in a, scene, a center position, and the width and height of a window.\n",
    "* We'll use numpy slice notation to make the crop.\n",
    "* Let's pick a spot and see what we get."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def crop_to_area(scene,x_c,y_c,w,h):\n",
    "    tlx = x_c-(w/2)\n",
    "    tly = y_c-(h/2)\n",
    "    brx = x_c+(w/2)\n",
    "    bry = y_c+(h/2)\n",
    "    return scene[tly:bry,tlx:brx,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(0,figsize=(3,4))\n",
    "plt.imshow(crop_to_area(merged,3000,3000,640,480))\n",
    "plt.title(\"merged\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now to figure out how our lat/long values map to pixels.\n",
    "* The next thing we need is a way to map from a lat and long in our slippy map to the pixel position in our image. \n",
    "* We'll use what we know about the lat/long of the corners of our image to do that. \n",
    "* We'll ask GDAL to tell us the extents of our scene and the geotransofrm.\n",
    "* We'll then apply the GeoTransform from GDAL to the coordinates that are the extents of our scene. \n",
    "* Now we have the corners of our scene in Lat/Long"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Liberally borrowed from this example\n",
    "# https://gis.stackexchange.com/questions/57834/how-to-get-raster-corner-coordinates-using-python-gdal-bindings\n",
    "def GetExtent(gt,cols,rows):\n",
    "    \"\"\"\n",
    "    Get the list of corners in our output image in the format \n",
    "    [[x,y],[x,y],[x,y]]\n",
    "    \"\"\"\n",
    "    ext=[]\n",
    "    # for the corners of the image\n",
    "    xarr=[0,cols]\n",
    "    yarr=[0,rows]\n",
    "\n",
    "    for px in xarr:\n",
    "        for py in yarr:\n",
    "            # apply the geo coordiante transform \n",
    "            # using the affine transform we got from GDAL\n",
    "            x=gt[0]+(px*gt[1])+(py*gt[2])\n",
    "            y=gt[3]+(px*gt[4])+(py*gt[5])\n",
    "            ext.append([x,y])\n",
    "        yarr.reverse()\n",
    "    return ext\n",
    "\n",
    "def ReprojectCoords(coords,src_srs,tgt_srs):\n",
    "    trans_coords=[]\n",
    "    # create a transform object from the source and target ref system\n",
    "    transform = osr.CoordinateTransformation( src_srs, tgt_srs)\n",
    "    for x,y in coords:\n",
    "        # transform the points\n",
    "        x,y,z = transform.TransformPoint(x,y)\n",
    "        # add it to the list. \n",
    "        trans_coords.append([x,y])\n",
    "    return trans_coords"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Here we'll call the functions we wrote.\n",
    "* First we open the scene and get the width and height.\n",
    "* Then from the geotransorm we'll reproject those points to lat and long. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TLDR: pixels => UTM coordiantes => Lat Long \n",
    "raster=output_mosaic\n",
    "# Load the GDAL File\n",
    "ds=gdal.Open(raster)\n",
    "# get the geotransform\n",
    "gt=ds.GetGeoTransform()\n",
    "# get the width and height of our image\n",
    "cols = ds.RasterXSize\n",
    "rows = ds.RasterYSize\n",
    "# Generate the coordinates of our image in utm\n",
    "ext=GetExtent(gt,cols,rows)\n",
    "# get the spatial referencec object \n",
    "src_srs=osr.SpatialReference()\n",
    "# get the data that will allow us to move from UTM to Lat Lon. \n",
    "src_srs.ImportFromWkt(ds.GetProjection())\n",
    "tgt_srs = src_srs.CloneGeogCS()\n",
    "extents = ReprojectCoords(ext,src_srs,tgt_srs)\n",
    "print extents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now we'll do a bit of hack.\n",
    "* That bit above is precise but complext, we are going to make everything easier to think about. \n",
    "* We are going to linearize our scene, which isn't perfect, but good enough for our application.\n",
    "* What this function does is take in a given lat,long, the size of the image, and the extents as lat,lon coordinates.\n",
    "* For a given pixel we map it's x and y values to the value between a given lat and long and return the results.\n",
    "* Now we can ask, for a given lat,long pair what is the corresponding pixel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def poor_mans_lat_lon_2_pix(lon,lat,w,h,extents):\n",
    "    # split up our lat and longs \n",
    "    lats = [e[1] for e in extents]\n",
    "    lons = [e[0] for e in extents]\n",
    "    # calculate our scene extents  max and min\n",
    "    lat_max = np.max(lats)\n",
    "    lat_min = np.min(lats)    \n",
    "    lon_max = np.max(lons)\n",
    "    lon_min = np.min(lons) \n",
    "    # calculate the difference between our start point\n",
    "    # and our minimum\n",
    "    lat_diff = lat-lat_min\n",
    "    lon_diff = lon-lon_min\n",
    "    # create the linearization\n",
    "    lat_r = float(h)/(lat_max-lat_min)\n",
    "    lon_r = float(w)/(lon_max-lon_min) \n",
    "    # generate the results. \n",
    "    return int(lat_r*lat_diff),int(lon_r*lon_diff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's check our work\n",
    "* First we'll create a draw point function that just puts a red dot at given pixel.\n",
    "* We'll get our scene, and map all of the lat/long points in our path to pixel values.\n",
    "* Finally we'll load our image, plot the points and show our results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_point(x,y,img,t=40):\n",
    "    h,w,d = img.shape\n",
    "    y = h-y\n",
    "    img[(y-t):(y+t),(x-t):(x+t),:] = [255,0,0]\n",
    "h,w,c = merged.shape\n",
    "waypoints = [poor_mans_lat_lon_2_pix(point[0],point[1],w,h,extents) for point in pathAOI[\"coordinates\"]]\n",
    "print waypoints\n",
    "merged = load_image3(output_mosaic)\n",
    "[draw_point(pt[1],pt[0],merged) for pt in waypoints]\n",
    "plt.figure(0,figsize=(18,18))\n",
    "plt.imshow(merged)\n",
    "plt.title(\"merged\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now things get interesting....\n",
    "* Our path is just a few waypoint but to make a video we need just about every point between our waypoints.\n",
    "* To get all of the points between our waypoints we'll have to write a little interpolation script. \n",
    "* Interpolation is just a fancy word for nicely space points bewteen or waypoints, we'll call the space between each point as our \"velocity.\"\n",
    "* If we were really slick we could define a heading vector and and build a spline so the camera faces the direction of heading. Our approach is fine as the top of the frame is always North, which makes reckoning easy.\n",
    "* Once we have our interpolation function all we need to do is to crop our large mosaic at each point in our interpolation point list and save it in a sequential file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def interpolate_waypoints(waypoints,velocity=10.0):\n",
    "    retVal = []\n",
    "    last_pt = waypoints[0]\n",
    "    # for each point in our waypoints except the first\n",
    "    for next_pt in waypoints[1:]:\n",
    "        # calculate distance between the points\n",
    "        distance = np.sqrt((last_pt[0]-next_pt[0])**2+(last_pt[1]-next_pt[1])**2)\n",
    "        # use our velocity to calculate the number steps.\n",
    "        steps = np.ceil(distance/velocity)\n",
    "        # linearly space points between the two points on our line\n",
    "        xs = np.array(np.linspace(last_pt[0],next_pt[0],steps),dtype='int64')\n",
    "        ys = np.array(np.linspace(last_pt[1],next_pt[1],steps),dtype='int64')\n",
    "        # zip the points together\n",
    "        retVal += zip(xs,ys)\n",
    "        # move to the next point\n",
    "        last_pt = next_pt\n",
    "    return retVal\n",
    "\n",
    "def build_scenes(src,waypoints,window=[640,480],path=\"./movie/\"):\n",
    "    count = 0 \n",
    "    # Use opencv to change the color space of our image.\n",
    "    src = cv2.cvtColor(src, cv2.COLOR_BGR2RGB)\n",
    "    # define half our sampling window. \n",
    "    w2 = window[0]/2\n",
    "    h2 = window[1]/2\n",
    "    # for our source image get the width and height\n",
    "    h,w,d = src.shape\n",
    "    for pt in waypoints:\n",
    "        # for each point crop the area out.\n",
    "        # the y value of our scene is upside down. \n",
    "        temp = crop_to_area(src,pt[1],h-pt[0],window[0],window[1])\n",
    "        # If we happen to hit the border of the scene, just skip\n",
    "        if temp.shape[0]*temp.shape[1]== 0:\n",
    "            # if we have an issue, just keep plugging along\n",
    "            continue\n",
    "        # Resample the image a bit, this just makes things look nice. \n",
    "        temp = cv2.resize(temp, (int(window[0]*0.75), int(window[1]*.75))) \n",
    "        # create a file name\n",
    "        fname = os.path.abspath(path+\"img{num:06d}.png\".format(num=count))\n",
    "        # Save it\n",
    "        cv2.imwrite(fname,temp)\n",
    "        count += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Before we generate our video frames, let's check our work\n",
    "* We'll load our image. \n",
    "* Build the interpolated waypoints list.\n",
    "* Draw the points on the image using our draw_point method.\n",
    "* Plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the image\n",
    "merged = load_image3(output_mosaic)\n",
    "# interpolate the waypoints\n",
    "interp = interpolate_waypoints(waypoints)\n",
    "# draw them on our scene\n",
    "[draw_point(pt[1],pt[0],merged) for pt in interp]\n",
    "# display the scene\n",
    "plt.figure(0,figsize=(18,18))\n",
    "plt.imshow(merged)\n",
    "plt.title(\"merged\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now let's re-load the image and run the scene maker."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "os.system(\"rm ./movie/*.png\")\n",
    "merged = load_image3(output_mosaic)\n",
    "build_scenes(merged,interp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finally, let's make a movie. \n",
    "* Our friend AVConv, which is like ffmpeg is a handy command line util for transcoding video.\n",
    "* AVConv can also convert a series of images into a video and vice versa.\n",
    "* We'll set up our command and use subprocess to make the call. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# avconv -framerate 30 -f image2 -i ./movie/img%06d.png -b 65536k out.mpg\n",
    "os.system(\"rm ./movie/*.png\")\n",
    "framerate = 30\n",
    "output = \"out.mpg\"\n",
    "command = [\"avconv\",\"-framerate\", str(framerate), \"-f\", \"image2\", \"-i\", \"./movie/img%06d.png\", \"-b\", \"65536k\", output]\n",
    "os.system(\" \".join(command))"
   ]
  },
  {
   "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": 2
}