[GRASS-SVN] r68142 - grass-addons/grass7/raster/r.forestfrag
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Mar 25 05:34:17 PDT 2016
Author: pvanbosgeo
Date: 2016-03-25 05:34:17 -0700 (Fri, 25 Mar 2016)
New Revision: 68142
Modified:
grass-addons/grass7/raster/r.forestfrag/r.forestfrag.html
grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py
Log:
r.forestfrag addon: better handling of edge effects / no-data cells
Modified: grass-addons/grass7/raster/r.forestfrag/r.forestfrag.html
===================================================================
--- grass-addons/grass7/raster/r.forestfrag/r.forestfrag.html 2016-03-25 08:32:44 UTC (rev 68141)
+++ grass-addons/grass7/raster/r.forestfrag/r.forestfrag.html 2016-03-25 12:34:17 UTC (rev 68142)
@@ -13,13 +13,20 @@
map refers to "between-pixel" fragmentation around the corresponding
forest location.
+<p>As input it requires a binary map with (1) forest and (0)
+non-forest. Obviously, one can replace forest any other land cover
+type. If one wants to exclude the influence of a specific land cover
+type, e.g., water bodies, it should be classified as no-data (NA) in
+the input map. See e.g.,
+<a href="https://pvanb.wordpress.com/2016/03/25/update-of-r-forestfrag-addon/"> blog post</a>.
+
<p>Let <em>Pf</em> be the proportion of pixels in the window that
are forested. Define <em>Pff</em> (strictly) as the proportion of
all adjacent (cardinal directions only) pixel pairs that include at
least one forest pixel, for which both pixels are forested. <em>Pff</em>
thus (roughly) estimates the conditional probability that,
given a pixel of forest, its neighbor is also forest. The
-classification model identifies six fragmentation categories:
+classification model then identifies six fragmentation categories as:
<div class="code"><pre>
interior: Pf = 1.0
@@ -36,26 +43,28 @@
<li>The moving window size is user-defined (default=3)
and must be an odd number. If an even number is given the function
will stop with an error message.
-<li>To avoid edge effects at the border of the raster layer the user
-can choose to have the output raster trimmed with a number of raster
-cells equal to 1/2 * the size of the moving window.
+<li>No-data cells are ignored. This means that statistics at the
+raster edges are based on fewer cells (smaller) moving windows. If this
+is a problem, the user can choose to have the output raster trimmed
+with a number of raster cells equal to 1/2 * the size of the moving
+window.
<li>The function respects the region. The user has however the option
to set the region to match the input layer.
</ul>
<h2>SEE ALSO</h2>
-This addon is based on the
-<a href="http://grasswiki.osgeo.org/wiki/AddOns/GRASS_6#r.forestfrag>">r.forestfrag.sh</a>
-script, with as extra options user-defined moving window size,
-option to trim the region (by default it respects the region).
+The addon is based on the
+<a href="http://grasswiki.osgeo.org/wiki/AddOns/GRASS_6#r.forestfrag>">
+r.forestfrag.sh</a> script, with as extra options user-defined
+moving window size, option to trim the region (by default it
+respects the region) and a better handling of no-data cells.
-<h2>AUTHORS</h2>
-<em>Emmanuel Sambale</em>: Author original script<br>
-<em>Stefan Sylla</em>: Revised / updated the r.forestfrag.sh for
+<h2>AUTHORS</h2>
+<em>Paulo van Breugel<em>: Rewritten in Python, added option to set
+moving window size, and improved handling no-data cells
+<em>Stefan Sylla</em>: Revised / updated the r.forestfrag.sh for
GRASS GIS 6<br>
-<em>Paulo van Breugel</em>: Rewritten in Python
-for use with GRASS GIS 7, and added extra options.
-Contact: paulo at ecodiv.org
+<em>Emmanuel Sambale</em>: Author original shell script<br>
<h2>REFERENCES</h2>
Riitters, K., J. Wickham, R. O'Neill, B. Jones,
Modified: grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py
===================================================================
--- grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py 2016-03-25 08:32:44 UTC (rev 68141)
+++ grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py 2016-03-25 12:34:17 UTC (rev 68142)
@@ -5,7 +5,7 @@
#
# MODULE: r.forestfrag
#
-# AUTHOR(S): Emmanuel Sambale, Stefan Sylla and Paulo van Breugel
+# AUTHOR(S): Paulo van Breugel, Emmanuel Sambale, Stefan Sylla
# Maintainer: Paulo van Breugel (paulo at ecodiv.org)
#
# PURPOSE: Creates forest fragmentation index map from a forest-non-forest
@@ -15,7 +15,7 @@
# [online] URL:http://www.consecol.org/vol4/iss2/art3/
#
#
-# COPYRIGHT: (C) 1997-2015 by the GRASS Development Team
+# COPYRIGHT: (C) 1997-2016 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
@@ -168,20 +168,18 @@
# let forested pixels be x and number of all pixels in moving window
# be y, then pf=x/y"
- # generate grid with pixel-value=number of forest-pixels in 3x3
- # moving-window:
- tmpA2 = tmpname('tmpA2_')
+ # 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_')
+ tmpC3 = tmpname('tmpA02_')
grass.run_command("r.neighbors", quiet=True, input=ipl,
- output=tmpA2, method="sum", size=wz)
+ output=[tmpA2, tmpC3], method=["sum", "count"], size=wz)
+ #grass.mapcalc("$tmpC3 = float($wz^2)",tmpC3=tmpC3a,wz=str(wz),quiet=True)
- # generate grid with pixel-value=number of pixels in moving window:
- tmpC3 = tmpname('tmpC3_')
- grass.mapcalc("$tmpC3 = float($wz^2)",tmpC3=tmpC3,wz=str(wz),quiet=True)
-
# create pf map
- pf = tmpname('pf_')
- grass.mapcalc("$pf = float($tmpA2) / float($tmpC3)", pf=pf,
- tmpA2=tmpA2, tmpC3=tmpC3,quiet=True)
+ pf = tmpname('tmpA03_')
+ grass.mapcalc("$pf = if(" + ipl + ">=0, float($tmpA2) / float($tmpC3))", pf=pf,
+ tmpA2=tmpA2, tmpC3=tmpC3, quiet=True)
#------------------------------------------------------------------------
# computing pff values
@@ -194,6 +192,11 @@
grass.info("step 2: computing pff values ...\n")
+ # Create copy of forest map and convert NULL to 0 (if any)
+ tmpC4 = tmpname('tmpA04_')
+ grass.run_command("g.copy", raster=[ipl, tmpC4], quiet=True)
+ grass.run_command("r.null", map=tmpC4, null=0, quiet=True)
+
# calculate number of 'forest-forest' pairs
#------------------------------------------------------------------------
@@ -202,7 +205,7 @@
# Write mapcalc expression to tmpf - rows
fd2, tmpfile2 = tempfile.mkstemp()
- tmpl4 = tmpname('pff1_')
+ tmpl4 = tmpname('tmpA05_')
text_file = open(tmpfile2, "w")
text_file.write(tmpl4 + " = ")
@@ -211,13 +214,13 @@
csub=range(-SWn, SWn, 1)
for k in rsub:
for l in csub:
- text_file.write(ipl + "[" + str(k) + "," + str(l) + "]*"
- + ipl + "[" + str(k) + "," + str(l+1) + "] + ")
+ text_file.write(tmpC4 + "[" + str(k) + "," + str(l) + "]*"
+ + tmpC4 + "[" + str(k) + "," + str(l+1) + "] + ")
rsub=range(-SWn, SWn+1, 1)
csub=range(SWn, -SWn, -1)
for k in rsub:
for l in csub:
- text_file.write(ipl + "[" + str(l) + "," + str(k) + "]*" + ipl
+ text_file.write(tmpC4 + "[" + str(l) + "," + str(k) + "]*" + tmpC4
+ "[" + str(l-1) + "," + str(k) + "] + ")
text_file.write("0")
text_file.close()
@@ -230,7 +233,7 @@
#------------------------------------------------------------------------
fd3, tmpfile3 = tempfile.mkstemp()
- tmpl5 = tmpname('pff2_')
+ tmpl5 = tmpname('tmpA06_')
text_file = open(tmpfile3, "w")
text_file.write(tmpl5 + " = ")
@@ -240,14 +243,14 @@
csub=range(-SWn, SWn, 1)
for k in rsub:
for l in csub:
- text_file.write("if((" + ipl + "[" + str(k) + "," + str(l) + "]+"
- + ipl + "[" + str(k) + "," + str(l+1) + "])>0,1) + ")
+ text_file.write("if((" + tmpC4 + "[" + str(k) + "," + str(l) + "]+"
+ + tmpC4 + "[" + str(k) + "," + str(l+1) + "])>0,1) + ")
rsub=range(-SWn, SWn+1, 1)
csub=range(SWn, -SWn, -1)
for k in rsub:
for l in csub:
- text_file.write("if((" + ipl + "[" + str(l) + "," + str(k) + "]+"
- + ipl + "[" + str(l-1) + "," + str(k) + "])>0,1) + ")
+ text_file.write("if((" + tmpC4 + "[" + str(l) + "," + str(k) + "]+"
+ + tmpC4 + "[" + str(l-1) + "," + str(k) + "])>0,1) + ")
text_file.write("0")
text_file.close()
@@ -258,8 +261,8 @@
# create pff map
#------------------------------------------------------------------------
- pff = tmpname('pff3_')
- grass.mapcalc("$pff = float($tmpl4) / float($tmpl5)", tmpl4=tmpl4,
+ pff = tmpname('tmpA07_')
+ grass.mapcalc("$pff = if(" + ipl + " >= 0, float($tmpl4) / float($tmpl5))", tmpl4=tmpl4,
tmpl5=tmpl5, pff=pff,quiet=True)
#------------------------------------------------------------------------
@@ -267,35 +270,29 @@
#------------------------------------------------------------------------
grass.info("step 3: computing fragmentation index ...\n")
- pf2 = tmpname('pf2_')
+ pf2 = tmpname('tmpA08_')
grass.mapcalc("$pf2 = $pf - $pff", pf2=pf2, pf=pf, pff=pff, quiet=True)
- f1 = tmpname('f1_') # patch
- grass.mapcalc("$f1 = if($pf<0.4,1,0)",
- f1=f1, pf=pf, quiet=True)
- f2 = tmpname('f2_') # transitional
- grass.mapcalc("$f2 = if($pf>=0.4 && $pf<0.6,2,0)",
- pf=pf, f2=f2, quiet=True)
- f3 = tmpname('f3_') # edge
- grass.mapcalc("$f3 = if($pf>=0.6 && $pf2<0,3,0)",
- pf=pf, pf2=pf2, f3=f3, quiet=True)
- f4 = tmpname('f4_') # perforate
- grass.mapcalc("$f4 = if($pf>0.6 && $pf2>0,4,0)",
- pf=pf, pf2=pf2, f4=f4, quiet=True)
- f5 = tmpname('f5_') # interior
- grass.mapcalc("$f5 = if($pf==1 && $pff==1,5,0)",
- pf=pf, pff=pff, f5=f5, quiet=True)
- f6 = tmpname('f6_') # undetermined
- grass.mapcalc("$f6 = if($pf>0.6 && $pf<1 && $pf2==0,6,0)",
- pf=pf, pf2=pf2, pff=pff, f6=f6, quiet=True) # undetermined
- Index = tmpname('index_')
+ f1 = tmpname('tmpA09_') # patch
+ grass.mapcalc("$f1 = if($pf<0.4,1,0)", f1=f1, pf=pf, quiet=True)
+ f2 = tmpname('tmpA10_') # transitional
+ grass.mapcalc("$f2 = if($pf>=0.4 && $pf<0.6,2,0)", pf=pf, f2=f2, quiet=True)
+ f3 = tmpname('tmpA11_') # edge
+ grass.mapcalc("$f3 = if($pf>=0.6 && $pf2<0,3,0)", pf=pf, pf2=pf2, f3=f3, quiet=True)
+ f4 = tmpname('tmpA12_') # perforate
+ grass.mapcalc("$f4 = if($pf>0.6 && $pf<1 && $pf2>0,4,0)", pf=pf, pf2=pf2, f4=f4, quiet=True)
+ f5 = tmpname('tmpA13_') # interior
+ grass.mapcalc("$f5 = if($pf==1,5,0)", pf=pf, pff=pff, f5=f5, quiet=True)
+ f6 = tmpname('tmpA14_') # undetermined
+ grass.mapcalc("$f6 = if($pf>0.6 && $pf<1 && $pf2==0,6,0)", pf=pf, pf2=pf2, pff=pff, f6=f6, quiet=True) # undetermined
+ Index = tmpname('tmpA15_')
grass.run_command("r.series", input=[f1,f2,f3,f4,f5,f6], output=Index,
method="sum", quiet=True)
- indexfin2 = tmpname('indexfin2_')
+ indexfin2 = tmpname('tmpA16_')
grass.mapcalc("$if2 = $ipl * $Index", if2=indexfin2, ipl=ipl,
Index=Index, quiet=True)
#create categories
- indexfin3 = tmpname('indexfin3_')
+ indexfin3 = tmpname('tmpA17_')
fd4, tmprul = tempfile.mkstemp()
text_file = open(tmprul, "w")
text_file.write("0 = 0 nonforest\n")
@@ -305,6 +302,7 @@
text_file.write("4 = 4 perforated\n")
text_file.write("5 = 5 interior\n")
text_file.write("6 = 6 undetermined\n")
+ text_file.write("* = NULL\n")
text_file.close()
grass.run_command("r.reclass", quiet=True, input=indexfin2, output=indexfin3,
title="fragmentation index", rules=tmprul)
@@ -313,7 +311,7 @@
# Shrink the region
if flag_a:
- regionoriginal = tmpname('region_')
+ regionoriginal = tmpname('tmpA18_')
grass.run_command("g.region", save=regionoriginal, quiet=True, overwrite=True)
reginfo = grass.parse_command("g.region", flags="gp")
NSCOR = SWn * float(reginfo['nsres'])
@@ -325,7 +323,6 @@
w=float(reginfo['w'])+EWCOR,
quiet=True)
grass.mapcalc("$opl = $if3", opl=opl, if3=indexfin3, quiet=True)
- grass.run_command("r.null", map=opl, null=0, quiet=True)
#create color codes
fd5, tmpcol = tempfile.mkstemp()
@@ -343,7 +340,6 @@
os.remove(tmpcol)
# Function call
-
if flag_r:rflag="\n\t-r"
else: rflag=""
if flag_t: tflag="\n\t-t"
More information about the grass-commit
mailing list