

In [ ]:

    
%matplotlib inline

We start by initializing the Earth Engine API.



In [ ]:

    
import ee
ee.Initialize()

Then we import in other Python packages we need.



In [ ]:

    
import datetime
from matplotlib import dates
import matplotlib.dates as mdates
from pylab import *

Next we define the Landsat bands that we would like to plot, along with the starting and ending times.



In [ ]:

    
xBand = 'time'
yBandList = [
        'B1',
        u'B2',
        u'B3',
        u'B4',
        u'B5',
        u'B6_VCID_1',
        u'B6_VCID_2',
        u'B7',
        u'B8',
    ]
startTime = datetime.datetime(2000, 1, 1)
endTime = datetime.datetime(2004, 1, 1)

Next we contruct a filtered ImageCollection for the date range, and extract band information for a specified point location.



In [ ]:

    
collection = ee.ImageCollection('LE7_L1T').filterDate(startTime, endTime)
point = {'type':'Point', 'coordinates':[ -116.88629,36.56122]};  # death valley (should be stable)
info = collection.getRegion(point,500).getInfo()

We separate the information returned into column headers and data.



In [ ]:

    
# extract the header column names
header = info[0]
# create a Numpy array of the data
data = array(info[1:])

Next we extract time information and convert it to at Python datatime data type.



In [ ]:

    
# extract the time information
iTime = header.index('time')
# convert to Python datetime objects
time = [datetime.datetime.fromtimestamp(i/1000) for i in (data[0:,iTime].astype(int))]

Extract the data columns what we want to display on the plot.



In [ ]:

    
iBands = [header.index(b) for b in yBandList]
yData = data[0:,iBands].astype(np.float)

Calculate NDVI



In [ ]:

    
band3 = yData[:,2]
band4 = yData[:,3]
ndvi = (band4 - band3) / (band4 + band3)

And finally, we create a plot of MODIS band values as a function of time.



In [ ]:

    
# matplotlib date format object

fig = figure(figsize=(12,8), dpi=80)

# plot the band values
ax1 = fig.add_subplot(211)
ax1.plot(time, yData[:,2], 'o', color="red", label="Band 3")
ax1.plot(time, yData[:,3], 'o', color="magenta",  label="Band 4")
ax1.legend(loc='best')
ax1.grid(True)

#plt.title('Band values as a function of time')
ax1.set_ylabel('Band Values')

# plot NDVI
ax2 = fig.add_subplot(212, sharex=ax1)
ax2.plot(time, ndvi, 'o', color="black", label="NDVI")
ax2.grid(True)
start, end = ax2.get_xlim()
ax2.xaxis.set_ticks(np.arange(start, end, 64.5))

# Format the ticks.
years    = mdates.YearLocator()   # every year
months   = mdates.MonthLocator()  # every month
yearsFmt = mdates.DateFormatter('%Y')

ax2.set_ylabel('NDVI')

ax2.xaxis.set_major_locator(years)
ax2.xaxis.set_major_formatter(yearsFmt)
ax2.xaxis.set_minor_locator(months)



In [ ]:

    
# Convert the timestamp to a numpy array
t = np.array([i.toordinal() for i in time])
t

$$ \begin{bmatrix} t_1 & 1 \\\ t_2 & 1 \\\ \vdots & \vdots \\\ t_n & 1 \end{bmatrix} \begin{bmatrix} x_0 \\\ x_1 \end{bmatrix}= \begin{bmatrix} y_1 \\\ y_2 \\\ \vdots \\\ y_n \end{bmatrix} $$$$ \mathbf{A}\mathbf{x} = \mathbf{b} $$$$ \mathbf{x} = \mathbf{A} \backslash \mathbf{b} $$



In [ ]:

    
A = array([ t, ones(len(t))]).transpose()
b = ndvi
x = linalg.lstsq(A,b)[0] # obtaining the parameters
x



In [ ]:

    
b_hat = A.dot(x)



In [ ]:

    
fig2 = figure(figsize=(12,4), dpi=80)
plot(time,b_hat.transpose(),'r-',t,b,'o')
fig2.fmt_xdata = mdates.DateFormatter('%Y-%m-%d')
fig2.autofmt_xdate()



In [ ]: