# 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 [2]:

    %matplotlib inline
    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

    def remove_bands(M):
        """
        Remove the bands with atmospheric
        scattering.
        Remove:
            [0..4]
            [102..110]
            [148..169]
            [211..end]
        """
        p1 = range(5,102)
        p2 = range(111,148)
        p3 = 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 = '../data1'
    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 0x16b1a940>

Extract twelve endmembers

In [3]:

    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 0x16b354e0>

Invert

In [4]:

    am = amp.FCLS()
    amaps = am.map(data_clean, U, normalize=False)
    am.display(colorMap='jet', suffix='Pine Creek')

    <matplotlib.figure.Figure at 0xb3ec518>

In [5]:

    # apply a threshold and add the two abundance maps
    #man_made = (amaps[:,:,3] > 0.02)+(amaps[:,:,7] > 0.1)
    man_made = (amaps[:,:,3] > 0.05)+(amaps[:,:,7] > 0.1)
    print 'Roads and buildings'
    display(man_made)

    Roads and buildings

    <matplotlib.figure.Figure at 0x174f4d68>

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

In [6]:

    # 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
    Image(filename=r'C:\dev\pysptools\examples\pine_creek_quadrangle_map.png')

Out[6]: