notebook.community

Edit and run 





In [1]:



    


import os
import time
import numpy as np
import pysptools.util as util
import pysptools.eea as eea #endmembers extraction algorithms
import pysptools.abundance_maps as amap
import pysptools.classification as cls
import pysptools.material_count as cnt



    











In [2]:



    


load_start_time = time.time()
data_path = os.getcwd()
header_file = 'NEON_D02_SERC_DP3_368000_4306000_reflectance.hdr'
data_file = os.path.join(data_path,header_file)
print('SERC Data Header:',data_file)
data,header = util.load_ENVI_file(data_file)
print('Data took ', time.time() - load_start_time, ' seconds to load.')



    


















    






SERC Data Header: C:\Users\bhass\Desktop\RSDI\RSDI_2018\NEON_D02_SERC_DP3_368000_4306000_reflectance.hdr
Data took  27.226665496826172  seconds to load.

















In [3]:



    


for key in sorted(header.keys()):
  print(key)



    


















    






bands
byte order
coordinate system string
data ignore value
data type
dataset names
description
file type
fwhm
header offset
interleave
lines
map info
reflectance scale factor
samples
wavelength
wavelength units

















In [4]:



    


def remove_bands(M):
    """Remove bands with atmospheric scattering.
    Water vapor bands: 
        Band_Window_1_Nanometers = 1340,1445
        Band_Window_2_Nanometers = 1790,1955
    [191:212] -> 1340.04 -1445.21 nm
    [281:314] -> 1790.76 - 1956.02 nm 
    [416:425] -> 2466.82 - 2511.89 nm"""
    
    #Valid bands: range(0,191), range(212, 281), range(315,415)
    
    vb1 = list(range(0,190)); 
    vb2 = list(range(213,281))
    vb3 = list(range(315,415))
    Mp = M[:,:,vb1+vb2+vb3]
    return Mp

#Display wavelengths corresponding to bands removed:
print('Wavelengths Removed:')
print(header['wavelength'][191],'-',header['wavelength'][212],'nm')
print(header['wavelength'][281],'-',header['wavelength'][314],'nm')
print(header['wavelength'][416],'-',header['wavelength'][425],'nm')

#Display wavelengths corresponding to bands kept:
print('\nWavelengths Kept:')
print(header['wavelength'][0],'-',header['wavelength'][190],'nm')
print(header['wavelength'][213],'-',header['wavelength'][280],'nm')
print(header['wavelength'][315],'-',header['wavelength'][415],'nm')



    


















    






Wavelengths Removed:
1340.044434 - 1445.210449 nm
1790.755981 - 1956.016846 nm
2466.823486 - 2511.894531 nm

Wavelengths Kept:
383.534302 - 1335.036499 nm
1450.218384 - 1790.755981 nm
1961.024902 - 2461.815430 nm

















In [5]:



    


print('Raw Data Dimensions:',data.shape)
#remove bad bands
data_clean = remove_bands(data)
print('Cleaned Data Dimensions:',data_clean.shape)



    


















    






Raw Data Dimensions: (1000, 1000, 426)
Cleaned Data Dimensions: (1000, 1000, 358)

















In [ ]:



    


#adjust number of bands
wavelength_center_bands = np.array(list(map(float, header['wavelength'])))
print('total # of bands:',len(wavelength_center_bands))
wavelength_center_bands_valid = wavelength_center_bands[np.r_[0:190,213:281,315:415]]
print('total # of valid bands:',len(wavelength_center_bands_valid))



    











In [ ]:



    


#remove bad bands in the header wavelengths
valid_band_range = [i for j in (range(0,190), range(213, 281), range(315,415)) for i in j]
header['wavelength'] = [header['wavelength'][i] for i in valid_band_range]
#len(header['wavelength'])



    











In [ ]:



    


print(header['map info'])



    











In [ ]:



    


#Extract Spatial Extent (UTM)
xmin = float(header['map info'][3]); 
ymax = float(header['map info'][4]); 
xmax = xmin + 1000;
ymin = ymax - 1000;

print('xmin:',xmin)
print('xmax:',xmax)
print('ymin:',ymin)
print('ymax:',ymax)



    











In [ ]:



    


#Display linearly stretched RGB image of the SERC tile
#type(data_clean)
#util.display_linear_stretch(data_clean,58,34,19,suffix='SERC')
#wavelength_center_bands_valid
#display(data_clean[:,:,34])

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





In [ ]:



    


#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'

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'





In [ ]:



    


#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.')



    











In [ ]:



    


#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.')



    











In [ ]:



    


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



    











In [ ]:



    


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



    











In [ ]:



    


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



    











In [ ]:



    


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



    











In [ ]:



    


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



    











In [ ]:



    


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



    











In [ ]:



    


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



    











In [ ]:



    


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



    











In [ ]:



    


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



    











In [ ]:



    


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

STOP HERE





In [ ]:



    


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



    











In [ ]:



    


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



    











In [ ]:



    


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



    











In [ ]:



    


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



    











In [ ]:



    


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

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,5,1); plt.title('EM1') amap1_hist = plt.hist(np.ndarray.flatten(amaps[:,:,0]),bins=50,range=[0,0.25])

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

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

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

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

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

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

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

ax9 = fig.add_subplot(2,5,9); plt.title('EM9') amap1_hist = plt.hist(np.ndarray.flatten(amaps[:,:,8]),bins=50,range=[0,0.01])

ax10 = fig.add_subplot(2,5,10); plt.title('EM10') amap1_hist = plt.hist(np.ndarray.flatten(amaps[:,:,9]),bins=50,range=[0,0.001])

Content source: mjones01/NEON-Data-Skills

Similar notebooks:

notebook.community | gallery | about