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

# Remove some bands
data_clean = remove_bands(data)

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

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

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