[GRASS-SVN] r70060 - in grass-addons/grass7/imagery: . i.cva
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Dec 11 13:29:33 PST 2016
Author: ciupava
Date: 2016-12-11 13:29:33 -0800 (Sun, 11 Dec 2016)
New Revision: 70060
Added:
grass-addons/grass7/imagery/i.cva/
grass-addons/grass7/imagery/i.cva/Makefile
grass-addons/grass7/imagery/i.cva/i.cva.html
grass-addons/grass7/imagery/i.cva/i.cva.py
grass-addons/grass7/imagery/i.cva/i.cva_1angle.png
grass-addons/grass7/imagery/i.cva/i.cva_2angleclass.png
grass-addons/grass7/imagery/i.cva/i.cva_3magnitude.png
grass-addons/grass7/imagery/i.cva/i.cva_4change.png
Log:
Committing folder i.cva in grass_addons for grass7
Added: grass-addons/grass7/imagery/i.cva/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.cva/Makefile (rev 0)
+++ grass-addons/grass7/imagery/i.cva/Makefile 2016-12-11 21:29:33 UTC (rev 70060)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR=/home/anna/grass7_trunk
+
+PGM = i.cva
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/imagery/i.cva/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/imagery/i.cva/i.cva.html
===================================================================
--- grass-addons/grass7/imagery/i.cva/i.cva.html (rev 0)
+++ grass-addons/grass7/imagery/i.cva/i.cva.html 2016-12-11 21:29:33 UTC (rev 70060)
@@ -0,0 +1,139 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.cva</em> calculates Change Vector Analysis (CVA) for two input
+variables. CVA is a remote sensing technique used for change detection
+analysis. As input for CVA, two maps for each date must be given: in general,
+on X axis an indicator of overall reflectance and on Y axis an indicator of
+vegetation conditions. A common choice for the indicators is Albedo and NDVI
+(Normalized Difference Vegetation Index) or the Brightness and Greenness
+features features of the Tasselled Cap (TC) transform.
+<br>
+For each pixel of the original image, CVA gives in output a map of the angle
+and a map of the magnitude of the vector of the change between two dates.
+<br>
+Read Malila et al. for a complete explanation of the technique. This module
+might require a first transformation of the data to Top Of Atmosphere
+Reflectance (TOAR); if the TC transform are chosen as indicators, the TC
+transform should be then performed as well before running CVA.
+<p>
+Four parameters are required in input:
+<ul>
+<li> <b>xaraster</b>: first date map for X axis,
+<li> <b>xbraster</b>: second date map for X axis,
+<li> <b>yaraster</b>: first date map for Y axis,
+<li> <b>ybraster</b>: second date map for Y axis.
+</ul>
+<p>
+The following maps can be generated in output:
+<ul>
+<li> <i>basename</i>_angle: map of the angles of the change vector between
+the two dates;
+<li> <i>basename</i>_angle_class: map of the angles, classified by the four
+quadrants (0-90, 90-180, ...);
+<li> <i>basename</i>_magnitude: map of the magnitudes of the change vector
+between the two dates;
+<li> <i>basename</i>_change: final map of the change
+</ul>
+<p>
+The change detection map is created using the classified angle map and applying
+a threshold to the magnitude: the change is given by the pixels that have
+values higher than the threshold, divided in four categories depending on
+the quadrant they belong to.
+<br>
+The threshold can be chosen manually (<em>custom value</em>, given by personal
+criteria) or using statistical criteria. In this case the mean of the magnitude
+values is used and the user can choose the multiples of <em>N</em> standard
+deviation to sum to the mean (threshold = mean + N * standard deviation).
+<br>
+One could consider of running the module at first without assigning a
+threshold, in order to have an idea of the range of the magnitude and to
+choose an appropriate custom threshold (for univariate statystical parameters
+run <a href="r.univar.html">r.univar</a>). In this case <em>i.cva</em> gives
+in output only three maps: the angle, angle classified and magnitude maps.
+<h2>EXAMPLE</h2>
+
+Calculation of CVA maps from North Carolina Landsat 7 ETM scene, between
+lsat5_1987 and lsat7_2002.
+<br>
+The Tasselled cap maps are calculated for TOAR data.
+<div class="code"><pre>
+i.cva xaraster=lsat5_1987_tasscap.1 xbraster=lsat7_2002_tasscap.1 \
+yaraster=lsat5_1987_tasscap.2 ybraster=lsat7_2002_tasscap.2 \
+output=CVA_87_02 stat_threshold=1
+
+Calculating DeltaX and DeltaY
+Writing angle map CVA_87_02_angle
+Writing magnitude map CVA_87_02_magnitude
+Mean of magnitude values is: 0.091335330260002
+Standard deviation of magnitude values is: 0.0671211630131731
+Writing change detection map CVA_87_02_change
+Threshold is 0.158456493273
+</pre></div>
+
+Results:
+
+<p>
+<center>
+ <table border=1>
+ <tr>
+ <td align=center>
+ <img src="i.cva_1angle.png" alt="Angle map">
+ <br>
+ <font size="-1">
+ <i>CVA angle map 1</i>
+ </font>
+ </td>
+ <td align=center>
+ <img src="i.cva_2angleclass.png" alt="Classified angle
+ map">
+ <br>
+ <font size="-1">
+ <i>CVA classified angle map</i>
+ </font>
+ </td>
+ </tr>
+ <tr>
+ <td align=center>
+ <img src="i.cva_3magnitude.png" alt="Magnitude map">
+ <br>
+ <font size="-1">
+ <i>CVA magnitude map</i>
+ </font>
+ </td>
+ <td align=center>
+ <img src="i.cva_4change.png">
+ <br>
+ <font size="-1">
+ <i>CVA change map</i>
+ </font>
+ </td>
+ </tr>
+ </table>
+</center>
+<br>
+
+
+<h2>REFERENCES</h2>
+
+<ul>
+<li>Malila W A, Lafayette W. Change Vector Analysis: An Approach
+for Detecting Forest Changes with Landsat.
+LARS Symp. 1980;326-335.
+</ul>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="i.albedo.html">i.albedo</a>,
+<a href="i.vi.html">i.vi</a>,
+<a href="i.aster.toar.html">i.aster.toar</a>,
+<a href="i.landsat.toar.html">i.landsat.toar</a>,
+<a href="r.univar.html">r.univar</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+Anna Zanchetta
+<p><i>Last changed: $Date: 2016-12-11 23:00:00 +0200 (Sun, 11 Dec
+2016) $</i>
+
Property changes on: grass-addons/grass7/imagery/i.cva/i.cva.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/imagery/i.cva/i.cva.py
===================================================================
--- grass-addons/grass7/imagery/i.cva/i.cva.py (rev 0)
+++ grass-addons/grass7/imagery/i.cva/i.cva.py 2016-12-11 21:29:33 UTC (rev 70060)
@@ -0,0 +1,184 @@
+#!/usr/bin/env python
+
+#%module
+#% description: Performs Change Vector Analysis (CVA) in two dimensions
+#% keyword: imagery
+#%end
+#%option G_OPT_R_INPUT
+#% key: xaraster
+#% description: Name of the first raster for X axis
+#%end
+#%option G_OPT_R_INPUT
+#% key: xbraster
+#% description: Name of the the second raster for X axis
+#%end
+#%option G_OPT_R_INPUT
+#% key: yaraster
+#% description: Name of the first raster for Y axis
+#%end
+#%option G_OPT_R_INPUT
+#% key: ybraster
+#% description: Name of the second raster for Y axis
+#%end
+#%option G_OPT_R_BASENAME_OUTPUT
+#% label: Name for output basename raster maps (angle and magnitude)
+#%end
+#%option
+#% key: custom_threshold
+#% description: Use a custom threshold
+#% guisection: Magnitude threshold
+#% type: double
+#% required: no
+#% descriptions: Insert numerical value for the threshold to perform the analysis
+#% multiple: no
+#%end
+#%option
+#% key: stat_threshold
+#% description: Use a statystical parameter for the threshold (mean + N * standard deviation)
+#% guisection: Magnitude threshold
+#% type: double
+#% required: no
+#% descriptions: Insert the integer value for a multiple of standard deviation (to be summed to the mean of the magnitude values )
+#% multiple: no
+#%end
+#%rules
+#% exclusive: custom_threshold, stat_threshold
+#%end
+
+from __future__ import print_function
+import atexit, sys
+import math
+import grass.script as grass
+from grass import script
+import subprocess
+from subprocess import PIPE
+import grass.temporal as tgis
+
+TMPRAST = []
+
+def delta_calculation(delta, second_band, first_band):
+ """
+ Calculating the Delta as difference between the second band and the first band
+ """
+ equation = ('$delta = $second_band - $first_band')
+ grass.mapcalc(equation, delta=delta, second_band=second_band, first_band=first_band)
+
+def angle_calculation(anglemap, deltaX, deltaY):
+ equation = ('$anglemap = atan($deltaX,$deltaY)')
+ grass.mapcalc(equation, anglemap=anglemap, deltaX=deltaX, deltaY=deltaY)
+
+def magnitude_calculation(magnitudemap, deltaX, deltaY):
+ equation = ('$magnitudemap = sqrt((($deltaX)^2)+(($deltaY)^2))')
+ grass.mapcalc(equation, magnitudemap=magnitudemap, deltaX=deltaX, deltaY=deltaY)
+
+def change_map_calculation(change_map, magnitude_map, threshold, angle_map_class):
+ equation = ('$change_map = if($magnitude_map>$threshold,$angle_map_class,null())')
+ grass.mapcalc(equation, change_map=change_map, magnitude_map=magnitude_map, threshold=threshold, angle_map_class=angle_map_class)
+
+def main():
+ options, flags = grass.parser()
+ xAmap = options['xaraster']
+ xBmap = options['xbraster']
+ yAmap = options['yaraster']
+ yBmap = options['ybraster']
+ output_basename = options['output']
+ custom_threshold = options['custom_threshold']
+ stat_threshold = options['stat_threshold']
+ Xdelta_name = 'deltaX'
+ Ydelta_name = 'deltaY'
+ anglemap_name = output_basename + '_angle'
+ anglemap_class = anglemap_name + '_class'
+ magnitudemap_name = output_basename + '_magnitude'
+ changemap_name = output_basename + '_change'
+
+ if not grass.find_file(name=xAmap, element='cell')['file']:
+ grass.fatal("xaraster map <%s> not found" % xAmap)
+ if not grass.find_file(name=xBmap, element='cell')['file']:
+ grass.fatal("xbraster map <%s> not found" % xBmap)
+ if not grass.find_file(name=yAmap, element='cell')['file']:
+ grass.fatal("yaraster map <%s> not found" % yAmap)
+ if not grass.find_file(name=xBmap, element='cell')['file']:
+ grass.fatal("ybraster map <%s> not found" % yBmap)
+
+ TMPRAST.append(Xdelta_name)
+ TMPRAST.append(Ydelta_name)
+
+ """
+ Calculating Delta for X and Y bands
+ """
+ grass.message(_("Calculating DeltaX and DeltaY"))
+ delta_calculation(Xdelta_name, xBmap, xAmap)
+ delta_calculation(Ydelta_name, yBmap, yAmap)
+
+ """
+ Calculating angle and magnitude maps
+ """
+ grass.message(_("Writing angle map %s") % anglemap_name)
+ angle_calculation(anglemap_name, Xdelta_name, Ydelta_name)
+
+ grass.message(_("Writing magnitude map %s") % magnitudemap_name)
+ magnitude_calculation(magnitudemap_name, Xdelta_name, Ydelta_name)
+
+ """
+ Reclassifing angle map to get a map with the four quadrants
+ """
+ keys = ['1', '2', '3', '4']
+ vals = [0, 90, 180, 270, 360]
+ rvals = [(int(vals[i-1]), int(vals[i]), keys[i-1], vals[i-1], vals[i]) for i in range(1, len(vals))]
+ rules = '\n'.join(['%3d thru %3d = %s %s-%s' % v for v in rvals])
+ script.write_command('r.reclass', input=anglemap_name, output=anglemap_class, rules='-', overwrite=True, stdin=rules.encode())
+
+ """
+ Going to generate the change detection map using the given threshold
+ """
+ if custom_threshold:
+ threshold = custom_threshold
+ grass.message(_("Threshold is %s") % threshold)
+ grass.message(_("Writing change detection map %s") % changemap_name)
+ """
+ Creating final map of the change, using a custom threshold
+ """
+ change_map_calculation(changemap_name, magnitudemap_name, threshold, anglemap_class)
+ elif stat_threshold:
+ """
+ Getting values of mean and standard dev of magnitude to calculate the change detection criteria (> mean + N*stdev)
+ """
+ univar = grass.read_command('r.univar', map=magnitudemap_name, flags='g')
+
+ found = 0
+ for line in univar.splitlines():
+ name,val = line.split('=')
+ if name == 'mean':
+ grass.message(_("Mean of magnitude values is: %s") % val)
+ mean = val
+ found += 1
+ if name == 'stddev':
+ grass.message(_("Standard deviation of magnitude values is: %s") % val)
+ stddev = val
+ found += 1
+ if found != 2:
+ grass.fatal("Couldn\'t find mean or stddev!")
+
+ adding_value = float(stat_threshold) * float(stddev)
+ threshold = float(mean) + float(adding_value)
+ grass.message(_("Threshold is %s") % threshold)
+ """
+ Creating final map of the change, using a statistical threshold
+ """
+ change_map_calculation(changemap_name, magnitudemap_name, threshold, anglemap_class)
+ else:
+ grass.message(_("No threshold given, only angle and magnitude maps have been created"))
+
+
+ return 0
+
+def cleanup():
+ """!Delete temporary maps"""
+ TMPRAST.reverse()
+ for i in TMPRAST:
+ script.run_command("g.remove", flags='f', type='raster', name=i, quiet=True)
+
+if __name__ == "__main__":
+ atexit.register(cleanup)
+ sys.exit(main())
+
Property changes on: grass-addons/grass7/imagery/i.cva/i.cva.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/imagery/i.cva/i.cva_1angle.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.cva/i.cva_1angle.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: grass-addons/grass7/imagery/i.cva/i.cva_2angleclass.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.cva/i.cva_2angleclass.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: grass-addons/grass7/imagery/i.cva/i.cva_3magnitude.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.cva/i.cva_3magnitude.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: grass-addons/grass7/imagery/i.cva/i.cva_4change.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.cva/i.cva_4change.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
More information about the grass-commit
mailing list