%matplotlib inline
We start by initializing the Earth Engine API.
# Load code to setup authorization for an IPython Notebook. Note that this is a temporary step
# that is required until the Earth Engine Python API is updated to include this logic.
%run 'authorize_earth_engine_in_notebook.ipynb'
# Initialize Earth Engine
# Note that we are calling a function defined in the previously run IPython Notebook, rather than
# the typical call to ee.Initialize()
ee_initialize()
Then we import in other Python packages we need.
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.
xBand = 'time'
yBandList = [
'10',
u'20',
u'30',
u'40',
u'50',
u'61',
u'62',
u'70',
u'80',
]
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.
collection = ee.ImageCollection('L7_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.
# 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.
# 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.
iBands = [header.index(b) for b in yBandList]
yData = data[0:,iBands].astype(np.float)
Calculate NDVI
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.
# 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)
# Convert the timestamp to a numpy array
t = np.array([i.toordinal() for i in time])
t
A = array([ t, ones(len(t))]).transpose()
b = ndvi
x = linalg.lstsq(A,b)[0] # obtaining the parameters
x
b_hat = A.dot(x)
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()