#ee.display?
wavelength_float = [float(i) for i in header['wavelength']]
ee_axes = {}
ee_axes['wavelength'] = wavelength_float
ee_axes['x']='Wavelength, nm'
ee_axes['y']='Reflectance'

"""
Explanation: Spectral Unmixing allows pixels to be composed of fractions or abundances of each class. 
Endmembers can be thought of as the basis spectra of an image. Once these endmember spectra are determined, the image cube can be 'unmixed' into the fractional abundance of each material in each pixel (Winter, 1999).
End of explanation
"""


#Endmember Extraction (Unmixing) - NFINDR Algorithm (Winter, 1999)
ee_start_time = time.time()
ee = eea.NFINDR()
U = ee.extract(data_clean,4,maxit=5,normalize=False,ATGP_init=True)
ee.display(axes=ee_axes,suffix='SERC')
print('Endmember extraction took ', time.time() - ee_start_time, ' seconds to run.')

#Abundance Maps
amap_start_time = time.time()
am = amap.FCLS()
amaps = am.map(data_clean,U,normalize=False)
am.display(colorMap='jet',columns=4,suffix='SERC')
print('Abundance maps took ', time.time() - amap_start_time, ' seconds to generate.')

#Print mean values of each abundance map to better estimate thresholds
print('Abundance Map Mean Values:')
print('EM1:',np.mean(amaps[:,:,0]))
print('EM2:',np.mean(amaps[:,:,1]))
print('EM3:',np.mean(amaps[:,:,2]))
print('EM4:',np.mean(amaps[:,:,3]))
#print('EM5:',np.mean(amaps[:,:,4]))
#print('EM6:',np.mean(amaps[:,:,5]))
#print('EM7:',np.mean(amaps[:,:,6]))
#print('EM8:',np.mean(amaps[:,:,7]))
#print('EM9:',np.mean(amaps[:,:,8]))
#print('EM10:',np.mean(amaps[:,:,9]))

#Look at histogram of each abundance map to determine ballpark for thresholds to use
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(18,8))

ax1 = fig.add_subplot(2,4,1); plt.title('EM1')
amap1_hist = plt.hist(np.ndarray.flatten(amaps[:,:,0]),bins=50,range=[0,1.0]) 

ax2 = fig.add_subplot(2,4,2); plt.title('EM2')
amap1_hist = plt.hist(np.ndarray.flatten(amaps[:,:,1]),bins=50,range=[0,0.001]) 

ax3 = fig.add_subplot(2,4,3); plt.title('EM3')
amap1_hist = plt.hist(np.ndarray.flatten(amaps[:,:,2]),bins=50,range=[0,0.05]) 

ax4 = fig.add_subplot(2,4,4); plt.title('EM4')
amap1_hist = plt.hist(np.ndarray.flatten(amaps[:,:,3]),bins=50,range=[0,1.0]) 

ax5 = fig.add_subplot(2,4,5); plt.title('EM5')
amap1_hist = plt.hist(np.ndarray.flatten(amaps[:,:,4]),bins=50,range=[0,0.01]) 

ax6 = fig.add_subplot(2,4,6); plt.title('EM6')
amap1_hist = plt.hist(np.ndarray.flatten(amaps[:,:,5]),bins=50,range=[0,0.05]) 

ax7 = fig.add_subplot(2,4,7); plt.title('EM7')
amap1_hist = plt.hist(np.ndarray.flatten(amaps[:,:,6]),bins=50,range=[0,0.00001]) 

ax8 = fig.add_subplot(2,4,8); plt.title('EM8')
amap1_hist = plt.hist(np.ndarray.flatten(amaps[:,:,7]),bins=50,range=[0,0.0001]) 

#Create functions to compute and display SAM and SID 

#Spectral Angle Mapper
def SAM(data,E,thrs=None):
    sam = cls.SAM()
    cmap = sam.classify(data,E,threshold=thrs)
    sam.display(colorMap='Paired',suffix='SERC')

#Spectral Information Divergence
def SID(data,E,thrs=None):
    sid = cls.SID()
    cmap = sid.classify(data,E,threshold=thrs)
    sid.display(colorMap='tab20b',suffix='SERC')

U2 = U[[0,1,2],:]

SID(data_clean, U2, [0.8,0.3,0.4])

U3 = U[[0,2,3],:]
SID(data_clean, U3, [0.4,0.5,0.03])

def display(image,plot_title=''):
    import matplotlib.pyplot as plt
    img = plt.imshow(image,extent=[xmin,xmax,ymin,ymax])
    img.set_cmap('gray')
    ax=plt.gca()
    plt.setp(ax.get_xticklabels(),rotation=90); #rotate x tick labels 90 degrees
    plt.title(plot_title)
    plt.show()
    plt.clf()

#man_made = (amaps[:,:,1]>0.10)+(amaps[:,:,2]>0.10)+(amaps[:,:,3]>0.6)
#display(man_made)
#display((amaps[:,:,5]>0.1))
display((amaps[:,:,2]>0.035))
#display((amaps[:,:,2]>0.03)+(amaps[:,:,5]>0.05),'SERC Roads and Buildings')

#Display USA Topo Base Map
#Reference: https://viewer.nationalmap.gov/
from IPython.core.display import Image
UStopo_filename = 'SERC_368000_4307000_UStopo.PNG'
UStopo_image = os.path.join(data_path,UStopo_filename)
print('Base Map')
Image(filename=UStopo_image,width=225,height=225)

"""
Explanation: create dictionary from header['wavelength'] list in order to plot using ee.display
ee_axes = {}
ee_axes['wavelength'] = {k: v for v,k in enumerate(header['wavelength'])}
ee_axes['x']='Wavelength, nm'
ee_axes['y']='Reflectance'
End of explanation
"""


# U2 is a library made of 4 selected endmembers
SAM(data_clean, U2, [0.6,0.02,0.02,0.5])

# Try another combination of 6 endmembers (adding 9)
U2 = U[[1,2,3,4,5,8],:]
SAM(data_clean, U2, [0.5,0.5,0.1,0.5,0.5,0.5])
SID(data_clean, U2, [0.5,0.5,0.05,0.5,0.3,0.5])

U2 = U[[3,4,5,8],:]
SID(data_clean, U2, [0.1,0.5,0.5,0.5])

# Try with just 3 endmembers
U2 = U[[3,4,5,8],:]
SAM(data_clean, U2, [0.5,0.1,0.5,0.5])

##Package Installation (Can also install required packages from command line - activate py35 env first)
#https://stackoverflow.com/questions/45156080/installing-modules-to-anaconda-from-tar-gz
#import sys
#!{sys.executable} -m pip install spectral
#!{sys.executable} -m pip install "C:\Users\bhass\Downloads\pysptools-0.14.2.tar.gz
#!conda install --yes --prefix {sys.prefix} numpy
#!conda install --yes --prefix {sys.prefix} matplotlib
#!conda install --yes --prefix {sys.prefix} scikit-learn
#!conda install --yes --prefix {sys.prefix} cvxopt

"""
Explanation: STOP HERE
End of explanation
"""


