3 - Example Time Series - Landsat.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by initializing the Earth Engine API."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import ee\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we import in other Python packages we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import datetime\n",
    "from matplotlib import dates\n",
    "import matplotlib.dates as mdates\n",
    "from pylab import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define the Landsat bands that we would like to plot, along with the starting and ending times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xBand = 'time'\n",
    "yBandList = [\n",
    "        'B1',\n",
    "        u'B2',\n",
    "        u'B3',\n",
    "        u'B4',\n",
    "        u'B5',\n",
    "        u'B6_VCID_1',\n",
    "        u'B6_VCID_2',\n",
    "        u'B7',\n",
    "        u'B8',\n",
    "    ]\n",
    "startTime = datetime.datetime(2000, 1, 1)\n",
    "endTime = datetime.datetime(2004, 1, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we contruct a filtered ImageCollection for the date range, and extract band information for a specified point location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "collection = ee.ImageCollection('LE7_L1T').filterDate(startTime, endTime)\n",
    "point = {'type':'Point', 'coordinates':[ -116.88629,36.56122]};  # death valley (should be stable)\n",
    "info = collection.getRegion(point,500).getInfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We separate the information returned into column headers and data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# extract the header column names\n",
    "header = info[0]\n",
    "# create a Numpy array of the data\n",
    "data = array(info[1:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we extract time information and convert it to at Python datatime data type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# extract the time information\n",
    "iTime = header.index('time')\n",
    "# convert to Python datetime objects\n",
    "time = [datetime.datetime.fromtimestamp(i/1000) for i in (data[0:,iTime].astype(int))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the data columns what we want to display on the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "iBands = [header.index(b) for b in yBandList]\n",
    "yData = data[0:,iBands].astype(np.float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate NDVI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "band3 = yData[:,2]\n",
    "band4 = yData[:,3]\n",
    "ndvi = (band4 - band3) / (band4 + band3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally, we create a plot of MODIS band values as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# matplotlib date format object\n",
    "\n",
    "fig = figure(figsize=(12,8), dpi=80)\n",
    "\n",
    "# plot the band values\n",
    "ax1 = fig.add_subplot(211)\n",
    "ax1.plot(time, yData[:,2], 'o', color=\"red\", label=\"Band 3\")\n",
    "ax1.plot(time, yData[:,3], 'o', color=\"magenta\",  label=\"Band 4\")\n",
    "ax1.legend(loc='best')\n",
    "ax1.grid(True)\n",
    "\n",
    "#plt.title('Band values as a function of time')\n",
    "ax1.set_ylabel('Band Values')\n",
    "\n",
    "# plot NDVI\n",
    "ax2 = fig.add_subplot(212, sharex=ax1)\n",
    "ax2.plot(time, ndvi, 'o', color=\"black\", label=\"NDVI\")\n",
    "ax2.grid(True)\n",
    "start, end = ax2.get_xlim()\n",
    "ax2.xaxis.set_ticks(np.arange(start, end, 64.5))\n",
    "\n",
    "# Format the ticks.\n",
    "years    = mdates.YearLocator()   # every year\n",
    "months   = mdates.MonthLocator()  # every month\n",
    "yearsFmt = mdates.DateFormatter('%Y')\n",
    "\n",
    "ax2.set_ylabel('NDVI')\n",
    "\n",
    "ax2.xaxis.set_major_locator(years)\n",
    "ax2.xaxis.set_major_formatter(yearsFmt)\n",
    "ax2.xaxis.set_minor_locator(months)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Convert the timestamp to a numpy array\n",
    "t = np.array([i.toordinal() for i in time])\n",
    "t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{bmatrix}\n",
    "t_1 & 1 \\\\\\\n",
    "t_2 & 1 \\\\\\\n",
    "\\vdots & \\vdots \\\\\\\n",
    "t_n & 1\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "x_0 \\\\\\\n",
    "x_1 \n",
    "\\end{bmatrix}=\n",
    "\\begin{bmatrix}\n",
    "y_1 \\\\\\\n",
    "y_2 \\\\\\\n",
    "\\vdots \\\\\\\n",
    "y_n\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "$$ \\mathbf{A}\\mathbf{x} = \\mathbf{b} $$\n",
    "\n",
    "$$ \\mathbf{x} = \\mathbf{A} \\backslash \\mathbf{b} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = array([ t, ones(len(t))]).transpose()\n",
    "b = ndvi\n",
    "x = linalg.lstsq(A,b)[0] # obtaining the parameters\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "b_hat = A.dot(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig2 = figure(figsize=(12,4), dpi=80)\n",
    "plot(time,b_hat.transpose(),'r-',t,b,'o')\n",
    "fig2.fmt_xdata = mdates.DateFormatter('%Y-%m-%d')\n",
    "fig2.autofmt_xdate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "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
}