[GRASS-SVN] r60180 - grass-addons/grass7/vector/v.habitat.dem
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun May 11 09:01:50 PDT 2014
Author: hellik
Date: 2014-05-11 09:01:50 -0700 (Sun, 11 May 2014)
New Revision: 60180
Modified:
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:
svn propset
Modified: grass-addons/grass7/vector/v.habitat.dem/Makefile
===================================================================
--- grass-addons/grass7/vector/v.habitat.dem/Makefile 2014-05-11 15:52:33 UTC (rev 60179)
+++ grass-addons/grass7/vector/v.habitat.dem/Makefile 2014-05-11 16:01:50 UTC (rev 60180)
@@ -1,7 +1,7 @@
-MODULE_TOPDIR = ../..
-
-PGM = v.habitat.dem
-
-include $(MODULE_TOPDIR)/include/Make/Script.make
-
-default: script
+MODULE_TOPDIR = ../..
+
+PGM = v.habitat.dem
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/vector/v.habitat.dem/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.html
===================================================================
--- grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.html 2014-05-11 15:52:33 UTC (rev 60179)
+++ grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.html 2014-05-11 16:01:50 UTC (rev 60180)
@@ -1,210 +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>
+<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$</i>
</p>
\ No newline at end of file
Property changes on: grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.py
===================================================================
--- grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.py 2014-05-11 15:52:33 UTC (rev 60179)
+++ grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.py 2014-05-11 16:01:50 UTC (rev 60180)
@@ -1,519 +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())
+#!/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())
Property changes on: grass-addons/grass7/vector/v.habitat.dem/v.habitat.dem.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list