[GRASS-SVN] r70073 - grass-addons/grass7/imagery/i.cva
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Dec 12 13:04:42 PST 2016
Author: ciupava
Date: 2016-12-12 13:04:41 -0800 (Mon, 12 Dec 2016)
New Revision: 70073
Modified:
grass-addons/grass7/imagery/i.cva/i.cva.py
Log:
Fixed Python conventions
Modified: grass-addons/grass7/imagery/i.cva/i.cva.py
===================================================================
--- grass-addons/grass7/imagery/i.cva/i.cva.py 2016-12-12 20:16:39 UTC (rev 70072)
+++ grass-addons/grass7/imagery/i.cva/i.cva.py 2016-12-12 21:04:41 UTC (rev 70073)
@@ -47,138 +47,131 @@
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)
+ """
+ Calculates 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)
+ """
+ Calculates the vector angle as the arctg of deltaY/deltaX
+ """
+ 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)
+ """
+ Calculates the vector length (magnitude) as sqrt((deltaX)^2+(deltaY)^2)
+ """
+ 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)
+ """
+ Generates the change map as the values of the classified angle map whose magnitude follows the criterion (values higher than the threshold)
+ """
+ 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'
+ 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)
+ # Checking that the input maps exist
+ 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)
+ 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 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)
+ #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)
+ 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())
+ # 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')
+ # Generating 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 the 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!")
+ 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"))
+ adding_value = float(stat_threshold) * float(stddev)
+ threshold = float(mean) + float(adding_value)
+ grass.message(_("Threshold is %s") % threshold)
+ #Creating the 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
+ 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)
+ # !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())
+ atexit.register(cleanup)
+ sys.exit(main())
More information about the grass-commit
mailing list