4729-Plot up our Endmembers.sagews

# Import Stuff
from numpy import array, append, dot
from scipy.optimize import nnls
from scipy.linalg import pinv2


bands = array(range(5))+2

endmembers = array([[  15.16388709,   66.81481205,   85.84328969],
       [   8.11539187,   75.71889195,  117.31059236],
       [  78.72301204,   68.19003156,  120.38751132],
       [  41.42404339,   94.27057946,  124.4524361 ],
       [  12.39576191,   33.45513443,   82.59387353]])

nnls?

File: /home/pete/bin/sage-5.0/local/lib/python2.7/site-packages/scipy/optimize/nnls.py

Type: <type ‘function’>

Definition: nnls(A, b)

Docstring:

Solve argmin_x || Ax - b ||_2 for x>=0.

A : ndarray

Matrix A as shown above.

b : ndarray

Right-hand side vector.

x : ndarray

Solution vector.

rnorm : float

The residual, || Ax-b ||_2.

This is a wrapper for NNLS.F.

Plot up our Endmembers

list_plot(zip(bands,endmembers[:,0]),plotjoined=True,color='green',frame=True,gridlines='automatic',axes_labels=['Landsat Band','DN'])+  list_plot(zip(bands,endmembers[:,1]),plotjoined=True,color='yellow')+list_plot(zip(bands,endmembers[:,2]),plotjoined=True,color='brown')

Now to Unmix

Make a sumulated pixel refelctance

fractions=array([0.1,0.3,0.6])

pixelSpectra=dot(endmembers, fractions)

list_plot(zip(bands,pixelSpectra),plotjoined=True,color='blue',frame=True,gridlines='automatic',axes_labels=['Landsat Band','DN'])

# Conventional unmixing
print dot(pinv2(endmembers),pixelSpectra)

# NNLS unmixing
print nnls(endmembers,pixelSpectra)

[-0.1 0.3 0.9] (array([ 0. , 0.32751853, 0.85135712]), 5.2579330023415265)

# Make a pixel reflectance that's wrong (< 0%)

fractions=array([-0.1,0.3,0.6])

pixelSpectra=dot(endmembers, fractions)

list_plot(zip(bands,pixelSpectra),plotjoined=True,color='blue',frame=True,gridlines='automatic',axes_labels=['Landsat Band','DN'])

# Conventional unmixing
print dot(pinv2(endmembers),pixelSpectra)

# NNLS unmixing
print nnls(endmembers,pixelSpectra)[0]

[-0.1 0.3 0.6] [ 0. 0.32751853 0.55135712]

# One that's wrong (> 100%)

fractions=array([0.3,0.3,0.6])

pixelSpectra=dot(endmembers, fractions)

list_plot(zip(bands,pixelSpectra),plotjoined=True,color='blue',frame=True,gridlines='automatic',axes_labels=['Landsat Band','DN'])

# Conventional unmixing
print dot(pinv2(endmembers),pixelSpectra)

# NNLS unmixing
print nnls(endmembers,pixelSpectra)[0]

[ 0.3 0.3 0.6] [ 0.3 0.3 0.6]

Weighting Example

# Solve the system using the weights

sum2oneWeight=array([35])

endmembersWeighted=append(endmembers,array([sum2oneWeight.repeat(endmembers[0].size)]),axis=0)

print endmembersWeighted

weightedPixelSpectra = append(pixelSpectra,sum2oneWeight)

print

# Conventional unmixing
print dot(pinv2(endmembersWeighted),weightedPixelSpectra)

# NNLS unmixing
print nnls(endmembersWeighted,weightedPixelSpectra)[0]

[[ 15.16388709 66.81481205 85.84328969] [ 8.11539187 75.71889195 117.31059236] [ 78.72301204 68.19003156 120.38751132] [ 41.42404339 94.27057946 124.4524361 ] [ 12.39576191 33.45513443 82.59387353] [ 35. 35. 35. ]] [ 0.25303764 0.21470612 0.66662582] [ 0.25303764 0.21470612 0.66662582]

# One that's wrong wrong

fractions=array([-0.1,0.3,0.9])

pixelSpectra=dot(endmembers, fractions)

list_plot(zip(bands,pixelSpectra),plotjoined=True,color='blue',frame=True,gridlines='automatic',axes_labels=['Landsat Band','DN'])

# Solve the system using the weights

sum2oneWeight=array([35.0])

endmembersWeighted=append(endmembers,array([sum2oneWeight.repeat(endmembers[0].size)]),axis=0)

print endmembersWeighted

weightedPixelSpectra = append(pixelSpectra,sum2oneWeight)

print

# Conventional unmixing
print dot(pinv2(endmembersWeighted),weightedPixelSpectra)

# NNLS unmixing
print nnls(endmembersWeighted,weightedPixelSpectra)[0]

[[ 15.16388709 66.81481205 85.84328969] [ 8.11539187 75.71889195 117.31059236] [ 78.72301204 68.19003156 120.38751132] [ 41.42404339 94.27057946 124.4524361 ] [ 12.39576191 33.45513443 82.59387353] [ 35. 35. 35. ]] [-0.12348118 0.25735306 0.93331291] [ 0. 0.24807895 0.89941565]