[GRASS-SVN] r60179 - grass-addons/grass7/vector/v.habitat.dem
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun May 11 08:52:33 PDT 2014
Author: hellik
Date: 2014-05-11 08:52:33 -0700 (Sun, 11 May 2014)
New Revision: 60179
Added:
grass-addons/grass7/vector/v.habitat.dem/Makefile
grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.html
grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.py
Log:
v.habitat.dem - DEM and solar derived characteristics of habitats
Added: grass-addons/grass7/vector/v.habitat.dem/Makefile
===================================================================
--- grass-addons/grass7/vector/v.habitat.dem/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.habitat.dem/Makefile 2014-05-11 15:52:33 UTC (rev 60179)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.habitat.dem
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.html
===================================================================
--- grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.html (rev 0)
+++ grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.html 2014-05-11 15:52:33 UTC (rev 60179)
@@ -0,0 +1,210 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.habitat.dem</em> calculates DEM and solar derived characteristics of habitat
+vector polygons. The user must specify the input <b>elevation raster</b> map,
+a <b>habitat vector</b> map with a <b>numeric unique ID</b> column and a <b>prefix</b>
+used for all results.
+
+<p>
+A preliminary visual delineation of habitats based upon digital orthophotos is a
+common task for an ecologist before fieldwork. Ecological site conditions of habitats
+are often influenced amongst others by terrain forms, solar irradiance and irradiation.
+<em>v.habitat.dem</em> gives some DEM derived characteristics for a quick validation
+of the preliminary visual habitat delineation.
+</p>
+
+<h2>NOTES</h2>
+
+The location has to be in a projected coordination system. Before running <em>v.habitat.dem</em>
+the region has to be aligned to the <b>elevation raster</b> map and the <b>habitat vector</b> map by
+<em>g.region</em>. During calculations, especially for solar characteristics, the region will be
+extended by a user input (default 5.000m). The results are as good as the DEM quality and resolution is.
+
+
+<h3>Terrain characteristics</h3>
+
+<b>Slope</b> and <b>aspect</b> are calculated by <em>r.slope.aspect</em>.
+
+<p>
+The <b>slope</b> output raster map contains slope values, stated in degrees of
+inclination from the horizontal.
+<br>
+The <b>aspect</b> output raster map indicates the direction that slopes are
+facing. The aspect categories represent the number degrees of east.
+</p>
+
+<b>Accumulation</b>, <b>drainage direction</b> and <b>topographic index</b> are
+calculated by <em>r.watershed</em>. The flag <b>-a</b> (use positive flow accumulation
+even for likely underestimates) is used as default.
+
+<p>
+The <b>accumulation</b> map contains the absolute value of each cell in this output
+map and is the amount of overland flow that traverses the cell. This value will be
+the number of upland cells plus one if no overland flow map is given.
+<br>
+The <b>drainage direction</b> map contains drainage direction. Provides the "aspect"
+for each cell measured CCW from East.
+<br>
+The <b>topographic index</b> raster map contains topographic index TCI and is
+computed as <tt>ln(α / tan(β))</tt> where α a is the cumulative
+upslope area draining through a point per unit contour length and <tt>tan(β)</tt>
+is the local slope angle. The TCI reflects the tendency of water to accumulate at
+any point in the catchment and the tendency for gravitaional forces to move that
+water downslope. This value will be negative if <tt>α / tan(β) < 1</tt>.
+</p>
+
+Terrain forms are calculated by <em>r.geomorphon</em>.
+
+<p>
+Geomorphon is a new concept of presentation and analysis of terrain forms using machine vision
+approach. This concept utilises 8-tuple pattern of the visibility neighbourhood and breaks
+well known limitation of standard calculus approach where all terrain forms cannot
+be detected in a single window size. The pattern arises from a comparison of a focus
+pixel with its eight neighbours starting from the one located to the east and continuing
+counterclockwise producing a ternary operator. All options in the <em>r.geomorphon</em>-calculation
+are set to default (<b>skip = 0</b>, <b>search = 3</b>, <b>flat = 1</b>, <b>dist = 0</b>) where
+<b>search</b> determines the length on the geodesic distances in all eight directions where
+line-of-sight is calculated, <b>skip</b> determines length on the geodesic distances at the
+beginning of calculation all eight directions where line-of-sight is yet calculated, <b>flat</b> defines
+the difference (in degrees) between zenith and nadir line-of-sight which indicate flat direction and <b>dist</b>
+determines > flat distance..
+</p>
+
+ The most common terrain forms calculated by <em>r.geomorphon</em> are:
+
+<ul>
+ <li>flat</li>
+ <li>summit</li>
+ <li>ridge</li>
+ <li>shoulder</li>
+ <li>spur</li>
+ <li>slope</li>
+ <li>hollow</li>
+ <li>footslope</li>
+ <li>valley</li>
+ <li>depression</li>
+</ul>
+
+The LS factor
+
+<p>
+The LS is the slope length-gradient factor. The LS factor represents a ratio of soil loss under given conditions to
+that at a site with the "standard" slope steepness of 9% and slope length of 22.13m. The steeper and longer the slope,
+the higher the risk for erosion.
+</p>
+
+The LS factor is calculated accordingly Neteler & Mitasova 2008 in <em>r.mapcalc</em> with flow accumulation of <em>r.flow</em>
+and slope of <em>r.slope.aspect</em>
+
+<div class="code">
+ <pre>
+ 1.4 * exp(flow_accumulation / 22.1, 0.4) * exp(sin(slope) 0.09, 1.2)
+ </pre>
+</div>
+
+The colors of the LS factor map are set to:
+
+<ul>
+ <li>0 white</li>
+ <li>3 yellow</li>
+ <li>6 orange</li>
+ <li>10 red</li>
+ <li>50 magenta</li>
+ <li>100 violet</li>
+</ul>
+
+<h3>Terrain characteristics uploaded to the habitat vector attribute table per polygon</h3>
+
+<ul>
+ <li><b>DEM altitude</b>: minimum, maximum, range, average, median</li>
+ <li><b>slope</b>: minimum, maximum, range, average, median</li>
+ <li><b>aspect</b>: minimum, maximum, range, average, median</li>
+ <li><b>geomorphons</b>: absolute area of flat, summit, ridge, shoulder, spur, slope, hollow, footslope, valley, depression</li>
+</ul>
+
+Additionally the mutual occurrence by <em>r.coin</em> of unique habitat ID and geomorphons in percent of the row is printed to the output.
+
+<h3>Simple check of terrain characteristics</h3>
+
+Simple checks regarding aspect and slope per unique habitat ID are evaluated and marked in the attribute table as follow:
+
+<ul>
+<li><b>simple check regarding aspect range:</b></li>
+<li>aspect range 100-200 *</li>
+<li>aspect range 201-300 **</li>
+<li>aspect range ≥ 300 ***</li>
+<li><b>simple checks regarding aspect range and slope median:</b></li>
+<li>aspect range 100-200 and median slope < 5 *</li>
+<li>aspect range 201-300 and median slope < 5 **</li>
+<li>aspect range ≥ 300 and median slope < 5 ***</li>
+</ul>
+
+These simple checks may indicate reconsidering of some preliminary visual habitat delineations.
+
+<h3>Solar characteristics</h3>
+
+<p>
+The solar characterstics (direct sunlight / shadows caused by terrain for a certain day in the year) are
+calculated by <em>r.sun.hourly</em> based upon <em>r.sun</em>. The <b>-b</b>-flag is used to create binary rasters
+instead of irradiation rasters. The user can define start time of interval, end time of interval, time step for
+running <em>r.sun</em>, number of day of the year and the year. As default is set summer solstice (21st June 2014, 8:00-18:00, 1 hour time step).
+</p>
+
+The results of the <em>r.sun.hourly</em>-analysis are automatically registered into a temporal database. The space time raster dataset
+can be easily animated in the <em>g.gui.animation</em>-tool.
+
+<h2>EXAMPLE</h2>
+
+<div class="code">
+ <pre>
+ # align region to DEM and habitat vector
+ g.region -a rast=DEM vect=myhabitats align=DEM
+
+ # run <em>v.habitat.dem</em>
+ v.habitat.dem elevation=DEM vector=myhabitats column=Id prefix=a dir=C:\wd
+
+ # do <em>r.null</em> to the <em>r.sun.hourly</em> output to get maps without direct beam
+ r.null map=a_beam_rad_08.00 setnull=1
+ [...]
+ r.null map=a_beam_rad_18.00 setnull=1
+
+ # animate the <em>r.sun.hourly</em> output by the <em>g.gui.animation</em>-tool
+ g.gui.animation strds=a_beam_rad
+ </pre>
+</div>
+
+<h2>DEPENDENCIES</h2>
+
+<ul>
+<li>r.geomorphon (addon)</li>
+<li>r.sun.hourly (addon)</li>
+</ul>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="g.gui.animation.html">g.gui.animation</a>
+<a href="g.region.html">g.region</a>,
+<a href="r.coin.html">r.coin</a>,
+<a href="r.geomorphon.html">r.geomorphon</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
+<a href="r.slope.aspect.html">r.slope.aspect</a>,
+<a href="r.sun.html">r.sun</a>,
+<a href="r.sun.hourly.html">r.sun.hourly</a>,
+<a href="r.stats.html">r.stats</a>,
+<a href="r.watershed.html">r.watershed</a>,
+<a href="v.rast.stats.html">v.rast.stats</a>,
+<a href="v.to.rast.html">v.to.rast</a>
+</em>
+
+<h2>REFERENCES</h2>
+
+Neteler, M. & Mitasova, H. 2008. Open Source GIS. A GRASS GIS Aproach. Third Edition. Springer.
+
+<h2>AUTHOR</h2>
+
+Helmut Kudrnovsky
+
+<p>
+<i>Last changed: $Date: 2014-05-11 12:00:00 +0200 (So, 11 May 2014) $</i>
+</p>
\ No newline at end of file
Added: grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.py
===================================================================
--- grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.py (rev 0)
+++ grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.py 2014-05-11 15:52:33 UTC (rev 60179)
@@ -0,0 +1,519 @@
+#!/usr/bin/env python
+
+"""
+MODULE: v.habitat.dem
+
+AUTHOR(S): Helmut Kudrnovsky <alectoria AT gmx at>
+
+PURPOSE: DEM and solar derived characteristics of habitats
+
+COPYRIGHT: (C) 2014 by the GRASS Development Team
+
+ This program is free software under the GNU General Public
+ License (>=v2). Read the file COPYING that comes with GRASS
+ for details.
+"""
+
+#%module
+#% description: DEM derived characteristics of habitats
+#% keywords: vector
+#% keywords: raster
+#% keywords: terrain
+#% keywords: statistics
+#% keywords: sun
+#% keywords: zonal statistics
+#%end
+
+#%option G_OPT_R_ELEV
+#% key: elevation
+#% description: Name of elevation raster map
+#% required: yes
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: vector
+#% description: Name of habitat vector map
+#% required: yes
+#%end
+
+#%option G_OPT_DB_COLUMN
+#% description: Name of attribute column with a unique habitat ID (must be numeric)
+#% required: yes
+#%end
+
+#%option
+#% key: prefix
+#% type: string
+#% key_desc: prefix
+#% description: output prefix (must start with a letter)
+#% required: yes
+#%end
+
+#%option G_OPT_M_DIR
+#% key: dir
+#% description: Directory where the output will be found
+#% required : yes
+#%end
+
+#%option
+#% key: region_extension
+#% type: double
+#% key_desc: float
+#% description: region extension
+#% required : no
+#% answer: 5000
+#%end
+
+#%option
+#% key: start_time
+#% type: double
+#% label: Start time of interval
+#% description: Use up to 2 decimal places
+#% options: 0-24
+#% answer: 8
+#%end
+
+#%option
+#% key: end_time
+#% type: double
+#% label: End time of interval
+#% description: Use up to 2 decimal places
+#% options: 0-24
+#% answer: 18
+#%end
+
+#%option
+#% key: time_step
+#% type: double
+#% label: Time step for running r.sun [decimal hours]
+#% description: Use up to 2 decimal places
+#% options: 0-24
+#% answer: 1
+#%end
+
+#%option
+#% key: day
+#% type: integer
+#% description: No. of day of the year
+#% options: 1-365
+#% answer: 172
+#%end
+
+#%option
+#% key: year
+#% type: integer
+#% label: Year used for map registration into temporal dataset or r.timestamp
+#% description: This value is not used in r.sun calculations
+#% options: 1900-9999
+#% answer: 2014
+#%end
+
+import sys
+import os
+import grass.script as grass
+import math
+from numpy import array
+from numpy import zeros
+import csv
+
+if not os.environ.has_key("GISBASE"):
+ grass.message( "You must be in GRASS GIS to run this program." )
+ sys.exit(1)
+
+def main():
+ r_elevation = options['elevation'].split('@')[0]
+ v_habitat = options['vector'].split('@')[0]
+ v_column = options['column']
+ regext = options['region_extension']
+ directory = options['dir']
+ prefix = options['prefix']
+ r_accumulation = prefix+'_accumulation'
+ r_drainage = prefix+'_drainage'
+ r_tci = prefix+'_topographicindex'
+ r_slope = prefix+'_slope'
+ r_aspect = prefix+'_aspect'
+ r_flow_accum = prefix+'_flow_accumulation'
+ r_LS = prefix+'_LS_factor'
+ r_habitat = prefix+'_'+v_habitat+'_rast'
+ r_geomorphon = prefix+'_'+r_elevation+'_geomorphon'
+ r_beam_rad = prefix+'_beam_rad'
+ f_habitat_small_areas = prefix+'_small_areas.csv'
+ f_habitat_geomorphons = prefix+'_habitat_geomorphons.csv'
+ f_habitat_geomorphons_pivot = prefix+'_habitat_geomorphons_pivot.csv'
+ f_LS_colorrules = prefix+'_LS_color.txt'
+ t_habitat_geomorphons = prefix+'_habgeo_tab'
+ t_habitat_geomorphons_pivot = prefix+'_habgeo_tab_pivot'
+ d_start_time = options['start_time']
+ d_end_time = options['end_time']
+ d_time_step = options['time_step']
+ d_day = options['day']
+ d_year = options['year']
+ saved_region = prefix+'_region_ori'
+ current_region = grass.region()
+ X = float(regext)
+ N = current_region["n"]
+ S = current_region["s"]
+ E = current_region["e"]
+ W = current_region["w"]
+ Yres = current_region['nsres']
+ Xres = current_region['ewres']
+ PixelArea = Xres * Yres
+ global tmp
+
+ ## check if r.geomorphon addon is installed
+ if not grass.find_program('r.geomorphon', '--help'):
+ grass.fatal(_("The 'r.geomorphon' addon was not found, install it first:") +
+ "\n" +
+ "g.extension r.geomorphon")
+
+ ## check if r.geomorphon addon is installed
+ if not grass.find_program('r.sun.hourly', '--help'):
+ grass.fatal(_("The 'r.sun.hourly' addon was not found, install it first:") +
+ "\n" +
+ "g.extension r.sun.hourly")
+
+ # Region settings
+ grass.message( "Current region will be saved and extended." )
+ grass.message( "----" )
+
+ # Print and save current region
+ grass.run_command('g.region', flags = 'p', save = saved_region, overwrite = True)
+ grass.message( "Current region saved." )
+ grass.message( "----" )
+
+ # Align region to elevation raster and habitat vector
+ # grass.message( "Align region to elevation raster and habitat vector ..." )
+ # grass.run_command('g.region', flags = 'ap',
+ # rast = r_elevation,
+ # vect = v_habitat,
+ # align = r_elevation)
+ #grass.message( "Alignment done." )
+ #grass.message( "----" )
+
+ # Extend region
+ grass.message( "Extend region ..." )
+ grass.message( "n, s, e, w" )
+ grass.message( [current_region[key] for key in "nsew"] )
+ grass.message( "by" )
+ grass.message( regext )
+ grass.message( "in all directions" )
+
+ grass.run_command('g.region', n = N+X,
+ s = S-X,
+ e = E+X,
+ w = W-X)
+ grass.message( "Region extension done." )
+
+ grass.message( "Extended region:" )
+ grass.run_command('g.region', flags = 'p')
+ grass.message( "----" )
+
+ # Watershed calculation: accumulation, drainage direction, topographic index
+ grass.message( "Calculation of accumulation, drainage direction, topographic index by r.watershed ..." )
+ grass.run_command('r.watershed', elevation = r_elevation,
+ accumulation = r_accumulation,
+ drainage = r_drainage,
+ tci = r_tci,
+ convergence = 5,
+ flags = 'am',
+ overwrite = True)
+ grass.message( "Calculation of accumulation, drainage direction, topographic done." )
+ grass.message( "----" )
+
+ # Calculation of slope and aspect maps
+ grass.message( "Calculation of slope and aspect by r.slope.aspect ..." )
+ grass.run_command('r.slope.aspect', elevation = r_elevation,
+ slope = r_slope,
+ aspect = r_aspect,
+ overwrite = True)
+ grass.message( "Calculation of slope and aspect done." )
+ grass.message( "----" )
+
+ # Calculate pixel area by nsres x ewres
+ grass.message( "Pixel area:" )
+ grass.message( PixelArea )
+ grass.message( "----" )
+
+ # Calculate habitat area and populate it to the attribute table
+ grass.message( "Calculate habitat's areas and populate it to the attribute table ..." )
+ grass.run_command("v.db.addcolumn", map = v_habitat,
+ layer = 1,
+ columns = "habarea double")
+
+ grass.run_command("v.to.db", map = v_habitat,
+ option = 'area',
+ layer = 1,
+ columns = 'habarea',
+ overwrite = True)
+
+ grass.message( "Calculate habitat's areas done." )
+ grass.message( "----" )
+
+ # Show habitat areas smaller than Pixel Area
+ grass.message( "Habitat areas smaller than pixel area." )
+ grass.run_command("v.db.select", map = v_habitat,
+ flags = 'v',
+ layer = 1,
+ columns = v_column,
+ where = "habarea < %s" % (PixelArea))
+
+ smallareacsv = os.path.join( directory, f_habitat_small_areas )
+
+ grass.run_command("v.db.select", map = v_habitat,
+ flags = 'v',
+ layer = 1,
+ columns = v_column,
+ file = smallareacsv,
+ where = "habarea < %s" % (PixelArea))
+
+ grass.message( "A list of habitat areas smaller than pixel area can be found in: " )
+ grass.message( smallareacsv )
+ grass.message( "----" )
+
+ # Mark habitats smaller than pixel area in attribute table
+ grass.message( "Mark habitats smaller than pixel area in attribute table ..." )
+ grass.run_command("v.db.addcolumn", map = v_habitat,
+ layer = 1,
+ columns = "%s_smallarea varchar(1)" % prefix )
+ grass.run_command("v.db.update", map = v_habitat, layer = 1, column = '%s_smallarea' % (prefix), value = '*', where = 'habarea < %s' % (PixelArea))
+ grass.message( "See column" )
+ grass.message( '%s_smallarea' % prefix )
+ grass.message( "marked by *." )
+ grass.message( "----" )
+
+ # Upload DEM zonal statistics to the attribute table
+ grass.message( "Upload DEM zonal statistics to the attribute table ..." )
+ grass.run_command("v.rast.stats", map = v_habitat,
+ flags = 'c',
+ layer = 1,
+ raster = r_elevation,
+ column_prefix = prefix+'_dem',
+ method = 'minimum,maximum,range,average,median')
+ grass.message( "Upload DEM zonal statistics done." )
+ grass.message( "----" )
+
+ # Upload slope zonal statistics to the attribute table
+ grass.message( "Upload slope zonal statistics to the attribute table ..." )
+ grass.run_command("v.rast.stats", map = v_habitat,
+ flags = 'c',
+ layer = 1,
+ raster = r_slope,
+ column_prefix = prefix+'_slope',
+ method = 'minimum,maximum,range,average,median')
+ grass.message( "Upload slope zonal statistics done." )
+ grass.message( "----" )
+
+ # Upload slope zonal statistics to the attribute table
+ grass.message( "Upload aspect zonal statistics to the attribute table ..." )
+ grass.run_command("v.rast.stats", map = v_habitat,
+ flags = 'c',
+ layer = 1,
+ raster = r_aspect,
+ column_prefix = prefix+'_aspect',
+ method = 'minimum,maximum,range,average,median')
+ grass.message( "Upload aspect zonal statistics done." )
+ grass.message( "----" )
+
+ # Do some simple checks regarding aspect range
+ grass.message( "Do some simple checks regarding aspect range and populate it to the attribute table..." )
+ grass.message( "aspect range 100-200 *" )
+ grass.message( "aspect range 201-300 **" )
+ grass.message( "aspect range >= 300 ***" )
+ grass.run_command("v.db.addcolumn", map = v_habitat,
+ layer = 1,
+ columns = "%s varchar(3)" % (prefix+'_check_aspect_range'))
+
+ grass.run_command("db.execute", sql = "UPDATE %s SET %s ='*' WHERE %s < 200 AND %s >= 100" % (v_habitat, prefix+'_check_aspect_range', prefix+'_aspect_range', prefix+'_aspect_range'))
+ grass.run_command("db.execute", sql = "UPDATE %s SET %s ='**' WHERE %s < 300 AND %s >= 200" % (v_habitat, prefix+'_check_aspect_range', prefix+'_aspect_range', prefix+'_aspect_range'))
+ grass.run_command("db.execute", sql = "UPDATE %s SET %s ='***' WHERE %s >= 300" % (v_habitat, prefix+'_check_aspect_range', prefix+'_aspect_range'))
+
+
+ grass.message( "Simple checks regarding aspect range done." )
+ grass.message( "----" )
+
+ # Do some simple checks regarding aspect and and slope
+ grass.message( "Do some simple checks regarding aspect range and slope median and populate it to the attribute table..." )
+ grass.message( "aspect range 100-200 and median slope < 5 *" )
+ grass.message( "aspect range 201-300 and median slope < 5 **" )
+ grass.message( "aspect range >= 300 and median slope < 5 ***" )
+ grass.run_command("v.db.addcolumn", map = v_habitat,
+ layer = 1,
+ columns = "%s varchar(3)" % (prefix+'_check_aspect_slope'))
+
+ grass.run_command("db.execute", sql = "UPDATE %s SET %s ='*' WHERE (%s < 200 AND %s >= 100) AND %s < 5" % (v_habitat, prefix+'_check_aspect_slope', prefix+'_aspect_range', prefix+'_aspect_range', prefix+'_slope_median'))
+ grass.run_command("db.execute", sql = "UPDATE %s SET %s ='**' WHERE (%s < 300 AND %s >= 200) AND %s < 5" % (v_habitat, prefix+'_check_aspect_slope', prefix+'_aspect_range', prefix+'_aspect_range', prefix+'_slope_median'))
+ grass.run_command("db.execute", sql = "UPDATE %s SET %s ='***' WHERE %s >= 300 AND %s < 5" % (v_habitat, prefix+'_check_aspect_slope', prefix+'_aspect_range', prefix+'_slope_median'))
+
+ grass.message( "Simple checks regarding aspect range and median slope done." )
+ grass.message( "----" )
+
+ # Do some simple checks regarding aspect and and slope
+ grass.message( "Convert habitat vector to raster ..." )
+ grass.run_command("v.to.rast", input = v_habitat,
+ layer = 1,
+ type = 'area',
+ use = 'attr',
+ attrcolumn = v_column,
+ output = r_habitat)
+ grass.message( "Conversion done." )
+ grass.message( "----" )
+
+ # Calculate the most common geomorphons
+ grass.message( "Calculate the most common geomorphons ..." )
+ grass.run_command("r.geomorphon", dem = r_elevation,
+ skip = 0,
+ search = 3,
+ flat = 1,
+ dist = 0,
+ forms = r_geomorphon)
+ grass.message( "Geomorphon calculations done." )
+ grass.message( "----" )
+
+ # Mutual occurrence (coincidence) of categories of habitats and geomorphons
+ grass.message( "Calculate mutual occurrences of habitats and geomorphons ..." )
+ grass.message( "1 - flat" )
+ grass.message( "2 - summit" )
+ grass.message( "3 - ridge" )
+ grass.message( "4 - shoulder" )
+ grass.message( "5 - spur" )
+ grass.message( "6 - slope" )
+ grass.message( "7 - hollow" )
+ grass.message( "8 - footslope" )
+ grass.message( "9 - valley" )
+ grass.message( "10 - depression" )
+ grass.message( " " )
+ grass.message( "Mutual occurrence in percent of the row" )
+ grass.run_command("r.coin", map1 = r_habitat,
+ map2 = r_geomorphon,
+ flags = 'w',
+ units = 'y')
+ grass.message( "Calculations of mutual occurrences done." )
+ grass.message( "----" )
+
+ # Join geomorphons to habitat attribute table
+ grass.message( "Join geomorphon information to habitat attribute table ...." )
+
+ habgeocsv = os.path.join( directory, f_habitat_geomorphons)
+
+ grass.run_command("r.stats", input = [r_habitat, r_geomorphon],
+ flags = 'aln',
+ separator = ';',
+ output = habgeocsv)
+
+ grass.run_command("db.in.ogr", dsn = habgeocsv,
+ output = t_habitat_geomorphons)
+
+ grass.run_command("db.dropcolumn", table = t_habitat_geomorphons,
+ column = 'field_2',
+ flags = 'f')
+
+ grass.run_command("db.dropcolumn", table = t_habitat_geomorphons,
+ column = 'field_3',
+ flags = 'f')
+
+ habgeocsv_pivot = os.path.join( directory, f_habitat_geomorphons_pivot)
+
+ grass.run_command("db.select", separator = ';',
+ output = habgeocsv_pivot,
+ sql = "SELECT field_1, sum(case when field_4 = 'flat' then field_5 end) as flat, sum(case when field_4 = 'summit' then field_5 end) as summit, sum(case when field_4 = 'ridge' then field_5 end) as ridge, sum(case when field_4 = 'shoulder' then field_5 end) as shoulder, sum(case when field_4 = 'spur' then field_5 end) as spur, sum(case when field_4 = 'slope' then field_5 end) as slope, sum(case when field_4 = 'hollow' then field_5 end) as hollow, sum(case when field_4 = 'footslope' then field_5 end) as footslope, sum(case when field_4 = 'valley' then field_5 end) as valley, sum(case when field_4 = 'depression' then field_5 end) as depression , sum(field_5) as SubTotal FROM %s GROUP BY field_1" % t_habitat_geomorphons)
+
+ grass.run_command("db.in.ogr", dsn = habgeocsv_pivot,
+ output = t_habitat_geomorphons_pivot)
+
+ grass.run_command("v.db.join", map = v_habitat,
+ column = v_column,
+ otable = t_habitat_geomorphons_pivot,
+ ocolumn = 'field_1')
+ grass.message( "Geomorphon information joint to habitat attribute table." )
+ grass.message( "----" )
+
+ # Give information where output files are
+ grass.message( "Geomorphon information:" )
+ grass.message( habgeocsv )
+ grass.message( "Geomorphon information in pivot format:" )
+ grass.message( habgeocsv_pivot )
+ grass.message( "----" )
+
+ # Calculate LS factor see Neteler & Mitasova 2008. Open Source GIS - A GRASS GIS Approach
+ grass.message( "Calculate LS factor ..." )
+ grass.run_command("r.flow", elevation = r_elevation,
+ aspect = r_aspect,
+ flowaccumulation = r_flow_accum)
+
+ grass.message( "..." )
+ grass.mapcalc("$outmap = 1.4 * exp($flowacc / 22.1, 0.4) * exp(sin($slope) / 0.09, 1.2)",
+ outmap = r_LS,
+ flowacc = r_flow_accum,
+ slope = r_slope)
+
+ # create and define color rules file for LS factor map
+ ls_color_rules_out = os.path.join( directory, f_LS_colorrules )
+ with open(ls_color_rules_out, 'w') as f:
+ writer = csv.writer(f)
+ writer.writerow(['0 white'])
+ writer.writerow(['3 yellow'])
+ writer.writerow(['6 orange'])
+ writer.writerow(['10 red'])
+ writer.writerow(['50 magenta'])
+ writer.writerow(['100 violet'])
+
+ grass.run_command("r.colors", map = r_LS,
+ rules = ls_color_rules_out)
+
+ grass.message( "Calculation LS factor done." )
+ grass.message( "----" )
+
+ # Run r.sun.hourly in binary mode for light/shadow
+ grass.message( "Run r.sun.hourly in binary mode for light/shadow for a certain day in the year ..." )
+ grass.run_command("r.sun.hourly", elev_in = r_elevation,
+ flags = 'tb',
+ asp_in = r_aspect,
+ slope_in = r_slope,
+ start_time = d_start_time,
+ end_time = d_end_time,
+ day = d_day,
+ year = d_year,
+ beam_rad_basename = r_beam_rad)
+
+ grass.message( "----" )
+ grass.message( "Light/shadow conditions calculated for year" )
+ grass.message( d_year )
+ grass.message( "and day" )
+ grass.message( d_day )
+ grass.message( 'from' )
+ grass.message( d_start_time )
+ grass.message( 'to' )
+ grass.message( d_end_time )
+ grass.message( 'done.' )
+ grass.message( "----" )
+ grass.run_command("t.info", flags = 'h',
+ input = r_beam_rad)
+ grass.message( "----" )
+
+ # Set region to original
+ grass.message( "Restore original region settings:" )
+ grass.run_command("g.region", flags = 'p',
+ region = saved_region)
+ grass.message( "----" )
+
+ # clean up some temporay files and maps
+ grass.message( "Some clean up ..." )
+ grass.run_command("g.remove", region = saved_region)
+ grass.run_command("g.remove", rast = r_flow_accum)
+ grass.run_command("g.remove", rast = r_habitat)
+ grass.run_command("db.droptable", flags = 'f',
+ table = t_habitat_geomorphons)
+ grass.run_command("db.droptable", flags = 'f',
+ table = t_habitat_geomorphons_pivot)
+ grass.run_command("db.dropcolumn", flags = 'f',
+ table = v_habitat,
+ column = 'field_1')
+ grass.message( "Clean up done." )
+ grass.message( "----" )
+
+ # v.habitat.dem done!
+ grass.message( "v.habitat.dem done!" )
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ sys.exit(main())
More information about the grass-commit
mailing list