AEOCCG-AccessingAusCoverDataOnTheNCI.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Accessing AusCover Data on the NCI and RDS <br>\n",
    "\n",
    "<center><table>\n",
    "<tr><td><img src='http://portal.tern.org.au/assets/core/img/logos/logo-auscover.png' width=\"200\"></td>\n",
    "<td><img src='https://www.wpcentral.com.au/wp-content/uploads/2013/08/nci-logo1.png' width=\"200\"></td>\n",
    "<td><img src='http://nci.org.au/wp-content/uploads/2015/05/Official-RDS-Logo-RGB.png' width=\"200\"></td></tr>\n",
    "</table></center>\n",
    "\n",
    "#### A Technical Capacity Building Webinar for AEOCCG\n",
    "\n",
    "_Presented by [Peter Scarth](mailto:p.scarth@uq.edu.au?subject=AEOCCG%20webinar%20information) (Joint Remote Sensing Research Program) with lots of help and input from Kelsey  Druken and Claire Trenham from the [NCI](http://nci.org.au/about-nci/contact/nci-staff-2/)._\n",
    "\n",
    "### Abstract\n",
    "\n",
    "Have you ever wondered how to access some of the AusCover remotely sensed data which is available on the NCI using some of the available web services directly within your analysis environment?\n",
    "\n",
    "This webinar will give several examples using direct access using opendap on thredds and gdal/rasterio on thredds and http where you'll be able to query and directly interact with Landsat, MODIS and Himawari datasets in both  iPython notebooks and in GIS packages. By the end of the webinar, you’ll be able to discover some of the many data sets openly available  on the NCI and query them using web services.\n",
    "\n",
    "#### Background on RDS <> AusCover <> NCI <> AGDC?\n",
    " - [AusCover](http://auscover.org.au) - provides a national expert network and a data delivery service for provision of Australian biophysical remote sensing data time-series, continental-scale map products, and selected high-resolution datasets over TERN sites. **It is a virtual network for earth observation data sets across Australia that are typically produced by universities, state government agencies _(e.g. NSW OEH and Qld DSITI)_, and federal  agencies _(e.g. GA, BOM and CSIRO)_ ).**\n",
    " - [Research Data Storage Infrastructure (RDSI)](https://www.rds.edu.au) - supports collaborative access to research data assets of national significance (including national reference collections) through distributed data centre development. **They provide the online storage capacity to host Australian EO data sets online.**\n",
    " - [National Computational Infrastructure (NCI)](http://nci.org.au/) -  is Australia’s national research computing facility including the Southern Hemisphere’s most highly-integrated supercomputer and filesystems, Australia’s highest performance research cloud, and one of the nation’s largest data catalogues. **They provide the compute capacity for some of the Auscover data sets, interactive environments, and cataloging and training services.** \n",
    " - [Australian Geoscience Data Cube (AGDC)](http://www.datacube.org.au/) - provides an integrated gridded data analysis environment for earth observation satellite and related data from multiple satellite and other acquisition systems. **This project facilitates the operational processing and interactive  analysis of national EO data sets and is used to produce the Landsat and MODIS collections on the NCI.**\n",
    "\n",
    "These groups work together to improve the processing, storage, cataloging, discovery and analysis of EO based data across Australia.\n",
    "\n",
    "#### The NCI VDI\n",
    "Data downloading and analysis by many users has potential risks _(apart from the data being too big for this to be feasible!)_. **Bringing scientists to the data can help mitigate these issues by ensuring everyone is working on the same data.** To support this workflow, the NCI runs a Virtual Desktop Environment with:\n",
    " - Tools to support climate data analysis & visualisation\n",
    " - Virtual laboratory to access, process & analyse data\n",
    " - Analyses input data in a consistent format \n",
    " - Workflow tools allow science community to                                        implement own analyses without dealing directly with filesystems & HPC\n",
    " - A range of standard software tools connected to the global Lustre filesystem and HPC\n",
    " -  Desktops: 32GB RAM, 140GB local scratch, 8vCPUs\n",
    " - [How can I get VDI access?](http://nci.org.au/access/getting-access-to-the-national-facility/allocation-scheme/ )\n",
    "\n",
    "### NCI and RDS AusCover Services\n",
    "\n",
    "NCI have a multi-element system for metadata catalogues and data services\n",
    " - [GeoNetwork](http://geonetwork.nci.org.au/geonetwork): Find metadata records (akin to CSIRO DAP)\n",
    " - [THREDDS](http://dap.nci.org.au) Data Service: download or remotely access or view data\n",
    "   - OPeNDAP is one of the protocols served, permits subsetting and remote access to files\n",
    "   - Other protocols include HTTP download and Open Geospatial Consortium Web Services to stream JPEG, TIFF etc.\n",
    " - Geoserver, ERDDAP, Hyrax, others… and filesystem\n",
    " - PROMS (provenance), DOI minting (citation)\n",
    "\n",
    "[AusCover](http://qld.auscover.org.au/public/html/index.html) also uses [THREDDS](http://qld.auscover.org.au/thredds/catalog.html),  [Geoserver](http://qld.auscover.org.au/geoserver/web/), [FTP](ftp://qld.auscover.org.au/) and [HTTP](http://qld.auscover.org.au/public/data/) services to deliver data in a variety of formats, and is working to refine online discovery,  subsetting and reformatting tools for both raster and vector data. Keep an eye out for the new website launching soon.\n",
    "\n",
    "\n",
    "## Presentation Outline\n",
    "\n",
    "For these examples I'll be using a [Jupyter Notebook](http://jupyter.org/) with code in Python.\n",
    " - _The Jupyter Notebook is a web application that allows you to create and share documents that contain live code, equations, visualizations and explanatory text and has support for over 40 programming languages, including Python, R, Julia and Scala_. \n",
    " - If you've never used Jupyter Notebooks before, I highly recommend installing [Anaconda](https://www.continuum.io/downloads)\n",
    " - _As an aside, many of the packages used by JRSRP and partners, such as [RIOS](http://rioshome.org/), [RSGISLIB](http://www.rsgislib.org/) and [PyLidar](http://pylidar.org/) can be installed into this environment from the [OSGEO Conda index](https://conda.anaconda.org/osgeo)_\n",
    "\n",
    "This notebook will outline some simple online interaction with some of the **JRSRP Landsat seasonal mosaics**. We'll treat the data hosted on http://qld.auscover.org.au as files and use the [RasterIO](https://www.mapbox.com/blog/rasterio-announce/) package to interact with the data and undertake some typical remote sensing tasks. Finally we'll build a simple example to extract and analyse a time series of imagery across an agricultural research property in the Burdekin (from ~5 TB of raster data hosted online).\n",
    "\n",
    "Then we'll look at how you'd **access some of the GA Landsat data** produced out of the AGDC and hosted on the NCI using the OPeNDAP protocol via THREDDS. [Link to Hosted Notebook](https://github.com/nci/Notebooks/blob/master/Python_Examples/Python_GDAL_NetCDF.ipynb)\n",
    "\n",
    "Finally we'll check out a pretty cool notebook that uses the NCI THREDDS Data Server and queries the **CSIRO Auscover MODIS** data sets to extract a time series of imagery. [Link to hosted Notebook](https://github.com/nci/Notebooks/blob/master/Data_Access/Using_Siphon/Python_Siphon_II.ipynb)\n",
    "\n",
    "\n",
    "### Additional NCI Resources\n",
    "\n",
    "We won't have time to look at how you query and explore the THREDDS catalog to [access data](https://github.com/nci/Notebooks/blob/master/Data_Access/Using_Thredds/THREDDS_DataAccess.ipynb) or find [WMS and WCS service endpoints](https://github.com/nci/Notebooks/blob/master/Data_Access/Using_Thredds/THREDDS_WMS_WCS.ipynb) so I'd strongly encourage you to follow these links if you want to find out more.\n",
    "\n",
    "Similarly, **accessing many of these data sets is easy using your desktop GIS package**. This will be the topic of another AusCover session later in the year, but for now have a look at the [NCI QGIS examples](https://github.com/nci/Notebooks/tree/master/QGIS_Examples).\n",
    "\n",
    "\n",
    "This all looks a little more tricky than firing up your desktop remote sensing package, but you do get a highly flexible open source analysis environment that give you the ability to perform reproducible research, and operationalise your algorithms nationally with ease.\n",
    "See the links below for training information, more Jupyter notebooks or NCI help:\n",
    " - https://training.nci.org.au \n",
    " - https://github.com/nci/nci-notebooks\n",
    " - http://nci.org.au/user-support/getting-help/\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### Import python packages\n",
    "There are a number of Python packages that help interacting with raster data. Here I use [Rasterio](https://github.com/mapbox/rasterio), a GDAL and Numpy-based Python library designed to make your work with geospatial raster data more productive, more fun — more [Zen](https://www.python.org/dev/peps/pep-0020/).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import requests\n",
    "import rasterio\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Opening a dataset and discovering some information about the file\n",
    "\n",
    "**Note 1:** This does not yet load/extract any data, just opens the file.\n",
    "\n",
    "**Note 2:** RasterIO, like GDAL, is perfectly happy with virtual file systems such as zip files, memory bufers, streaming data or data hosted on a HTTP or FTP server. Here we are connecting to an 18GB file on a HTTP file system so we use the **/vsicurl/** driver\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The paths to the South Australian Seasonal Landsat Mosaics on the Queensland RDSI Node\n",
    "refDataPath = '/vsicurl/http://qld.auscover.org.au/public/data/landsat\\\n",
    "/surface_reflectance/sa/l8olre_sa_m201412201502_dbia2.tif'\n",
    "\n",
    "# Open the terrain corrected surface reflectance data\n",
    "refDataSet = rasterio.open(refDataPath)\n",
    "\n",
    "# The dataset object contains metadata about the raster data\n",
    "print \"Bands:\\t \", refDataSet.count\n",
    "print \"Height:\\t \", refDataSet.height\n",
    "print \"Width:\\t \", refDataSet.width\n",
    "print \"BoundingBox: \", refDataSet.bounds\n",
    "print \"CRS:\\t \", refDataSet.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A note on converting between real world and transformed coordinates\n",
    " - Most of the time you'll need to work out where the image pixel corresponding to a real world coordinate is.\n",
    " - That means working out how to do coordinate transformations.\n",
    " - Here we use the **affine** transformation attached to the dataset object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Some dummy coordinates as an example\n",
    "col, row = 0, 0\n",
    "easting, northing = (-460755.0, -2716865.0)\n",
    "\n",
    "# Convert from image to world coordinates\n",
    "print \"Forward Transformation:\\n\", refDataSet.affine\n",
    "print \"\\nConvert (%s %s) to Coordinates:\\t\" % (col, row) ,refDataSet.affine * (col, row)\n",
    "\n",
    "# Convert from world to image coordinates\n",
    "print \"\\nInverse Transformation:\\n\", ~refDataSet.affine\n",
    "print \"\\nConvert (%s %s) to Row and Column:\\t\" % (easting, northing) ,~refDataSet.affine * (easting, northing)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading in a data subset from the remote file system\n",
    "\n",
    " - First we compute the area we want to extract from the image\n",
    " - Then we pull the data subset from the server\n",
    " - In all these examples, **subsets are extracted without having to read the entire file** so you only pull the data you are working with back to your browser. In the case of this example, the (compressed) file is 18GB in size but we just read the portion we want so the subset is available in under a second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This is a subset over Port Augusta, SA with the bounding box in EPSG:3577 (Australian Albers)\n",
    "topLeft = (527643.414,-3541553.603)\n",
    "bottomRight = (574184.998,-3565421.083)\n",
    "\n",
    "# Compute the image coordinates from the real world coordinates using an inverse transform\n",
    "y1, x1 =  ~refDataSet.affine * topLeft\n",
    "y2, x2 =  ~refDataSet.affine * bottomRight\n",
    "\n",
    "# Read in the SWIR, NIR and Red bands within the specified bounding box\n",
    "refDataSubset = refDataSet.read([5,4,3], window=((x1,x2), (y1,y2)))\n",
    "\n",
    "# Print some information on the subset\n",
    "print \"Data extract shape: {0}, dtype: {1}\".format(refDataSubset.shape, refDataSubset.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Displaying a RGB Image\n",
    "\n",
    " - Like any remote sensing package, we have to stretch the data into 8 bit space and set up an informative band combination.\n",
    " - For more info on common band combinations: http://landsat.usgs.gov/L8_band_combos.php\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Stretch and scale the data from 0 to 255 as an unsigned 8 bit integer for display\n",
    "refDataSubsetScaled=numpy.clip(refDataSubset / 5000.0 * 255.0, 0, 255).astype('uint8')\n",
    "\n",
    "# Plot the image using matplotlib\n",
    "plt.figure(figsize=(16,6))\n",
    "\n",
    "# Rearrange the axes so that it works with matplotlib that likes BIP not BSQ data\n",
    "plt.imshow(numpy.rollaxis(refDataSubsetScaled,0,3))\n",
    "\n",
    "# Add a title to the plot\n",
    "plt.title('Port Augusta Area - Landsat Surface Reflectance (Bands 5,4,3)', fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing NDVI from the extracted Subset\n",
    "\n",
    " - Everyone uses NDVI as a 'Hello World\" remote sensing example so I will as well.\n",
    " - Note that these data use 32767 as a nodata value so we need to mask these out from the calculation.\n",
    " - We also embed a simple histogram of the NDVI to show the spread of values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Extract the bands from the nD array\n",
    "b4=refDataSubset[1].astype(float)\n",
    "b3=refDataSubset[2].astype(float)\n",
    "\n",
    "# Allow division by zero\n",
    "numpy.seterr(divide='ignore', invalid='ignore')\n",
    "\n",
    "# Calculate NDVI\n",
    "ndvi = (b4 - b3) / (b4 + b3)\n",
    "\n",
    "# Make the NoData values NaN so they get masked from the analysis\n",
    "ndvi[b4 == 32767] = numpy.nan\n",
    "\n",
    "# Plot the result using matplotlib\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "\n",
    "# Select an apprppriate colourmap\n",
    "plt.imshow(ndvi,cmap='brg')\n",
    "\n",
    "# Display a colorbar\n",
    "plt.colorbar()\n",
    "\n",
    "# Add a title to the plot\n",
    "plt.title('Port Augusta Area - Landsat NDVI Map and Histogram', fontsize=20)\n",
    "\n",
    "# Inset a histogram of the NDVI values\n",
    "a = plt.axes([.62, .65, .10, .20])\n",
    "\n",
    "# Compute the histogram\n",
    "n, bins, patches = plt.hist(ndvi.ravel(), bins=200, range=(-1, 1),facecolor='grey',histtype='stepfilled')\n",
    "\n",
    "# Add a title to the histogram plot\n",
    "plt.title('Distribution', fontsize=16, color='black')\n",
    "\n",
    "# Remove the Y axis labels for clarity\n",
    "a.yaxis.set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### Clipping that file locally for additional analysis\n",
    " - The [GDAL utilities](http://www.gdal.org/gdal_utilities.html) are perfect for this type of application\n",
    " - First we'll use [gdal_translate](http://www.gdal.org/gdal_translate.html) on the command line to clip out this same subset from the 18GB parent file on the server, resample it to 250m by pixel averaging (let's say we want to compare to MODIS) and change the format of the file to [KEA](http://kealib.org/)\n",
    " - Then we'll check the file information using [gdalinfo](http://www.gdal.org/gdalinfo.html)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Run the gdal_translate utility\n",
    "! gdal_translate -of KEA -projwin 527643 -3541553 574184 -3565421 -tr 250 250 -r average \\\n",
    "'/vsicurl/http://qld.auscover.org.au/public/data/landsat/surface_reflectance/sa/l8olre_sa_m201412201502_dbia2.tif' \\\n",
    "l8olre_portaugusta_m201412201502_dbia2.kea"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Run the gdalinfo utility\n",
    "! gdalinfo l8olre_portaugusta_m201412201502_dbia2.kea"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "***\n",
    "### Run an analysis of fractional cover data over the same area\n",
    " - These data contain values representing the percentage of a pixel that is covered by green vegetation, dead vegetation and bare ground. \n",
    " - A common task is to ascertain how much bare ground above a certain threshold exists in an area\n",
    " - For this example we simply repeat the extraction process again using the corresponding fractional cover data set\n",
    " - It's then a simple matter to plot a **cumulative histogram** of bare ground values for the reporting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Open the fractional cover data on the server\n",
    "fcDataPath = '/vsicurl/http://qld.auscover.org.au/public/data/landsat\\\n",
    "/seasonal_fractional_cover/fractional_cover/sa/lztmre_sa_m201412201502_dima2.tif'\n",
    "fcDataSet = rasterio.open(fcDataPath)\n",
    "\n",
    "# Extract the first band (Bare Ground) from the data on the Auscover server\n",
    "fcDataSubset = fcDataSet.read(1, window=((x1,x2), (y1,y2)))\n",
    "\n",
    "# Print some information on the extracted data\n",
    "print \"Data extract shape: {0}, dtype: {1}\".format(fcDataSubset.shape, fcDataSubset.dtype)\n",
    "\n",
    "# Scale the bare ground data as percent from 0 to 100\n",
    "barePercent = numpy.clip(fcDataSubset - 100.0, 0, 100)\n",
    "\n",
    "# Setup the matplotlib figure\n",
    "fig = plt.figure(figsize=(9,5))\n",
    "\n",
    "# Compute the cumulative histogram and plot\n",
    "n, bins, patches = plt.hist\\\n",
    "    (barePercent.ravel(), bins=100, range=(0, 100),color='brown',\\\n",
    "     cumulative=True, normed=1,histtype='step',linewidth=3)\n",
    "\n",
    "# Add in labels and grids to the histogram\n",
    "plt.xlabel('Bare Ground')\n",
    "plt.ylabel('Percent Area')\n",
    "plt.title('Cumulative Distribution of Bare Ground', fontsize=20)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lets have a look at the bare ground image to locate those highly bare areas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plot the result using matplotlib\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "\n",
    "# Select an apprppriate colourmap\n",
    "plt.imshow(barePercent,cmap='gist_earth')\n",
    "\n",
    "# Display a colorbar\n",
    "plt.colorbar()\n",
    "\n",
    "# Add a title to the plot\n",
    "plt.title('Port Augusta Area - Landsat Bare Ground', fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "##  Time series analysis example\n",
    "**We'll look at the Spyglass research property in Queensland**\n",
    "\n",
    " - This requires a catalog service or a good filename convention\n",
    " - We need to:\n",
    "   - work out what image data we need\n",
    "   - open the subsets\n",
    "   - calculate statistics for each subset\n",
    "   - Plot the results.\n",
    " - This can also be slow over the web depending on your internet connection speed and the size of the subset due to the face we download over 100 subsets.\n",
    " - It's still faster than downloading 2TB of data then trying to analyise that on your desktop!\n",
    "\n",
    "\n",
    "\n",
    "_Aside: This image is embedded using Geoserver to render on the fly_\n",
    "<img src='http://qld.auscover.org.au/geoserver/aus/wms?&time=2015-12&service=WMS&version=1.1.0&request=GetMap&layers=aus:fractional_cover&styles=&bbox=1420000,-2170000,1450000,-2150000&width=300&height=300&srs=EPSG:3577&format=image%2Fjpeg'>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "# Import some additional python packages needed for this example\n",
    "from xml.etree import ElementTree\n",
    "from IPython.display import display, HTML\n",
    "from scipy.interpolate import interp1d\n",
    "\n",
    "# Download and parse the THREDDS catalog XML to find all the seasonal fractional cover images in Queensland.\n",
    "catalogUrl = 'http://qld.auscover.org.au/thredds/catalog/auscover/landsat\\\n",
    "/seasonal_fractional_cover/fractional_cover/qld/catalog.xml'\n",
    "root = ElementTree.XML(requests.get(catalogUrl).text.encode('utf-8'))\n",
    "\n",
    "# Build a list of datasets in the catalog XML\n",
    "datasets=root.find(\"{http://www.unidata.ucar.edu/namespaces/thredds/InvCatalog/v1.0}dataset\")\n",
    "\n",
    "# Define the property bounding box in EPSG:3577 (Australian Albers)\n",
    "# This corresponds to the subset displayed above\n",
    "paddockTopLeft = (1420000,-2150000)\n",
    "paddockBottomRight = (1450000,-2170000)\n",
    "\n",
    "# Print some statistics about property size\n",
    "print \"Analysing % 0.0f Ha of Imagery\" % \\\n",
    "    ((paddockBottomRight[0]-paddockTopLeft[0]) * (paddockTopLeft[1]-paddockBottomRight[1]) / 10000)            \n",
    "\n",
    "\n",
    "# Setup empty lists to hold the property statistics\n",
    "seasonalDates= []\n",
    "seasonalBareGroundMedian = []\n",
    "seasonalBareGroundMean = []\n",
    "seasonalBareGroundStd = []\n",
    "seasonalBareGroundMin = []\n",
    "seasonalBareGroundMax = []\n",
    "\n",
    "# Loop through all the datasets found in the catalog\n",
    "for dataset in datasets:\n",
    "    \n",
    "    # Get the name of the file\n",
    "    fileName = dataset.get('name')\n",
    "    if fileName is not None:\n",
    "        # Open the dataset on the server using RasterIO\n",
    "        fileUrl =  \"/vsicurl/http://qld.auscover.org.au/thredds/fileServer/\" + str(dataset.get('urlPath'))\n",
    "        tsDataSet = rasterio.open(fileUrl)\n",
    "\n",
    "        # Compute the image coordinates from the real world coordinates using the inverse transform\n",
    "        y3, x3 =  ~tsDataSet.affine * paddockTopLeft\n",
    "        y4, x4 =  ~tsDataSet.affine * paddockBottomRight\n",
    "        \n",
    "        # Read the property subset into a numpy array\n",
    "        tsDataSubset = tsDataSet.read(1, window=((x3,x4), (y3,y4)))\n",
    "        \n",
    "        # Remove the nodata values (0) from the flattened array and rescale from 0 to 100\n",
    "        tsDataSubset = numpy.ravel(tsDataSubset[tsDataSubset <> 0] - 100)\n",
    "        \n",
    "        # Check that we have some data available for the time period\n",
    "        if tsDataSubset.size > 1:\n",
    "            \n",
    "            # Compute some statistics for the region and add to the lists\n",
    "            # Median\n",
    "            seasonalBareGroundMedian.append(numpy.median(tsDataSubset))\n",
    "            # Mean\n",
    "            seasonalBareGroundMean.append(numpy.mean(tsDataSubset))\n",
    "            # Standard Deviation\n",
    "            seasonalBareGroundStd.append(numpy.std(tsDataSubset))\n",
    "            #  5th percentile, a robust minimum\n",
    "            seasonalBareGroundMin.append(numpy.percentile(tsDataSubset,5))\n",
    "            # 95th percentile, a robust maximum\n",
    "            seasonalBareGroundMax.append(numpy.percentile(tsDataSubset,95))\n",
    "\n",
    "            # Lazy date calculation. Should really be using datetime here!\n",
    "            \n",
    "            # Starting year from the filename\n",
    "            year = float(fileName.split('_')[2][1:5])\n",
    "            # Starting month from the filename\n",
    "            month = float(fileName.split('_')[2][5:7])\n",
    "            # Add the fractional year to the data list\n",
    "            seasonalDates.append(year + month / 12)\n",
    "\n",
    "\n",
    "# Print a handy table of all of our statistics\n",
    "propertyStats = [seasonalDates,seasonalBareGroundMin,seasonalBareGroundMedian,\\\n",
    "                 seasonalBareGroundMean,seasonalBareGroundMax,seasonalBareGroundStd]\n",
    "\n",
    "# This hack inserts a formatted HTML table inline into our notebook\n",
    "display(HTML('<table><tr>{}</tr></table>'.format('</tr><tr>'.join(\n",
    "            '<td>{}</td>'.format('</td><td>'.join(str(_) for _ in row)) for row in propertyStats))))\n",
    "\n",
    "\n",
    "\n",
    "# Instead of just plotting the points, lets try some interpolation to make a smooth curve\n",
    "\n",
    "# Make a linearly spaced array of fractional years spanning the range of our data\n",
    "interpDates = numpy.linspace(min(seasonalDates), max(seasonalDates), num=500, endpoint=True)\n",
    "\n",
    "# Build a function to do linear interpolation of the maximum bare ground\n",
    "linearFunction = interp1d(seasonalDates, seasonalBareGroundMax, kind='linear')\n",
    "\n",
    "# Build a function to do cubic interpolation of the maximum bare ground\n",
    "cubicFunction = interp1d(seasonalDates, seasonalBareGroundMax, kind='cubic')\n",
    "\n",
    "# Fit a polynomial to the maximum bare ground\n",
    "fittedPolynomial = numpy.poly1d(numpy.polyfit(seasonalDates, seasonalBareGroundMax, 2))\n",
    "\n",
    "# Setup the figure\n",
    "fig = plt.figure(figsize=(16,8))\n",
    "\n",
    "# Plot the points, the linear interpolation and the cubic interpolation\n",
    "plt.plot(seasonalDates, seasonalBareGroundMax, 'o',\\\n",
    "         interpDates, cubicFunction(interpDates),'-',\\\n",
    "         interpDates, linearFunction(interpDates), '--',\\\n",
    "         seasonalDates, fittedPolynomial(seasonalDates),'.')\n",
    "\n",
    "# Add annotation to the plot\n",
    "plt.legend(['Bare Ground Observations', 'Cubic interpolation', 'Linear interpolation', 'Fitted Polynomial'], loc='best')\n",
    "plt.xlabel('Season Start Date')\n",
    "plt.ylabel('Bare Ground (%)')\n",
    "plt.title('Maximum (95th percentile) Bare Ground across Spyglass', fontsize=16)\n",
    "plt.grid(True)\n",
    "plt.ticklabel_format(useOffset=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# Aside - Accessing field data via services is possible as well!\n",
    "  - Several ways using WMS and WCS\n",
    "  - This super simple method used WCS geojson and the Folium library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Import folium, a leaflet based python \n",
    "import folium\n",
    "\n",
    "# The URL from Geoserver for the GeoJSON TLS Sites\n",
    "geojsonUrl = 'http://qld.auscover.org.au/geoserver/aus/ows?\\\n",
    "service=WFS&version=1.0.0&request=GetFeature&\\\n",
    "typeName=aus:tls_sites_display&maxFeatures=100&outputFormat=application%2Fjson'\n",
    "\n",
    "# GET the data from the server as text\n",
    "tlsSites = requests.get(geojsonUrl).text\n",
    "\n",
    "# Build the base Map object\n",
    "tlsMap = folium.Map(tiles='OpenStreetMap',location=[-27, 135], zoom_start=4)\n",
    "\n",
    "# Add the GeoJSON from the server\n",
    "tlsMap.choropleth(geo_str=tlsSites)\n",
    "\n",
    "# Show the map\n",
    "tlsMap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# NCI Data Access: Python NetCDF Landsat8\n",
    "\n",
    "**The following will go through how to:** <br \\>\n",
    "   1. Access published netCDF data through NCI's THREDDS Data Server (using OPeNDAP)\n",
    "   2. Extract/view data\n",
    "   3. Save data subset to new file\n",
    "\n",
    "\n",
    "### Import python libraries\n",
    "\n",
    "There are several Python libraries available to work with netCDF and HDF file formats. This tutorial will use `netCDF4` but others, such as `h5py`, `cdms2`, and `gdal` can also be used. For more information on these other libraries, please see the main tutorial page. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy\n",
    "from netCDF4 import Dataset\n",
    "import matplotlib.pyplot as plt \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Open/read file\n",
    "**Note:** This does not yet load/extract any data, just opens the file.\n",
    "\n",
    "### The 'Dataset' function is used to open a file with Python's netCDF4 library. \n",
    "For local files, this will be the filepath (i.e., /g/data...) while for remote access, this will be the OPeNDAP data URL. For instructions on how to find the OPeNDAP URL, please see: [THREDDS Data Access](https://nbviewer.jupyter.org/github/kdruken/Notebooks/blob/master/THREDDS_DataAccess.ipynb)\n",
    "\n",
    "#### Accessing data remotely (OPeNDAP)\n",
    "#### After opening the file with the OPeNDAP address, the file can be handled in the same manner as a local file. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "url = 'http://dapds00.nci.org.au/thredds/dodsC/rs0/tiles/EPSG3577/LS8_OLI_TIRS_NBAR/LS8_OLI_TIRS_NBAR_3577_-10_-27_2013.nc'\n",
    "f = Dataset(url, 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Browse information about the file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### File dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for item in f.dimensions:\n",
    "    print f.dimensions[item].name, f.dimensions[item].size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### File variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vars = f.variables.keys()\n",
    "for item in vars:\n",
    "    print 'Variable: \\t', item\n",
    "    print 'Dimensions: \\t', f[item].dimensions\n",
    "    print 'Shape:    \\t', f[item].shape, '\\n'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracting data (using index values)\n",
    "A nice feature of netCDF/HDF file formats is that you can extract subsets without having to read the entire file (or variable). The example below demonstrates the simplest subsetting example by directly specifying the subset indices. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Read variables (but not yet extract)\n",
    "band2 = f['band_2']\n",
    "y = f['y']\n",
    "x = f['x']\n",
    "t = f['time']\n",
    "\n",
    "# Subset indices\n",
    "x1, x2 = 1000,3999\n",
    "y1, y2 = 0,3000\n",
    "t1 = 9\n",
    "\n",
    "# Extract\n",
    "band2_subset = band2[t1, y1:y2, x1:x2]\n",
    "y_subset = y[y1:y2]\n",
    "x_subset = x[x1:x2]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set figure size\n",
    "plt.figure(figsize=(12,6))\n",
    "\n",
    "# Plot data subset with equal axes and colorbar\n",
    "plt.contourf(x_subset, y_subset, band2_subset)\n",
    "plt.axis('equal')\n",
    "cbar = plt.colorbar()\n",
    "\n",
    "# Add figure title and labels\n",
    "# We can make use of the defined variable attributes to do this\n",
    "plt.title(band2.long_name+'\\n', fontsize=18)\n",
    "plt.xlabel(x.long_name+' ('+x.units+') ', fontsize=16)\n",
    "plt.ylabel(y.long_name+' ('+y.units+') ', fontsize=16)\n",
    "\n",
    "# Adjust tick mark size\n",
    "cbar.ax.tick_params(labelsize=16) \n",
    "plt.tick_params(labelsize=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot subset as RGB image\n",
    "For more info on common band combinations: http://landsat.usgs.gov/L8_band_combos.php\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Read in bands\n",
    "band4_subset = f['band_4'][t1, y1:y2, x1:x2]\n",
    "band6_subset = f['band_6'][t1, y1:y2, x1:x2]\n",
    "band7_subset = f['band_7'][t1, y1:y2, x1:x2]\n",
    "\n",
    "# Bands must be clipped (value of 6000 was chosen in this case) and scaled to values between (0, 255) to plot as RGB image.\n",
    "b4 = band4_subset.clip(0, 6000) / 6000. * 255\n",
    "b6 = band6_subset.clip(0, 6000) / 6000. * 255\n",
    "b7 = band7_subset.clip(0, 6000) / 6000. * 255\n",
    "\n",
    "## Combine the bands of interest into numpy NxNx3 dimensional array\n",
    "# Note: The data type must be converted to 'uint8' to plot as image\n",
    "rgb = numpy.stack((b7, b6, b4), axis=2).astype('uint8')\n",
    "print \"New array shape: {0}, dtype: {1}\".format(rgb.shape, rgb.dtype)\n",
    "\n",
    "# Set figure size\n",
    "plt.figure(figsize=(12,12))\n",
    "\n",
    "# Plot image\n",
    "plt.imshow(rgb, extent=[x_subset[0], y_subset[-1], x_subset[-1], y_subset[0]])\n",
    "\n",
    "# Add figure title and labels\n",
    "# We can make use of the defined variable attributes to do this\n",
    "plt.title('Landsat 8 False Colour: Bands (7, 6, 4) \\n', fontsize=20)\n",
    "plt.xlabel(x.long_name+' ('+x.units+') ', fontsize=16)\n",
    "plt.ylabel(y.long_name+' ('+y.units+') ', fontsize=16)\n",
    "\n",
    "\n",
    "# Adjust tick mark size\n",
    "plt.tick_params(labelsize=16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Close the file\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# Using the  NCI THREDDS Data Server with Siphon\n",
    "\n",
    "Siphon is a collection of Python utilities for downloading data from Unidata data technologies. More information on installing and using Unidata's Siphon can be found: \n",
    "https://github.com/Unidata/siphon\n",
    "\n",
    "**The following will go through how to:**\n",
    "- Use Siphon to find and query Landsat data hosted at NCI\n",
    "- Use Siphon to find and query MODIS data hosted on NCI THREDDS Data Server\n",
    "\n",
    "\n",
    "## Browse for data\n",
    "\n",
    "Begin by going to NCI's Geonetwork page: http://geonetwork.nci.org.au/\n",
    "\n",
    "This page contains the metadata records for NCI Data Collections as well as information on where to find the data. \n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/nci/nci-notebooks/master/Data_Access/Using_Siphon/images/gn1.png\">\n",
    "\n",
    "In this example, we will search for Landsat data:\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/nci/nci-notebooks/master/Data_Access/Using_Siphon/images/gn2.png\">\n",
    "\n",
    "If we click on the first result, we see a brief overview of the metadata record. **Note:** For the full record, navigate to the upper-right corner of your browser to change to the \"Full view\" (eyeball icon). \n",
    "\n",
    "One of the options under **Download and links** is the NCI THREDDS Data Server Catalog page: \n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/nci/nci-notebooks/master/Data_Access/Using_Siphon/images/gn3.png\">\n",
    "\n",
    "By navigating to this link, the available (public) data subcollections and datasets will be visible:\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/nci/nci-notebooks/master/Data_Access/Using_Siphon/images/thredds10.png\">\n",
    "\n",
    "\n",
    "In this example, let's navigate to the ** LANDSAT data: scenes and tiles/ tiles/ EPSG3577/ LS8_OLI_TIRS_NBAR/ ** dataset: \n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/nci/nci-notebooks/master/Data_Access/Using_Siphon/images/thredds11.png\">\n",
    "\n",
    "## Using Siphon\n",
    "\n",
    "Once selecting a parent dataset directory, Siphon can be used to search and use the data access methods and services provided by THREDDS. For example, Siphon will return a list of data endpoints for the OPeNDAP data URL, NetCDF Subset Service (NCSS), Web Map Service (WMS), Web Coverage Service (WCS), and the HTTP link for direct download. \n",
    "\n",
    "### Import python packages\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from netCDF4 import Dataset \n",
    "from siphon import catalog, ncss\n",
    "import datetime\n",
    "import matplotlib.pyplot as plt "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Provide the top-level URL from the THREDDS page above\n",
    "Note: You can leave the '.html' but you will receive a message saying it was changed to '.xml'. \n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/nci/nci-notebooks/master/Data_Access/Using_Siphon/images/thredds13.png\">\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "catalog_url = 'http://dapds00.nci.org.au/thredds/catalog/rs0/tiles/EPSG3577/LS8_OLI_TIRS_NBAR/catalog.xml'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now we use Siphon to list all the available datasets under this catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tds = catalog.TDSCatalog(catalog_url)\n",
    "datasets = list(tds.datasets)\n",
    "endpts = tds.datasets.values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Some of the datasets...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "datasets[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### And their associated endpoints for data services:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for key, value in endpts[0].access_urls.items():\n",
    "    print key, value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now we can use Siphon along with some form of query method to find some data\n",
    "\n",
    "#### This example will use Shapely to find intersecting Polygon shapes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from shapely.geometry import Polygon\n",
    "from shapely.wkt import loads\n",
    "\n",
    "query = (136, 138, -29.3, -27.8)\n",
    "query = Polygon([[136, -29.3], [136, -27.8], [138, -27.8], [138, -29.3]])\n",
    "\n",
    "# What this query looks like in WKT\n",
    "print query.wkt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loop through the datasets and check if the Landsat's geospatial bounds (which is in a WKT polygon format) intersects with the query\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time \n",
    "\n",
    "matches = []\n",
    "for dataset in endpts:\n",
    "    dap = dataset.access_urls['OPENDAP']\n",
    "    with Dataset(dap) as f:\n",
    "        bounds = loads(f.geospatial_bounds.encode())\n",
    "        if bounds.intersects(query) == True:\n",
    "            print dap\n",
    "            matches.append(dap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's take a quick look at what was found \n",
    "\n",
    "(Because we are accessing data remotely through OPeNDAP, let's look at a lower resolution so it doesn't exceed memory limits.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,8))\n",
    "\n",
    "for match in matches:\n",
    "    with Dataset(match) as f:\n",
    "        x = f.variables['x'][::50]\n",
    "        y = f.variables['y'][::50]\n",
    "        t = f.variables['time'][::5]\n",
    "        \n",
    "        for i in range(0, len(t)):\n",
    "            b2 = f.variables['band_2'][i,::50,::50]\n",
    "            plt.pcolormesh(x, y, b2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Now we'll use the same process to look at a MODIS time series\n",
    "\n",
    "#### Start by defining the parent catalog URL from NCI's THREDDS Data Server\n",
    "#### Then use Siphon to explore the available datasets and data services end points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "url = 'http://dapds00.nci.org.au/thredds/catalog/u39/public/data/modis/fractionalcover-clw/v2.2/netcdf/catalog.xml'\n",
    "\n",
    "tds = catalog.TDSCatalog(url)\n",
    "datasets = list(tds.datasets)\n",
    "endpts = tds.datasets.values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We can create a small function that uses Siphon's Netcdf Subset Service (NCSS) to extract a spatial request (defined by a lat/lon box)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_data(dataset, bbox):    \n",
    "    nc = ncss.NCSS(dataset.access_urls['NetcdfSubset'])\n",
    "    query = nc.query()\n",
    "    query.lonlat_box(north=bbox[3],south=bbox[2],east=bbox[1],west=bbox[0])\n",
    "    query.variables('bs')\n",
    "    \n",
    "    data = nc.get_data(query)\n",
    "    \n",
    "    lon = data['longitude'][:]\n",
    "    lat = data['latitude'][:]\n",
    "    bs = data['bs'][0,:,:]\n",
    "    t = data['time'][:]\n",
    "    \n",
    "    time_base = datetime.date(year=1800, month=1, day=1)\n",
    "    time = time_base + datetime.timedelta(t[0])\n",
    "    \n",
    "    return lon, lat, bs, time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test the function on a single file and view result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bbox = (135, 140, -31, -27)\n",
    "lon, lat, bs, t = get_data(endpts[0], bbox)\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(bs, extent=bbox, cmap='gist_earth', origin='upper')\n",
    "\n",
    "plt.xlabel('longitude (degrees)', fontsize=14)\n",
    "plt.ylabel('latitude (degrees)', fontsize=14)\n",
    "print \"Date: \", t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loop and query over the collection and save every plot as a PNG image file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bbox = (135, 140, -31, -27)\n",
    "plt.figure(figsize=(10,10))\n",
    "\n",
    "for endpt in endpts[:15]:\n",
    "    try:\n",
    "        lon, lat, bs, t = get_data(endpt, bbox)\n",
    "\n",
    "        plt.imshow(bs, extent=bbox, cmap='gist_earth', origin='upper')\n",
    "        plt.clim(vmin=-2, vmax=100)\n",
    "\n",
    "        plt.tick_params(labelsize=14)\n",
    "        plt.xlabel('longitude (degrees)', fontsize=14)\n",
    "        plt.ylabel('latitude (degrees)', fontsize=14)\n",
    "\n",
    "        plt.title(\"Date: \"+str(t), fontsize=16, weight='bold')\n",
    "        plt.savefig(endpt.name+\".png\")\n",
    "        plt.cla()\n",
    "    except:\n",
    "        pass\n",
    "\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Convert the series of PNG files above into a GIF\n",
    "Uses Imagemagick convert on the command line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Run the command line conversion\n",
    "! convert -fuzz 50% -layers optimize -delay 50 -loop 0  *.png animated.gif\n",
    "# Delete the individial PNG files after creation of the GIF\n",
    "! rm *.png"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Show the animation of the temporal imagery\n",
    "<img src=\"animated.gif\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We can also use Siphon to extract a single point by creating another function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_point(dataset, lat, lon):\n",
    "    nc = ncss.NCSS(dataset.access_urls['NetcdfSubset'])\n",
    "    query = nc.query()\n",
    "    query.lonlat_point(lon, lat)\n",
    "    query.variables('bs')\n",
    "    \n",
    "    data = nc.get_data(query)\n",
    "    bs = data['bs'][0]\n",
    "    date = data['date'][0]\n",
    "    \n",
    "    return bs, date"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test this function on a point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bs, date = get_point(endpts[4], -27.75, 137)\n",
    "print bs, date"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time series example\n",
    " - We use our function to drill the data at the point for each timestep\n",
    " - Then we can show the timeseries on a plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = []\n",
    "for endpt in endpts[::20]:\n",
    "    bs, date = get_point(endpt, -27.75, 137)\n",
    "    data.append([date, bs])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sortOrder = numpy.argsort(numpy.array(data)[:,0])\n",
    "BS = numpy.array(data)[sortOrder,1]\n",
    "Date = numpy.array(data)[sortOrder,0]\n",
    "\n",
    "plt.figure(figsize=(12,6))\n",
    "plt.plot(Date, BS, '-o', linewidth=2, markersize=8)\n",
    "\n",
    "plt.tick_params(labelsize=14)\n",
    "plt.xlabel('date', fontsize=14)\n",
    "plt.ylabel('fractional cover of bare soil (%)', fontsize=14)\n",
    "plt.title('Lat, Lon: -27.75, 137', fontsize=16)"
   ]
  },
  {
   "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": 0
}