[GRASS-SVN] r69852 - grass-addons/grass7/raster/r.forestfrag
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Nov 19 19:51:42 PST 2016
Author: wenzeslaus
Date: 2016-11-19 19:51:42 -0800 (Sat, 19 Nov 2016)
New Revision: 69852
Modified:
grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py
Log:
r.forestfrag: prepare for 3D, rename non-forest
Modified: grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py
===================================================================
--- grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py 2016-11-20 00:54:05 UTC (rev 69851)
+++ grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py 2016-11-20 03:51:42 UTC (rev 69852)
@@ -5,8 +5,10 @@
#
# MODULE: r.forestfrag
#
-# AUTHOR(S): Emmanuel Sambale, Stefan Sylla,
+# AUTHOR(S): Emmanuel Sambale (original shell version)
+# Stefan Sylla (original shell version)
# Paulo van Breugel (Python version, paulo at ecodiv.org)
+# Vaclav Petras (major code clean up, wenzeslaus gmail com)
#
# PURPOSE: Creates forest fragmentation index map from a
# forest-non-forest raster; The index map is based on
@@ -93,7 +95,6 @@
#% required: no
#%end
-# import libraries
import os
import sys
import uuid
@@ -101,10 +102,12 @@
import tempfile
import string
import grass.script as gs
+# neutral naming for better compatibility between 2D and 3D version
+from grass.script.raster import mapcalc
LABELS = """\
-0 nonforest
+0 exterior
1 patch
2 transitional
3 edge
@@ -123,7 +126,6 @@
6 145:207:96
"""
-
# create set to store names of temporary maps to be deleted upon exit
CLEAN_RAST = []
@@ -159,7 +161,7 @@
"""Generate window (matrix) expression
:param map_name: name of the map to index
- :param max_index: the maximum positive index to use
+ :param max_index: the maximum positive index to use;
usually (window_size - 1) / 2
:param combine_op: operator used to combine values in the pair
:param aggregate_op: operator used to combine all pairs together
@@ -218,10 +220,10 @@
# set to current input map region if requested by the user
# default is (and should be) the current region
# we could use tmp region for this but if the flag is there
- # it makes sense to use it from now on
+ # it makes sense to use it from now on (but we should reconsider)
if flag_r:
gs.message(_("Setting region to input map..."))
- gs.run_command('g.region', quiet=True, raster=ipl)
+ gs.run_command('g.region', raster=ipl, quiet=True)
# check if map values are limited to 1 and 0
input_info = gs.raster_info(ipl)
@@ -237,11 +239,11 @@
% (input_info['min'], input_info['max']))
# computing pf values
- gs.info(_("Step 1: Computing Pf values..."))
-
# let forested pixels be x and number of all pixels in moving window
# be y, then pf=x/y"
+ gs.info(_("Step 1: Computing Pf values..."))
+
# generate grid with pixel-value=number of forest-pixels in window
# generate grid with pixel-value=number of pixels in moving window:
tmpA2 = tmpname('tmpA01_')
@@ -254,17 +256,17 @@
pf = user_pf
else:
pf = tmpname('tmpA03_')
- gs.mapcalc("$pf = if( $ipl >=0, float($tmpA2) / float($tmpC3))",
- ipl=ipl, pf=pf, tmpA2=tmpA2, tmpC3=tmpC3)
+ mapcalc("$pf = if( $ipl >=0, float($tmpA2) / float($tmpC3))",
+ ipl=ipl, pf=pf, tmpA2=tmpA2, tmpC3=tmpC3)
# computing pff values
- gs.info(_("Step 2: Computing Pff values..."))
-
# Considering pairs of pixels in cardinal directions in
# a 3x3 window, the total number of adjacent pixel pairs is 12.
# Assuming that x pairs include at least one forested pixel, and
# y of those pairs are forest-forest pairs, so pff equals y/x.
+ gs.info(_("Step 2: Computing Pff values..."))
+
# create copy of forest map and convert NULL to 0 (if any)
tmpC4 = tmpname('tmpA04_')
gs.run_command("g.copy", raster=[ipl, tmpC4], quiet=True)
@@ -284,15 +286,15 @@
else:
pff = tmpname('tmpA07_')
# potentially this can be split and parallelized
- gs.mapcalc("$pff = if($ipl >= 0, float($tmpl4) / float($tmpl5))",
- ipl=ipl, tmpl4=expr1, tmpl5=expr2, pff=pff)
+ mapcalc("$pff = if($ipl >= 0, float($tmpl4) / float($tmpl5))",
+ ipl=ipl, tmpl4=expr1, tmpl5=expr2, pff=pff)
# computing fragmentation index
# (a b) name, condition
# where a is a number used by Riitters et al. in ERRATUM (2)
# and b is a number used in the sh script by Sambale and Sylla
# b also defines 0 for non-forested which is consistent with input
- # (1 1) edge, if Pf > 0.6 and Pf - Pff < 0
+ # (1 3) edge, if Pf > 0.6 and Pf - Pff < 0
# (2 6) undetermined, if Pf > 0.6 and Pf = Pff
# (3 4) perforated, if Pf > 0.6 and Pf - Pff > 0
# (4 5) interior, if Pf = 1.0
@@ -305,7 +307,7 @@
indexfin2 = tmpname('tmpA16_')
else:
indexfin2 = opl
- gs.mapcalc(
+ mapcalc(
"eval("
"dpf = $pf - $pff,"
# individual classes
@@ -344,7 +346,7 @@
e=float(reginfo['e']) - ewcor,
w=float(reginfo['w']) + ewcor,
quiet=True)
- gs.mapcalc("$opl = $if3", opl=opl, if3=indexfin2, quiet=True)
+ mapcalc("$opl = $if3", opl=opl, if3=indexfin2, quiet=True)
# create categories
# TODO: parametrize classes (also in r.mapcalc, r.colors and desc)?
@@ -404,12 +406,13 @@
if flag_s:
gs.run_command("r.report", map=opl, units=["h", "p"],
flags="n", page_width=50, quiet=True)
+
gs.info(_("The following layers were created"))
gs.info(_("The fragmentation index: %s") % opl)
if user_pf:
- gs.info(_("The proportion forested (pf): %s") % pf)
+ gs.info(_("The proportion forested (Pf): %s") % pf)
if user_pff:
- gs.info(_("The proportion forested pixel pairs (pff): %s") % pff)
+ gs.info(_("The proportion forested pixel pairs (Pff): %s") % pff)
if __name__ == "__main__":
More information about the grass-commit
mailing list