Pine Creek Example 1

This example follows the Hematite and Methanol examples. With one major difference, Hematite and Methanol cubes are generated in laboratory conditions and are easier to analyse. For the Pine Creek example, the HSI cube comes from a satellite.

Before the analysis, four bands ranges were removed. They correspond to atmospheric scattering. Twelve endmembers was extracted and their corresponding abundance maps generated by NNLS. After inspecting the abundances maps, five endmembers are kept and are used to create the classes’ maps. For the classification step SAM and SID is used. For this particular case, SID gives a better result than SAM.

Finally, at the end of this document you can compare the endmembers classification to a supervised ECHO clustering made by David Landgrebe with the MultiSpec software. ECHO clustering do a right classification for 16 features, ECHO is a clear winner here.

In conclusion, endmembers extraction carries useful information and can be use to do a 'preview or quick' classification. However, it does not beat a standard supervised clustering algorithm. Abundance maps can be analysed individually. By example, EM4 and EM8 show many of the buildings and roads and EM2 give a good distribution of the forest. The former is the subject for the next example on Pine Creek.

In [1]:
# Run on Python 2.7 and 3.3

from __future__ import print_function

%matplotlib inline

import os
import os.path as osp
import numpy as np

import pysptools.util as util
import pysptools.eea as eea
import pysptools.abundance_maps as amp
import pysptools.classification as cls
import pysptools.material_count as cnt

def remove_bands(M):
    """
    Remove the bands with atmospheric
    scattering.
    Remove:
        [0..4]
        [102..110]
        [148..169]
        [211..end]
    """
    p1 = list(range(5,102))
    p2 = list(range(111,148))
    p3 = list(range(170,211))
    Mp = M[:,:,p1+p2+p3]
    return Mp

data_path = os.environ['PYSPTOOLS_DATA']

sample = '92AV3C.hdr'

data_file = osp.join(data_path, sample)
data, header = util.load_ENVI_file(data_file)

# Remove some bands
data_clean = remove_bands(data)
# ajust the number of bands
header['wavelength'] = range(1, data_clean.shape[2]+1)

# Display a linear stretched RGB image of the Pine Creek cube
#util.display_linear_stretch(data, 98, 86, 22)
util.display_linear_stretch(data, 169, 83, 45)

<matplotlib.figure.Figure at 0x7fe2ac716f50>

HySime do a good guess.

In [2]:
hy = cnt.HySime()
kf, Ek = hy.count(data_clean)
print('Testing HySime')
print('  Virtual dimensionality is: k =', kf)

Testing HySime
  Virtual dimensionality is: k = 14
In [3]:
def SAM(data, E, thrs=None):
    sam = cls.SAM()
    cmap = sam.classify(data, E, threshold=thrs)
    sam.display(colorMap='Paired', suffix='Pine Creek')

def SID(data, E, thrs=None):
    sid = cls.SID()
    cmap = sid.classify(data, E, threshold=thrs)
    sid.display(colorMap='Paired', suffix='Pine Creek')

We extract 12 endmembers.

In [4]:
ee = eea.NFINDR()
U = ee.extract(data_clean, 12, maxit=5, normalize=False, ATGP_init=True)
ee.display(header, suffix='Pine Creek')

<matplotlib.figure.Figure at 0x7fe2ac690fd0>

And generate the corresponding abundance maps. Abundance maps 2, 3, 6, 7 and 10 are used for the classification step.

In [5]:
am = amp.FCLS()
amaps = am.map(data_clean, U, normalize=False)
am.display(colorMap='jet', columns=3, suffix='Pine Creek')

<matplotlib.figure.Figure at 0x7fe2a2855e50>
In [6]:
# U2 is a library made of 5 selected endmembers
U2 = U[[1,2,5,6,9],:]
SAM(data_clean, U2, [0.1,0.3,0.2,0.2,0.15])
SID(data_clean, U2, [0.1,0.3,0.05,0.1,0.03])
In [10]:
# Display an ECHO classification
# Reference:
#   Landgrebe, David, 1997, Multispectral Data Analysis: A Signal Theory Perspective, 
#   School of Electrical Engineering, Purdue University.
from IPython.core.display import Image
home = os.environ['HOME']
Image(filename=osp.join(home, 'dev/pysptools/examples/pine_creek_echo_class.png'), width=400, height=400)
Out[10]: