import ee
ee.Initialize()

"""
Explanation: We start by initializing the Earth Engine API.
End of explanation
"""


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

"""
Explanation: Then we import in other Python packages we need.
End of explanation
"""


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)

"""
Explanation: Next we define the Landsat bands that we would like to plot, along with the starting and ending times.
End of explanation
"""


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()

"""
Explanation: Next we contruct a filtered ImageCollection for the date range, and extract band information for a specified point location.
End of explanation
"""


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

"""
Explanation: We separate the information returned into column headers and data.
End of explanation
"""


# 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))]

"""
Explanation: Next we extract time information and convert it to at Python datatime data type.
End of explanation
"""


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

"""
Explanation: Extract the data columns what we want to display on the plot.
End of explanation
"""


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

"""
Explanation: Calculate NDVI
End of explanation
"""


# 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

"""
Explanation: And finally, we create a plot of MODIS band values as a function of time.
End of explanation
"""


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()

"""
Explanation: $$
\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} $$
End of explanation
"""


