'''
Adjustments to the pysptools library to support linear spectral mixture
analysis (LSMA). Includes classes and functions:
* `PPI`
* `NFINDR`
* `FCLSAbundanceMapper`
* `combine_endmembers_and_normalize()`
* `convex_hull_graham()`
* `endmembers_by_maximum_angle()`
* `endmembers_by_maximum_area()`
* `endmembers_by_maximum_volume()`
* `endmembers_by_query()`
* `hall_rectification()`
* `iterate_endmember_combinations()`
* `normalize_reflectance_within_image()`
* `predict_spectra_from_abundance()`
* `point_to_pixel_geometry()`
* `ravel()`
* `ravel_and_filter()`
* `report_raster_dynamic_range()`
* `subtract_endmember_and_normalize()`
'''
import itertools
import json
import os
import re
import numpy as np
from concurrent.futures import ProcessPoolExecutor
from functools import reduce, partial, wraps
from unmixing.transform import mnf_rotation
from unmixing.utils import array_to_raster, as_array, as_raster, dump_raster, partition, xy_to_pixel, pixel_to_xy, spectra_at_xy, rmse
from lxml import etree
from osgeo import gdal, ogr, osr
from pykml.factory import KML_ElementMaker as KML
import pysptools.eea as sp_extract
import pysptools.abundance_maps as sp_abundance
import pysptools.classification as sp_classify
import pysptools.material_count as sp_matcount
class AbstractAbundanceMapper(object):
def __init__(self, mixed_raster, gt, wkt, nodata=-9999, processes=1):
assert np.all(np.greater(mixed_raster.shape, 0)), 'Raster array cannot have any zero-length axis'
self.shp = mixed_raster.shape
self.mixed_raster = mixed_raster
self.hsi = mixed_raster.T # HSI form: (p x m x n) as (n x m x p)
self.gt = gt
self.wkt = wkt
self.nodata = nodata
self.num_processes = processes
class AbstractExtractor(object):
def get_idx_as_kml(self, path, gt, wkt, data_dict=None):
'''
Exports a KML file containing the locations of the extracted endmembers
as point markers.
'''
# Despite that the HSI cube is the transpose of our raster array, the
# coordinates returned by `get_idx()` are already in the right order
# (longitude, latitude) because the (m by n) == (y by x) order
# transposed is (n by m) == (x by y); the row index is the latitude
# and the column index is the longitude.
coords = pixel_to_xy(self.get_idx(), gt=gt, wkt=wkt, dd=True)
if any(map(lambda x: x[0] == 0 and x[1] == 0, self.get_idx())):
print('Warning: Target endmember chosen at (0,0)')
print('One or more endmembers may be photometric shade')
if data_dict is None:
data_dict = {
'wavelength': range(1, len(coords) + 1),
'wavelength units': 'MNF Component',
'z plot titles': ['', '']
}
ico = 'http://maps.google.com/mapfiles/kml/paddle/%i.png'
pmarks = []
for i, pair in enumerate(coords):
pmarks.append(KML.Placemark(
KML.Style(
KML.IconStyle(
KML.Icon(KML.href(ico % (i + 1))))),
KML.name(data_dict['wavelength units'] + ' %d' % (i + 1)),
KML.Point(KML.coordinates('%f,%f' % pair))))
doc = KML.kml(KML.Folder(*pmarks))
with open(path, 'wb') as source:
source.write(etree.tostring(doc, pretty_print=True))
def get_idx_as_shp(self, path, gt, wkt):
'''
Exports a Shapefile containing the locations of the extracted
endmembers. Assumes the coordinates are in decimal degrees.
'''
coords = pixel_to_xy(self.get_idx(), gt=gt, wkt=wkt, dd=True)
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.CreateDataSource(path)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
layer = ds.CreateLayer(path.split('.')[0], srs, ogr.wkbPoint)
for pair in coords:
feature = ogr.Feature(layer.GetLayerDefn())
# Create the point from the Well Known Text
point = ogr.CreateGeometryFromWkt('POINT(%f %f)' % pair)
feature.SetGeometry(point) # Set the feature geometry
layer.CreateFeature(feature) # Create the feature in the layer
feature.Destroy() # Destroy the feature to free resources
# Destroy the data source to free resources
ds.Destroy()
class PPI(sp_extract.PPI, AbstractExtractor):
pass
class NFINDR(sp_extract.NFINDR, AbstractExtractor):
pass
class FCLSAbundanceMapper(AbstractAbundanceMapper):
'''
A class for generating an abundance map, containing both the raw spectral
(mixed) data and the logic to unmix the data into an abundance map with
the fully constrained least-squares (FCLS) approach. The "full"
constraints are the sum-to-one and non-negativity constraints. Given
q endmembers and p spectral bands, the mapper is forced to find
abundances within a simplex in a (q-1)-dimensional subspace.
NOTE: The mixed_raster provided in instantation is assumed to correspond
exactly to the mixing space used to induce endmembers; if mixed_raster
was MNF- or PCA-transformed, for instance, the endmember spectra provided
in, e.g., map_abundances(), must match this transformation exactly.
Arguments:
mixed_raster The raster to be unmixed; should NOT be in HSI form
but should be MNF- or PCA-transformed to match any
endmember spectra provided
gt A GDAL GeoTransform tuple for the mixed_raster
wkt Projection information as Well-Known Text for the
mixed_raster
nodata The NoData value for mixed_raster
processes The number of processes to create in mixture analysis
'''
def __init__(self, *args, **kwargs):
super(FCLSAbundanceMapper, self).__init__(*args, **kwargs)
self.mapper = sp_abundance.FCLS()
def __lsma__(self, cases, endmembers):
# For regular LSMA with single endmember spectra
# c is number of pixels, k is number of bands
c, k = cases.shape if len(cases.shape) > 1 else (1, cases.shape[0])
return self.mapper.map(cases.reshape((1, c, k)), endmembers,
normalize = False)
def __mesma__(self, array_pairs):
# For multiple endmember spectra, in chunks
cases, endmembers = array_pairs
# c is number of pixels, k is number of bands
c, k = cases.shape if len(cases.shape) > 1 else (1, cases.shape[0])
return [
self.mapper.map(cases[i,...].reshape((1, 1, k)), endmembers[i,...],
normalize = False) for i in range(0, c)
]
def __mesma2__(self, array_pairs):
# For multiple endmember spectra, pixel-wise
# NOTE: This pixel-wise implementation might be slower than __mesma__
# for large arrays
cases, endmembers = array_pairs
# c is number of pixels, k is number of bands
c, k = cases.shape if len(cases.shape) > 1 else (1, cases.shape[0])
return self.mapper.map(cases.reshape((1, c, k)), endmembers,
normalize = False)
def map_abundance(self, endmembers, pixelwise=False):
'''
Arguments:
endmembers A numpy.ndarray of endmembers; either (q x p) array
of q endmembers and p bands (for regular LSMA) or a
(c x q x p) array, where c = m*n, for multiple
endmember spectra for each pixel.
Returns: An (m x n x q) numpy.ndarray (in HSI form) that contains
the abundances for each of q endmember types.
'''
q = endmembers.shape[-2]
# FCLS with the sum-to-one constraint has an extra degree of freedom so it
# is able to form a simplex of q corners in (q-1) dimensions:
# q <= n (Settle and Drake, 1993)
k = q - 1 # Find q corners of simplex in (q-1) dimensions
endmembers = endmembers[...,0:k]
shp = self.hsi.shape
base_array = self.hsi[:,:,0:k].reshape((shp[0] * shp[1], k))
# Get indices for each process' work range
work = partition(base_array, self.num_processes, axis=0)
with ProcessPoolExecutor(max_workers = self.num_processes) as executor:
# We're working with multiple endmembers
if endmembers.ndim == 3 and pixelwise:
result = executor.map(self.__mesma2__, [ # Work done pixel-wise
(base_array[i,...], endmembers[i,...]) for i in range(0, base_array.shape[0])
])
elif endmembers.ndim == 3:
result = executor.map(self.__mesma__, [
(base_array[i:j,...], endmembers[i:j,...]) for i, j in work
])
# We're working with a single endmember per class
else:
# Curry an unmixing function with the present endmember array
unmix = partial(self.__lsma__, endmembers = endmembers)
result = executor.map(unmix, [
base_array[i:j,...] for i, j in work
])
combined_result = list(result) # Executes the multiprocess suite
if endmembers.ndim == 3 and not pixelwise:
# When chunking with multiple endmembers, we get list of lists
ext_array = [y for x in combined_result for y in x] # Flatten once
return np.concatenate(ext_array, axis = 1)\
.reshape((shp[0], shp[1], q))
return np.concatenate(combined_result, axis = 1)\
.reshape((shp[0], shp[1], q))
def validate_by_forward_model(
self, ref_image, abundances, ref_spectra=None,
ref_em_locations=None, dd=False, nodata=-9999, r=10000,
as_pct=True, convert_nodata=False):
'''
Validates LSMA result in the forward model of reflectance, i.e.,
compares the observed reflectance in the original (mixed) image to the
abundance predicted by a forward model of reflectance using the
provided endmember spectra. NOTE: Does not apply in the case of
multiple endmember spectra; requires only one spectral profile per
endmember type.
Arguments:
ref_image A raster array of the reference spectra (not MNF-
transformed data).
abundances A raster array of abundances; a (q x m x n) array for
q abundance types (q endmembers).
ref_spectra With single endmember spectra, user can provide the
reference spectra, e.g., the observed reflectance for
each endmember (not MNF spectra).
ref_em_locations With single endmember spectra, user can provide
the coordinates of each endmember, so that reference
spectra can be extracted for validation.
dd True if ref_em_locations provided and the coordinates
are in decimal degrees.
nodata The NoData value to use.
r The number of random samples to take in calculating
RMSE.
as_pct Report normalized RMSE (as a percentage).
convert_nodata
True to convert all NoData values to zero (as in zero
reflectance)
'''
rastr = ref_image.copy()
assert (ref_spectra is not None) or (ref_em_locations is not None), 'When single endmember spectra are used, either ref_spectra or ref_em_locations must be provided'
if ref_spectra is not None:
assert ref_spectra.shape[0] == abundances.shape[0], 'One reference spectra must be provided for each endmember type in abundance map'
else:
# Get the spectra for each endmember from the reference dataset
ref_spectra = spectra_at_xy(ref_image, ref_em_locations,
self.gt, self.wkt, dd = dd)
# Convert the NoData values to zero reflectance
if convert_nodata:
rastr[rastr == nodata] = 0
ref_spectra[ref_spectra == nodata] = 0
shp = rastr.shape # Reshape the arrays
arr = rastr.reshape((shp[0], shp[1]*shp[2]))
# Generate random sampling indices
idx = np.random.choice(np.arange(0, arr.shape[1]), r)
# Get the predicted reflectances
preds = predict_spectra_from_abundance(ravel(abundances), ref_spectra)
assert preds.shape == arr.shape, 'Prediction and observation matrices are not the same size'
# Take the mean RMSE (sum of RMSE divided by number of pixels), after
# the residuals are normalized by the number of endmembers
rmse_value = rmse(arr, preds, idx, n = ref_spectra.shape[0], nodata = nodata).sum() / r
norm = 1
if as_pct:
# Divide by the range of the measured data; minimum is zero
norm = arr.max()
return str(round(rmse_value / norm * 100, 2)) + '%'
return round(rmse_value / norm, 2)
def combine_endmembers_and_normalize(
abundances, es=(1, 2), at_end=True, nodata=-9999):
'''
Combines two endmembers from a fraction image into a single endmember.
If the original endmember abundances summed to one, they will sum to one
in the resulting image as well. Arguments:
abundances The raster array of endmember abundances
es A two-element tuple of the endmembers indices to combine
at_end Place the combined endmembers at the end of the array?
nodata The NoData value to ignore
'''
shp = abundances.shape
rast = abundances.copy() # Copy raster array
rast[rast == nodata] = 0 # Replace NoData values
c0 = rast[es[0], ...] # Get the endmembers to be combined
c1 = rast[es[1], ...]
# Stack the remaining bands
abunds = []
for e in range(0, shp[0]):
if e not in es:
abunds.append(rast[e, ...])
if at_end:
comps = (abunds, c0 + c1.reshape(1, shp[1], shp[2]))
else:
comps = (c0 + c1.reshape(1, shp[1], shp[2]), abunds)
rast = None
return np.vstack(comps)
def convex_hull_graham(points, indices=False):
'''
Returns points on convex hull of an array of points in CCW order according
to Graham's scan algorithm. By Tom Switzer <thomas.switzer@gmail.com>.
Arguments:
points The points for which a convex hull is sought
indices True to return a tuple of (indices, hull)
'''
TURN_LEFT, TURN_RIGHT, TURN_NONE = (1, -1, 0)
def cmp(a, b):
return (a > b) - (a < b)
def turn(p, q, r):
return cmp((q[0] - p[0])*(r[1] - p[1]) - (r[0] - p[0])*(q[1] - p[1]), 0)
def keep_left(hull, r):
while len(hull) > 1 and turn(hull[-2], hull[-1], r) != TURN_LEFT:
hull.pop()
if not len(hull) or hull[-1] != r:
hull.append(r)
return hull
pts_sorted = sorted(points)
l = reduce(keep_left, pts_sorted, [])
u = reduce(keep_left, reversed(pts_sorted), [])
hull = l.extend(u[i] for i in range(1, len(u) - 1)) or l
if indices:
return ([points.index(h) for h in hull], hull)
return hull
def endmembers_by_maximum_angle(
rast, targets, ref_target, gt=None, wkt=None, dd=False):
'''
Locates endmembers in (2-dimensional) feature space as the triad (3-corner
simplex) that maximizes the angle formed with a reference endmember target.
Returns the endmember coordinates in feature (not geographic) space.
Arguments:
rast The raster that describes the feature space
ref_target The coordinates (in feature space) of a point held fixed
targets The coordinates (in feature space) of all other points
gt The GDAL GeoTransform
wkt The GDAL WKT projection
dd True for coordinates in decimal degrees
Angle calculation from:
http://stackoverflow.com/questions/2827393/angles-between-two-n-dimensional-vectors-in-python/13849249#13849249
'''
def unit_vector(vector):
# Returns the unit vector of the vector.
return vector / np.linalg.norm(vector)
def angle_between(v1, v2):
# Returns the angle in radians between vectors 'v1' and 'v2'
v1_u = unit_vector(v1)
v2_u = unit_vector(v2)
return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))
# Can accept either a gdal.Dataset or numpy.array instance
if not isinstance(rast, np.ndarray):
rastr = rast.ReadAsArray()
gt = rast.GetGeoTransform()
wkt = rast.GetProjection()
else:
assert gt is not None and wkt is not None, 'gt and wkt arguments required'
rastr = rast.copy()
# Get the spectra for these targets; this works in two dimensions only
ref_spec = spectra_at_xy(rast, (ref_target,), gt, wkt)[...,0:2].reshape((2,))
target_specs = spectra_at_xy(rast, targets, gt, wkt)[...,0:2]
# All combinations of 2 of the targets
combos = list(itertools.combinations(range(max(target_specs.shape)), 2))
spec_map = [
[target_specs[i,:] for i in triad] for triad in combos
]
coord_map = [
[targets[i] for i in triad] for triad in combos
]
# Find vectors from ref_spec, not from origin (by vector subtraction)
# If (cx) is the ref_spec vector (line from origin to ref_spec),
# and (ca) and (cb) are the vectors to the points that form the angle
# (axb), then [(cx) - (ca)] and [(cx) - (cb)] are the vectors from point
# x to the points a and b, respectively.
vectors = [(ref_spec - a, ref_spec - b) for a, b in spec_map]
angles = [angle_between(v1, v2) for v1, v2 in vectors]
idx = angles.index(max(angles))
specs = spec_map[idx] # The optimized spectra
locs = coord_map[idx] # The optimized coordinates
specs.insert(0, ref_spec) # Add the reference target
locs.insert(0, ref_target) # Add the reference coordinates
return (np.array(specs), locs)
def endmembers_by_maximum_area(
rast, targets, area_dim=2, gt=None, wkt=None, dd=False):
'''
Find up to four (4) endmembers by findng the maximum volume of the mixing
space defined by every possible combination of endmembers where each
endmember type appears only once. Arguments:
rast The raster whose mixing space is to be investigated
targets The pseudo-invariant feature sets, a dictionary
area_dim The number of dimensions to use in area calculation;
only 2 or 3 dimensions are supported
gt The GDAL GeoTransform
wkt The GDAL WKT projection
dd True for coordinates in decimal degrees
'''
def area(a, b, c) :
# Courtesy of Corentin Lapeyre
# (http://code.activestate.com/recipes/576896-3-point-area-finder/)
return 0.5 * np.linalg.norm(np.cross(b-a, c-a))
spec_map, coord_map = iterate_endmember_combinations(rast, targets,
ref_target=None, ndim=3, gt=gt, wkt=wkt, dd=dd)
area_totals = [area(*[f[0:area_dim] for f in each]) for each in spec_map]
idx = area_totals.index(max(area_totals)) # Which set is optimal?
specs = spec_map[idx] # The optimized spectra
locs = list(coord_map[idx])
return (np.array(specs), locs)
def endmembers_by_maximum_volume(
rast, targets, ref_target=None, ndim=3, gt=None, wkt=None, dd=False):
'''
Arguments:
rast The raster that describes the feature space
targets The coordinates (in feature space) of all other points
ref_target (Optional) Can constrain the volume formed by one point
ndim The number of dimensions in which to calculate the solution
gt The GDAL GeoTransform
wkt The GDAL WKT projection
dd True for coordinates in decimal degrees
Angle calculation from:
http://stackoverflow.com/questions/2827393/angles-between-two-n-dimensional-vectors-in-python/13849249#13849249
NOTE: Argument ndim determines how many endmembers are returned; the
production of combinations depends on a square matrix, so ndim must also
be equal to or less than the dimensionality of the image array.
'''
def calc_volume(spec_map):
# Calculate the absolute volume of the determinant of the volume spanned
# by the construction vectors; this is proportional to the true volume
# of the mixing space spanned by these endmembers
return list(map(np.abs, map(np.linalg.det, spec_map)))
spec_map, coord_map = iterate_endmember_combinations(
rast, targets, ref_target, ndim, gt, wkt, dd)
volumes = calc_volume(spec_map)
idx = volumes.index(max(volumes))
specs = spec_map[idx] # The optimized spectra
locs = list(coord_map[idx]) # The optimized coordinates
if ref_target is not None:
locs = list(locs)
locs.insert(0, ref_target)
return (np.array(specs), locs)
def endmembers_by_query(rast, query, gt, wkt, dd=False):
'''
Returns a list of endmember locations based on a provided query, e.g.:
> query = rast[1,...] < -25 # Band 2 should be less than -25
> endmembers_by_query(rast, query, gt, wkt)
Arguments:
rast The raster array to find endmembers within
query A NumPy boolean array representing a query in the feature space
gt The GDAL GeoTransform
wkt The GDAL WKT projection
dd True for coordinates in decimal degrees
'''
assert isinstance(rast, np.ndarray), 'Requires a NumPy array'
shp = rast.shape
idx = np.indices((shp[-2], shp[-1]))
# Execute query on the indices (pixel locations), then return the coordinates
return list(pixel_to_xy([
(x, y) for y, x in idx[:,query].T
], gt, wkt, dd=dd))
def hall_rectification(
reference, subject, out_path, ref_set, sub_set, dd=False,
nodata=-9999, dtype=np.int32, keys=('High/Bright', 'Low/Dark'),
verbose=False):
'''
Performs radiometric rectification after Hall et al. (1991) in Remote
Sensing of Environment. Assumes first raster is the reference image and
that none of the targets are NoData pixels in the reference image (they
are filtered out in the subject images). Arguments:
reference The reference image, a gdal.Dataset
subject The subject image, a gdal.Dataset
out_path Path to a directory where the rectified images should be stored
ref_set Sequence of two sequences: "bright" radiometric control set,
then "dark" radiometric control set for reference image
sub_set As with ref_set, a sequence of sequences (e.g., list of two
lists): [[<bright targets>], [<dark targets]]
dd Coordinates are in decimal degrees?
dtype Date type (NumPy dtype) for the array; default is 32-bit Int
nodata The NoData value to use fo all the rasters
keys The names of the dictionary keys for the bright, dark sets,
respectively
'''
# Unpack bright, dark control sets for subject image
bright_targets, dark_targets = (sub_set[keys[0]], sub_set[keys[1]])
# Calculate the mean reflectance in each band for bright, dark targets
brights = spectra_at_xy(reference, ref_set[keys[0]], dd=dd)
bright_ref = brights[brights[:,0] != nodata].mean(axis = 0)
darks = spectra_at_xy(reference, ref_set[keys[1]], dd=dd)
dark_ref = darks[darks[:,0] != nodata].mean(axis = 0)
if verbose:
print('-- Bright reference:', str(np.round(bright_ref, 0).tolist()))
print('-- Dark reference:', str(np.round(dark_ref, 0).tolist()))
# Calculate transformation for the target image
brights = spectra_at_xy(subject, bright_targets, dd=dd) # Prepare to filter NoData pixels
darks = spectra_at_xy(subject, dark_targets, dd=dd)
# Get the "subject" image means for each radiometric control set
mean_bright = brights[brights[:,0] != nodata].mean(axis = 0)
mean_dark = darks[darks[:,0] != nodata].mean(axis = 0)
if verbose:
print('-- Bright subject mean:', str(np.round(mean_bright, 0).tolist()))
print('-- Dark subject mean:', str(np.round(mean_dark, 0).tolist()))
# Calculate the coefficients of the linear transformation
m = (bright_ref - dark_ref) / (mean_bright - mean_dark)
b = (dark_ref * mean_bright - mean_dark * bright_ref) / (mean_bright - mean_dark)
arr = subject.ReadAsArray()
shp = arr.shape # Remember the original shape
mask = arr.copy() # Save the NoData value locations
m = m.reshape((shp[0], 1))
b = b.reshape((shp[0], 1)).T.repeat(shp[1] * shp[2], axis=0).T
arr2 = ((arr.reshape((shp[0], shp[1] * shp[2])) * m) + b).reshape(shp)
arr2[mask == nodata] = nodata # Re-apply NoData values
# Dump the raster to a file
out_path = os.path.join(out_path, 'rect_%s' % os.path.basename(subject.GetDescription()))
dump_raster(
array_to_raster(
arr2, subject.GetGeoTransform(), subject.GetProjection(),
dtype=dtype), out_path)
def iterate_endmember_combinations(
rast, targets, ref_target=None, ndim=3, gt=None, wkt=None, dd=False):
'''
Creates all possible combinations of endmembers from a common pool or from
among groups of possible endmembers (when `targets` is a dictionary). When
a dictionary is provided, endmember combinations will contain (only) one
endmember from each group, where the groups are defined by dictionary keys,
i.e., `{'group1': [(x, y), ...], 'group2': [(x, y), ...], ...}`.
Arguments:
rast Input raster array
targets Possible endmembers; as a list or dictionary
ref_target (Optional) Constrain the optimization to always include
this endmember
ndim Number of dimensions to limit the search to
gt (Optional) A GDAL GeoTransform, required for array input
wkt (Optional) A GDAL WKT projection, required for array input
dd True if the target coordinates are in decimal degrees
'''
# Can accept either a gdal.Dataset or numpy.array instance
if not isinstance(rast, np.ndarray):
rastr = rast.ReadAsArray()
gt = rast.GetGeoTransform()
wkt = rast.GetProjection()
else:
assert gt is not None and wkt is not None, 'gt and wkt arguments required'
rastr = rast.copy()
# `targets` is a dictionary
if isinstance(targets, dict):
# Get the spectra for these targets; this works in two dimensions only
target_specs = {}
for label, cases in targets.items():
target_specs[label] = spectra_at_xy(rast, targets[label], gt, wkt)[...,0:ndim]
ncom = ndim # Determinant only defined for square matrices
if ref_target is not None:
assert ndim == len(targets.keys()) + 1, 'Number of groups among target endmembers should be one less than the dimensionality when ref_target is used'
ncom -= 1 # If reference target used, form combinations with 1 fewer
ref_spec = spectra_at_xy(rast, (ref_target,), gt, wkt,
dd=dd)[...,0:ndim].reshape((ndim,))
# Find all possible combinations of (ncom) of these spectra
spec_map = list(itertools.product(*[target_specs[label] for label in target_specs.keys()]))
coord_map = list(itertools.product(*[t[1] for t in targets.items()]))
else:
# Get the spectra for these targets; this works in two dimensions only
target_specs = spectra_at_xy(rast, targets, gt, wkt)[...,0:ndim]
ncom = ndim # Determinant only defined for square matrices
if ref_target is not None:
ncom -= 1 # If reference target used, form combinations with 1 fewer
ref_spec = spectra_at_xy(rast, (ref_target,), gt, wkt,
dd=dd)[...,0:ndim].reshape((ndim,))
# Find all possible combinations of (ncom) of these spectra
combos = list(itertools.combinations(range(max(target_specs.shape)), ncom))
spec_map = [[target_specs[i,:] for i in triad] for triad in combos]
coord_map = [[targets[i] for i in triad] for triad in combos]
# Add the reference target to each combination
if ref_target is not None:
spec_map = list(map(list, spec_map))
for spec in spec_map:
# FIXME Cannot use insert with tuples when dictionary input is provided
spec.insert(0, ref_spec)
return (spec_map, coord_map)
def normalize_reflectance_within_image(
rast, band_range=(0, 5), nodata=-9999, scale=100):
'''
Following Wu (2004, Remote Sensing of Environment), normalizes the
reflectances in each pixel by the average reflectance *across bands.*
This is an attempt to mitigate within-endmember variability. Arguments:
rast A gdal.Dataset or numpy.array instance
nodata The NoData value to use (and value to ignore)
scale (Optional) Wu's definition scales the normalized reflectance
by 100 for some reason; another reasonable value would
be 10,000 (approximating scale of Landsat reflectance units);
set to None for no scaling.
'''
# Can accept either a gdal.Dataset or numpy.array instance
if not isinstance(rast, np.ndarray):
rastr = rast.ReadAsArray()
else:
rastr = rast.copy()
shp = rastr.shape
b0, b1 = band_range # Get the beginning, end of band range
b1 += 1 # Ranges in Python are not inclusive, so add 1
rastr_normalized = np.divide(
rastr.reshape((shp[0], shp[1]*shp[2])),
rastr[b0:b1,...].mean(axis=0).reshape((1, shp[1]*shp[2])).repeat(shp[0], axis=0))
# Recover original shape; scale if necessary
rastr_normalized = rastr_normalized.reshape(shp)
if scale is not None:
rastr_normalized = np.multiply(rastr_normalized, scale)
# Fill in the NoData areas from the original raster
np.place(rastr_normalized, rastr == nodata, nodata)
return rastr_normalized
def point_to_pixel_geometry(
points, source_epsg=None, target_epsg=None, pixel_side_length=30):
'''
Where points is a list of X,Y tuples and X and Y are coordinates in
meters, returns a series of OGR Polygons where each Polygon is the
pixel extent with a given point at its center. Assumes square pixels.
Arguments:
points Sequence of X,Y numeric pairs or OGR Point geometries
source_epsg The EPSG code of the source projection (Optional)
target_epsg The EPSG code of the target projection (Optional)
pixel_side_length The length of one side of the intended pixel
'''
polys = []
# Convert points to atomic X,Y pairs if necessary
if isinstance(points[0], ogr.Geometry):
points = [(p.GetX(), p.GetY()) for p in points]
source_ref = target_ref = None
if all((source_epsg, target_epsg)):
source_ref = osr.SpatialReference()
target_ref = osr.SpatialReference()
source_ref.ImportFromEPSG(source_epsg)
target_ref.ImportFromEPSG(target_epsg)
transform = osr.CoordinateTransformation(source_ref, target_ref)
for p in points:
r = pixel_side_length / 2 # Half the pixel width
ring = ogr.Geometry(ogr.wkbLinearRing) # Create a ring
vertices = [
(p[0] - r, p[1] + r), # Top-left
(p[0] + r, p[1] + r), # Top-right, clockwise from here...
(p[0] + r, p[1] - r),
(p[0] - r, p[1] - r),
(p[0] - r, p[1] + r) # Add top-left again to close ring
]
# Coordinate transformation
if all((source_ref, target_ref)):
vertices = [transform.TransformPoint(*v)[0:2] for v in vertices]
for vert in vertices:
ring.AddPoint(*vert)
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
polys.append(poly)
return polys
def predict_spectra_from_abundance(abundances, endmembers):
'''
Predicts the mixed image from the endmember spectra and the fractional
abundances, i.e. R = AS + e, where R is the reflectance of a given pixel,
S is the matrix of endmember spectra, and A is the matrix of abundances.
Endmembers should be in (p x q) form, where each column corresponds to the
spectra of a given endmember. Return either a (p x 1) vector of the
predicted (mixed) spectra or a (p x (m*n)) matrix of the predicted (mixed)
spectra for (m*n) total abundance estimates. Arguments:
abundances (c x q) vector of abundances for q abundance types for
c total pixels (1 or more)
endmembers (q x p) vector of endmember spectra for p spectral bands
'''
return np.dot(abundances, endmembers).swapaxes(0, 1)
def ravel(arr):
'''
Reshapes a (p, m, n) array to ((m*n), p) where p is the number of
dimensions. Assumes the first axis is the shortest. Arguments:
arr A NumPy array with shape (p, m, n)
'''
return ravel_and_filter(arr, filter=False)
def ravel_and_filter(arr, filter=True, nodata=-9999):
'''
Reshapes a (p, m, n) array to ((m*n), p) where p is the number of
dimensions and, optionally, filters out the NoData values. Assumes the
first axis is the shortest. Arguments:
arr A NumPy array with shape (p, m, n)
filter True to filter out NoData values (otherwise, only ravels)
nodata The NoData value; only used in filtering
'''
shp = arr.shape
# If the array has already been raveled
if len(shp) == 1 and filter:
return arr[arr != nodata]
# If a "single-band" image
if len(shp) == 2:
arr = arr.reshape(1, shp[-2]*shp[-1]).swapaxes(0, 1)
if filter:
return arr[arr != nodata]
# For multi-band images
else:
arr = arr.reshape(shp[0], shp[1]*shp[2]).swapaxes(0, 1)
if filter:
return arr[arr[:,0] != nodata]
return arr
def report_raster_dynamic_range(
path, bands=(1,2,3,4,5,7), tpl='HDF4_EOS:EOS_GRID:"%s":Grid:sr_band%d', lj=40):
'''
Prints out the dynamic range of a given raster, averaged across the bands.
Arguments:
path The file path to a raster
bands The bands within which dynamic range should be calculated
tpl The template for the GDAL subdataset names
lj The left justification (for text output)
'''
# ComputeStatistics() returns the (0) minimum, (1) maximum, (2) mean and
# (3) standard deviation; DR is max - min
dr = lambda stats: stats[1] - stats[0]
stddev = lambda stats: stats[3]
dr_by_band = [] # Dynamic range by band
stddev_by_band = []
rast, gt0, wkt = as_raster(path)
# If the dataset is an HDF or some other linked dataset
if rast.GetRasterBand(1) is None:
for b in bands:
subrast = gdal.Open(tpl % (path, b))
try:
stats = subrast.GetRasterBand(1).ComputeStatistics(False)
except AttributeError:
pass
dr_by_band.append(dr(stats))
stddev_by_band.append(dr(stats))
# If, instead, the dataset is a flat stack like a GeoTIFF
else:
for i, b in enumerate(bands):
stats = rast.GetRasterBand(i + 1).ComputeStatistics(False)
dr_by_band.append(dr(stats))
stddev_by_band.append(dr(stats))
print('{:.2f} ({:.0f} max std. dev.) -- {:s}'.format(
np.mean(dr_by_band), max(stddev_by_band),
os.path.basename(path).ljust(lj)))
def subtract_endmember_and_normalize(abundances, e):
'''
Subtract the endmember at index e from the raster cube and normalize
the remaining endmembers so that they sum to one. Arguments:
abundances A raster array cube of abundance estimates
e The index of the endmember to subtract
'''
f = e + 1 # Index after endmember index
shp = abundances.shape
stack = np.vstack((abundances[0:e, ...], abundances[f:shp[0], ...]))\
.reshape((shp[0] - 1, shp[1]*shp[2]))
return np.apply_along_axis(lambda x: x / x.sum(), 0, stack)\
.reshape((shp[0] - 1, shp[1], shp[2]))
