[GRASS-SVN] r72212 - in grass-addons/grass7/imagery: . i.pysptools.unmix
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Feb 6 14:31:51 PST 2018
Author: sbl
Date: 2018-02-06 14:31:51 -0800 (Tue, 06 Feb 2018)
New Revision: 72212
Added:
grass-addons/grass7/imagery/i.pysptools.unmix/
grass-addons/grass7/imagery/i.pysptools.unmix/Makefile
grass-addons/grass7/imagery/i.pysptools.unmix/i.pysptools.unmix.html
grass-addons/grass7/imagery/i.pysptools.unmix/i.pysptools.unmix.py
Log:
New addon i.pysptools.unmix
Added: grass-addons/grass7/imagery/i.pysptools.unmix/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.pysptools.unmix/Makefile (rev 0)
+++ grass-addons/grass7/imagery/i.pysptools.unmix/Makefile 2018-02-06 22:31:51 UTC (rev 72212)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.pysptools.unmix
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/imagery/i.pysptools.unmix/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/imagery/i.pysptools.unmix/i.pysptools.unmix.html
===================================================================
--- grass-addons/grass7/imagery/i.pysptools.unmix/i.pysptools.unmix.html (rev 0)
+++ grass-addons/grass7/imagery/i.pysptools.unmix/i.pysptools.unmix.html 2018-02-06 22:31:51 UTC (rev 72212)
@@ -0,0 +1,104 @@
+<em>i.pysptools.unmix</em> extracts endmembers from imagery group and performs
+spectral unmixing using
+<a href="https://pysptools.sourceforge.io/">pysptools</a>.
+
+<p>
+The module is a wrapper around the <b>pysptools</b> Python library, that integrates
+its functionality for
+<a href="https://pysptools.sourceforge.io/eea.html">Endmember Extraction</a>
+and <a href="https://pysptools.sourceforge.io/abundance_maps.html">Spectral
+Unmixing</a> into GRASS GIS.
+</p>
+<p>
+It requires that the Python libraries <b>pysptools</b> and <b>scikit-learn</b>
+are installed.
+</p>
+<p>
+Supported algorithms for
+<a href="https://pysptools.sourceforge.io/eea.html">Endmember Extraction</a>
+are:
+<ul>
+<li><i>NFINDR</i>: N-FINDR endmembers induction algorithm after Winter (1999),
+that also makes use of an Automatic Target Generation Process (ATGP) (Plaza &
+Chang 2006). (<i>Default</i>)</li>
+<li><i>FIPPI</i>: Fast Iterative Pixel Purity Index after Chang (2006)</li>
+<li><i>PPI</i>: Pixel Purity Index</li>
+</ul>
+</p>
+<p>
+Supported algorithms for
+<a href="https://pysptools.sourceforge.io/abundance_maps.html">Spectral Unmixing</a>
+ are:
+<ul>
+<li><i>FCLS</i>: Fully Constrained Least Squares (FCLS): Estimates endmember
+abundance per pixel with the constraint that values are non-negative and sum up
+to one per pixel (<i>Default</i>)</li>
+<li><i>UCLS</i>: Unconstrained Least Squares (UCLS): Estimates endmember
+abundance per pixel in an unconstrained way</li>
+<li><i>NNLS</i>: Non-negative Constrained Least Squares (NNLS): Estimates endmember
+abundance per pixel with the constraint that values are non-negative</li>
+</ul>
+</p>
+
+<h2>NOTES</h2>
+
+<p>
+Number of endmembers to extract (<i>endmember_n</i>) is supposed to be lower
+than the number of bands in the imagery group. Only the <i>PPI</i> method can
+extract more endmembers than there are bands in the imagery group.
+</p>
+
+<h2>EXAMPLES</h2>
+
+<div class="code"><pre>
+# List bands
+bands=`g.list type=raster pattern=lsat7_2002* separator=','`
+
+# Create imagery group
+i.group group=lsat_2002 input="$bands"
+
+# Extract Endmembers and perform spectral unmixing using pysptools
+i.pysptools.unmix input=lsat_2002 endmembers=endmembers endmember_n=5 \
+output=spectrum.txt prefix=lsat_spectra --v
+
+# Compare to result from i.spec.unmix
+i.spec.unmix group=lsat7_2002 matrix=sample/spectrum.dat result=lsat7_2002_unmix \
+error=lsat7_2002_unmix_err iter=lsat7_2002_unmix_iterations
+</pre></div>
+
+
+<h2>REQUIREMENTS</h2>
+
+<ul>
+ <li><a href="https://pypi.python.org/pypi/pysptools">pysptools library</a></li>
+ <li><a href="https://pypi.python.org/pypi/scikit-learn">scikit-learn library</a></li>
+</ul>
+
+<h2>REFERENCES</h2>
+<p>
+Chang, C.-I. 2006: A fast iterative algorithm for implementation of pixel
+purity index. Geoscience and Remote Sensing Letters, IEEE, 3(1): 63-67.
+</p>
+<p>
+Plaza, A. & Chang, C.-I. 2006: Impact of Initialization on Design of
+Endmember Extraction Algorithms. Geoscience and Remote Sensing,
+IEEE Transactions. 44(11): 3397-3407.
+</p>
+<p>
+Winter, M. E. 1999: N-FINDR: an algorithm for fast autonomous spectral
+end-member determination in hyperspectral data. Presented at the
+Imaging Spectrometry V, Denver, CO, USA, (3753): 266-275.
+</p>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="i.spec.unmix.html">i.spec.unmix</a>
+</em>
+
+<h2>AUTHORS</h2>
+
+Stefan Blumentrath, <a href="http://www.nina.no">
+Norwegian Institute for Nature Research (NINA), Oslo, Norway</a><br>
+Zofie Cimburova, <a href="http://www.nina.no">
+Norwegian Institute for Nature Research (NINA), Oslo, Norway</a>
Property changes on: grass-addons/grass7/imagery/i.pysptools.unmix/i.pysptools.unmix.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/imagery/i.pysptools.unmix/i.pysptools.unmix.py
===================================================================
--- grass-addons/grass7/imagery/i.pysptools.unmix/i.pysptools.unmix.py (rev 0)
+++ grass-addons/grass7/imagery/i.pysptools.unmix/i.pysptools.unmix.py 2018-02-06 22:31:51 UTC (rev 72212)
@@ -0,0 +1,346 @@
+#!/usr/bin/env python
+
+"""
+MODULE: i.pysptools.unmix
+
+AUTHOR(S): Stefan Blumentrath < stefan.blumentrath AT nina.no>
+ Zofie Cimburova < zofie.cimburova AT nina.no>
+
+PURPOSE: Extract endmembers from imagery group and perform spectral unmixing
+ using pysptools
+ Depemds on pysptools and scikit-learn
+
+COPYRIGHT: (C) 2018 by the GRASS GIS Development Team,
+ Norwegian Institute for Nature Research (NINA)
+
+ This program is free software under the GNU General Public
+ License (>=v2). Read the file COPYING that comes with GRASS
+ for details.
+"""
+
+#%module
+#% description: Extract endmembers from imagery group and perform spectral unmixing using pysptools
+#% keyword: imagery
+#% keyword: endmember
+#% keyword: spectral unmixing
+#%end
+
+#%option G_OPT_I_GROUP
+#% key: input
+#% description: Input imagery group
+#% required : yes
+#%end
+
+#%option G_OPT_F_OUTPUT
+#% key: output
+#% guisection: output
+#% description: Text file storing endmember information for i.spec.unmix
+#% required : no
+#%end
+
+#%option
+#% key: prefix
+#% description: Prefix for resulting raster maps
+#% guisection: output
+#% required : no
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: endmembers
+#% description: Vector map representing identified endmembers
+#% guisection: output
+#% required : no
+#%end
+
+#%option
+#% key: endmember_n
+#% type: integer
+#% description: Number of endmembers to identify
+#% required: yes
+#%end
+
+#%option
+#% key: extraction_method
+#% type: string
+#% description: Method for endmember extraction
+#% options: FIPPI,PPI,NFINDR
+#% answer: NFINDR
+#%end
+
+#%option
+#% key: unmixing_method
+#% type: string
+#% description: Algorithm for spectral unmixing
+#% options: FCLS,UCLS,NNLS
+#% answer: FCLS
+#%end
+
+#%option
+#% key: maxit
+#% type: integer
+#% description: Maximal number of iterations for endmember extraction (default=3*number of bands)
+#% required: no
+#%end
+
+#%flag
+#% key: n
+#% description: Do not use Automatic Target Generation Process (ATGP)
+#%end
+
+#%rules
+#% required: output,prefix
+#%end
+
+import os
+import sys
+import numpy as np
+from grass.pygrass import raster as r
+from grass.pygrass.utils import getenv
+import grass.script as gs
+
+if not "GISBASE" in os.environ.keys():
+ gs.message("You must be in GRASS GIS to run this program.")
+ sys.exit(1)
+
+
+def get_rastertype(raster):
+
+ if type(raster[0,0]) != np.float32 and type(raster[0,0]) != np.float64:
+ map_type = u'INTEGER'
+ else:
+ map_type = u'REAL'
+
+ return map_type
+
+def mask_rasternd(raster):
+ if type(raster[0,0]) != np.float32 and type(raster[0,0]) != np.float64:
+ mask = raster != -2147483648
+ else:
+ mask = np.isnan(raster) == False
+
+ return mask
+
+
+def main():
+
+ try:
+ import pysptools.eea as eea
+ except ImportError:
+ gs.fatal(_("Cannot import pysptools \
+ (https://pypi.python.org/pypi/pysptools) library."
+ " Please install it (pip install pysptools)"
+ " or ensure that it is on path"
+ " (use PYTHONPATH variable)."))
+
+ try:
+ # sklearn is a dependency of used pysptools functionality
+ import sklearn
+ except ImportError:
+ gs.fatal(_("Cannot import sklearn \
+ (https://pypi.python.org/pypi/scikit-learn) library."
+ " Please install it (pip install scikit-learn)"
+ " or ensure that it is on path"
+ " (use PYTHONPATH variable)."))
+
+ # Parse input options
+ input = options['input']
+ output = options['output']
+ prefix = options['prefix']
+ endmember_n = int(options['endmember_n'])
+ endmembers = options['endmembers']
+ maxit = int(options['maxit'])
+ extraction_method = options['extraction_method']
+ unmixing_method = options['unmixing_method']
+ atgp_init = True if not flags['n'] else False
+
+ # List maps in imagery group
+ try:
+ maps = gs.read_command('i.group', flags='g', group=input,
+ quiet=True).rstrip('\n').split('\n')
+ except:
+ pass
+
+ # Validate input
+ # q and maxit can be None according to manual, but does not work in curtrent pysptools version
+ if endmember_n <= 0:
+ gs.fatal('Number of endmembers has to be > 0')
+ """if (extraction_method == 'PPI' or
+ extraction_method == 'NFINDR'):
+ gs.fatal('Extraction methods PPI and NFINDR require endmember_n >= 2')
+ endmember_n = None"""
+
+ if maxit <= 0:
+ maxit = 3 * len(maps)
+
+ if endmember_n > len(maps) + 1:
+ gs.warning('More endmembers ({}) requested than bands in \
+ input imagery group ({})'.format(endmember_n, len(maps)))
+ if extraction_method != 'PPI':
+ gs.fatal('Only PPI method can extract more endmembers than number \
+ of bands in the imagery group')
+
+ if not atgp_init and extraction_method != 'NFINDR':
+ gs.verbose('ATGP is only taken into account in \
+ NFINDR extraction method...')
+
+ # Get metainformation from input bands
+ band_types = {}
+ img = None
+ n = 0
+ gs.verbose('Reading imagery group...')
+ for m in maps:
+ map = m.split('@')
+
+ # Build numpy stack from imagery group
+ raster = r.raster2numpy(map[0], mapset=map[1])
+ if raster == np.float64:
+ raster = float32(raster)
+ gs.warning('{} is of type Float64.\
+ Float64 is currently not supported.\
+ Reducing precision to Float32'.format(raster))
+
+ # Determine map type
+ band_types[map[0]] = get_rastertype(raster)
+
+ # Create cube and mask from GRASS internal NoData value
+ if n == 0:
+ img = raster
+ # Create mask from GRASS internal NoData value
+ mask = mask_rasternd(raster)
+ else:
+ img = np.dstack((img, raster))
+ mask = np.logical_and((mask_rasternd(raster)), mask)
+
+ n = n + 1
+
+ # Read a mask if present and give waringing if not
+ # Note that otherwise NoData is read as values
+ gs.verbose('Checking for MASK...')
+ try:
+ MASK = r.raster2numpy('MASK', mapset=getenv('MAPSET')) == 1
+ mask = np.logical_and(MASK, mask)
+ MASK = None
+ except:
+ pass
+
+ if extraction_method == 'NFINDR':
+ # Extract endmembers from valid pixels using NFINDR function from pysptools
+ gs.verbose('Extracting endmembers using NFINDR...')
+ nfindr = eea.NFINDR()
+ E = nfindr.extract(img, endmember_n, maxit=maxit, normalize=False,
+ ATGP_init=atgp_init, mask=mask)
+ elif extraction_method == 'PPI':
+ # Extract endmembers from valid pixels using PPI function from pysptools
+ gs.verbose('Extracting endmembers using PPI...')
+ ppi = eea.PPI()
+ E = ppi.extract(img, endmember_n, numSkewers=10000, normalize=False,
+ mask=mask)
+ elif extraction_method == 'FIPPI':
+ # Extract endmembers from valid pixels using FIPPI function from pysptools
+ gs.verbose('Extracting endmembers using FIPPI...')
+ fippi = eea.FIPPI()
+ # q and maxit can be None according to manual, but does not work
+ """if not maxit and not endmember_n:
+ E = fippi.extract(img, q=None, normalize=False, mask=mask)
+ if not maxit:
+ E = fippi.extract(img, q=endmember_n, normalize=False, mask=mask)
+ if not endmember_n:
+ E = fippi.extract(img, q=int(), maxit=maxit, normalize=False,
+ mask=mask)
+ else:
+ E = fippi.extract(img, q=endmember_n, maxit=maxit, normalize=False,
+ mask=mask)"""
+ E = fippi.extract(img, q=endmember_n, maxit=maxit, normalize=False,
+ mask=mask)
+
+ # Write output file in format required for i.spec.unmix addon
+ if output:
+ gs.verbose('Writing spectra file...')
+ n = 0
+ with open(output, 'w') as o:
+ o.write('# Channels: {}\n'.format('\t'.join(band_types.keys())))
+ o.write('# Wrote {} spectra line wise.\n#\n'.format(endmember_n))
+ o.write('Matrix: {0} by {1}\n'.format(endmember_n, len(maps)))
+ for e in E:
+ o.write('row{0}: {1}\n'.format(n, '\t'.join([str(i) for i in e])))
+ n = n + 1
+
+ # Write vector map with endmember information if requested
+ if endmembers:
+ gs.verbose('Writing vector map with endmembers...')
+ from grass.pygrass import utils as u
+ from grass.pygrass.gis.region import Region
+ from grass.pygrass.vector import Vector
+ from grass.pygrass.vector import VectorTopo
+ from grass.pygrass.vector.geometry import Point
+
+ # Build attribute table
+ # Deinfe columns for attribute table
+ cols = [(u'cat', 'INTEGER PRIMARY KEY')]
+ for b in band_types.keys():
+ cols.append((b.replace('.','_'), band_types[b]))
+
+ # Get region information
+ reg = Region()
+
+ # Create vector map
+ new = Vector(endmembers)
+ new.open('w', tab_name=endmembers, tab_cols=cols)
+
+ cat = 1
+ for e in E:
+ # Get indices
+ idx = np.where((img[:,:]==e).all(-1))
+
+ # Numpy array is ordered rows, columns (y,x)
+ if len(idx[0]) == 0 or len(idx[1]) == 0:
+ gs.warning('Could not compute coordinated for endmember {}. \
+ Please consider rescaling your data to integer'.format(cat))
+ cat = cat + 1
+ continue
+
+ coords = u.pixel2coor((idx[1][0], idx[0][0]), reg)
+ point = Point(coords[1] + reg.ewres / 2.0,
+ coords[0] - reg.nsres / 2.0)
+
+ # Get attributes
+ n = 0
+ attr = []
+ for b in band_types.keys():
+ if band_types[b] == u'INTEGER':
+ attr.append(int(e[n]))
+ else:
+ attr.append(float(e[n]))
+ n = n + 1
+
+ # Write geometry with attributes
+ new.write(point, cat=cat,
+ attrs=tuple(attr))
+ cat = cat + 1
+
+ # Close vector map
+ new.table.conn.commit()
+ new.close(build=True)
+
+ if prefix:
+ # Run spectral unmixing
+ import pysptools.abundance_maps as amaps
+ if unmixing_method == 'FCLS':
+ fcls = amaps.FCLS()
+ result = fcls.map(img, E, normalize=False, mask=mask)
+ elif unmixing_method == 'NNLS':
+ nnls = amaps.NNLS()
+ result = nnls.map(img, E, normalize=False, mask=mask)
+ elif unmixing_method == 'UCLS':
+ ucls = amaps.UCLS()
+ result = ucls.map(img, E, normalize=False, mask=mask)
+
+ # Write results
+ for l in range(endmember_n):
+ rastname = '{0}_{1}'.format(prefix, l + 1)
+ r.numpy2raster(result[:,:,l], 'FCELL', rastname)
+
+# Run the module
+if __name__ == "__main__":
+ options, flags = gs.parser()
+ sys.exit(main())
Property changes on: grass-addons/grass7/imagery/i.pysptools.unmix/i.pysptools.unmix.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list