Pine Creek Example 2

Man made like buildings and drives can be characterized by the endmembers EM4 and EM8. We simply add the two related abundance maps and show the result below.

In [1]:
# Run on Python 2.7 and 3.x

from __future__ import print_function

%matplotlib inline

import os
import os.path as osp
import sys
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

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

def display(image):
    import matplotlib.pyplot as plt
    img = plt.imshow(image)
    img.set_cmap('gray')
    plt.show()
    plt.clf()

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)
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, 75, 34, 0)

<matplotlib.figure.Figure at 0x7f817b60dad0>

Extract twelve endmembers

In [2]:
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 0x7f817b5aeb50>

Invert

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

<matplotlib.figure.Figure at 0x7f8175b91f90>
In [4]:
# apply a threshold and add the two abundance maps
#man_made = (amaps[:,:,3] > 0.02)+(amaps[:,:,7] > 0.1)
man_made = (amaps[:,:,2] > 0.05)+(amaps[:,:,8] > 0.05)+(amaps[:,:,9] > 0.05)
print('Roads and buildings')
display(man_made)

Roads and buildings

<matplotlib.figure.Figure at 0x7f8175ca4990>

We can compare with a map of the region. The result is not bad but there is false positive.

In [5]:
# Display a USGS Quadrangle map
# 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_quadrangle_map.png'), width=250, height=250)
Out[5]: