[GRASS-SVN] r64306 - grass-addons/grass7/raster/r.forestfrag
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Jan 24 15:02:23 PST 2015
Author: pvanbosgeo
Date: 2015-01-24 15:02:23 -0800 (Sat, 24 Jan 2015)
New Revision: 64306
Added:
grass-addons/grass7/raster/r.forestfrag/Makefile
grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py
Removed:
grass-addons/grass7/raster/r.forestfrag/r.forestfrag
Modified:
grass-addons/grass7/raster/r.forestfrag/r.forestfrag.html
Log:
Rewrite of the plugin in Python, bugfixes and now name and description are written to output layer.
Added: grass-addons/grass7/raster/r.forestfrag/Makefile
===================================================================
--- grass-addons/grass7/raster/r.forestfrag/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.forestfrag/Makefile 2015-01-24 23:02:23 UTC (rev 64306)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.forestfrag
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.forestfrag/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Deleted: grass-addons/grass7/raster/r.forestfrag/r.forestfrag
===================================================================
--- grass-addons/grass7/raster/r.forestfrag/r.forestfrag 2015-01-24 20:41:54 UTC (rev 64305)
+++ grass-addons/grass7/raster/r.forestfrag/r.forestfrag 2015-01-24 23:02:23 UTC (rev 64306)
@@ -1,377 +0,0 @@
-#!/bin/sh
-############################################################################
-#
-# MODULE: r.forestfrag
-#
-# AUTHOR(S): Emmanuel Sambale, Stefan Sylla, Paulo van Breugel
-#
-# PURPOSE: Creates forest fragmentation index map from a forest-/
-# non-forest raster; The index map is based on Riitters,
-# K., J. Wickham, R. O'Neill, B. Jones, and E. Smith. 2000.
-# in: Global-scale patterns of forest fragmentation.
-# Conservation Ecology 4(2): 3. [online]
-# URL:http://www.consecol.org/vol4/iss2/art3/
-#
-# NOTES This addon is an adaptation of the r.forestfrag.sh addon, to
-# work on GRASS 7.0 and to add the option to have the
-# size of the moving window determined by the user. For now
-# it does not work on GRASS 6.4. Use r.forestfrag.sh instead.
-#
-# COPYRIGHT: (C) 1997-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: creates forest fragmentation index from a GRASS raster map (where forest=1, non-forest=0) based on a method developed by Riitters et. al (2000). The index is computed using a moving window of user-defined size (default=3).
-#% keyword: raster
-#% keyword: forest
-#% keyword: fragmentation index
-#% keyword: Riitters
-#%End
-
-#%Option
-#% key: input
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Name of forest raster map (where forest=1, non-forest=0)
-#% required : yes
-#%End
-
-#%option
-#% key: output
-#% type: string
-#% gisprompt: new
-#% description: Name output layer
-#% key_desc: name
-#% required: yes
-#%end
-
-#%option
-#% key: window
-#% type: integer
-#% description: Window size (odd number)
-#% key_desc: number
-#% answer : 3
-#% required: yes
-#%end
-
-#%flag
-#% key: r
-#% description: Set region to raster?
-#%end
-
-#%flag
-#% key: t
-#% description: keep temporary maps
-#%END
-
-#%flag
-#% key: s
-#% description: Run r.report on output map
-#%END
-
-#%flag
-#% key: a
-#% description: shrink the output map?
-#%END
-
-
-# For testing purposes
-#-----
-#GIS_OPT_INPUT=ForestMosaic at Testing
-#GIS_OPT_OUTPUT=FMtest
-#GIS_OPT_WINDOW=3
-#GIS_FLAG_R=0
-#GIS_FLAG_T=1
-#GIS_FLAG_S=1
-#GIS_FLAG_A=1
-#-----
-
-#=======================================================================
-## GRASS team recommandations
-#=======================================================================
-
-if [ -z "$GISBASE" ] ; then
- echo "You must be in GRASS GIS to run this program." 1>&2
- exit 1
-fi
-
-if [ "$1" != "@ARGS_PARSED@" ] ; then
- exec g.parser "$0" "$@"
-fi
-
-eval `g.findfile el=cell file=$GIS_OPT_INPUT`
-if [ ! "$file" ] ; then
- echo "Raster map '$GIS_OPT_input' not found in mapset search path"
- exit 1
-fi
-
-#=======================================================================
-## Prepare data
-#=======================================================================
-
-#set to current input map region (user option, default=current region)
-#--------------------------------------------------------------------
-if [ $GIS_FLAG_R -eq 1 ];
-then
- echo "setting region to input map ..."
- g.region raster=$GIS_OPT_INPUT
-fi
-
-# Set root name of temporary output files
-#--------------------------------------------------------------------
-tmpl=forestfrag$$_
-
-# get map (assuming input-map is a cell-raster with forest=1,
-# and setting all other values to 0 (assuming to be non-forest=)
-#--------------------------------------------------------------------
-g.copy raster=${GIS_OPT_INPUT},${tmpl}
-# recode forest-map: nonforest = 0, forest = 1
-r.reclass input=${tmpl} output=${tmpl}A rules=- <<EOF
-1 = 1
-* = 0
-end
-EOF
-
-#=======================================================================
-# computing pf values
-#=======================================================================
-
-echo "computing pf values ..."
-
-## 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:
-r.neighbors input=${tmpl}A output=${tmpl}A2 method=sum size=$GIS_OPT_WINDOW
-
-# generate grid with pixel-value=number of pixels in moving window:
-r.mapcalc "${tmpl}C3 = 1.0*$GIS_OPT_WINDOW^2" --overwrite
-
-# create pf map
-r.mapcalc file=- << EOF
-${tmpl}A3 = 1.0 * ${tmpl}A2
-${tmpl}pf = (${tmpl}A3/${tmpl}C3)
-EOF
-
-#=======================================================================
-# computing pff values
-#=======================================================================
-
-echo "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"
-
-# calculate number of 'forest-forest' pairs
-#=======================================================================
-
-echo "... calculate number of 'forest-forest' pairs"
-
-SW=$GIS_OPT_WINDOW
-tmpf=`g.tempfile pid=$$`
-echo -n "${tmpl}F1 = " >> $tmpf
-
-# Compute matrix dimensions
-SWn=$(echo "($SW - 1) / 2" | bc)
-
-# Write mapcalc expression to tmpf - rows
-rsub=$(seq $SWn -1 -$SWn)
-csub=$(seq -$SWn 1 $(echo "($SWn-1)" | bc))
-
-for k in $rsub
-do
- for l in $csub
- do
- echo -n "${tmpl}A[$k,$l]*${tmpl}A[$k,$(echo "$l+1" | bc)] + " >> $tmpf
- done
-done
-
-# Write mapcalc expression to tmpf - columns
-rsub=$(seq -$SWn 1 $SWn)
-csub=$(seq $SWn -1 $(echo "(-$SWn+1)" | bc))
-
-for k in $rsub
-do
- for l in $csub
- do
- echo -n "${tmpl}A[$l,$k]*${tmpl}A[$(echo "$l-1" | bc),$k] + " >> $tmpf
- done
-done
-
-echo -n "0" >> $tmpf
-r.mapcalc file=$tmpf
-unlink $tmpf
-
-# number of 'forest-forest' pairs
-#=======================================================================
-
-tmpf2=`g.tempfile pid=$$$$`
-echo -n "${tmpl}F2 = " >> $tmpf2
-
-# Compute matrix dimensions
-SWn=$(echo "($SW - 1) / 2" | bc)
-
-# Write mapcalc expression to tmpf - rows
-rsub=$(seq $SWn -1 -$SWn)
-csub=$(seq -$SWn 1 $(echo "($SWn-1)" | bc))
-
-for k in $rsub
-do
- for l in $csub
- do
- echo -n "if((${tmpl}A[$k,$l]+${tmpl}A[$k,$(echo "$l+1" | bc)])>0,1) + " >> $tmpf2
- done
-done
-
-# Write mapcalc expression to tmpf - columns
-rsub=$(seq -$SWn 1 $SWn)
-csub=$(seq $SWn -1 $(echo "(-$SWn+1)" | bc))
-
-for k in $rsub
-do
- for l in $csub
- do
- echo -n "if((${tmpl}A[$l,$k]+${tmpl}A[$(echo "$l-1" | bc),$k])>0,1) + " >> $tmpf2
- done
-done
-
-echo -n "0" >> $tmpf2
-r.mapcalc file=$tmpf2
-unlink $tmpf2
-
-# create pff map
-r.mapcalc --overwrite file=- << EOF
-${tmpl}F1 = 1.0 * ${tmpl}F1
-${tmpl}F2 = 1.0 * ${tmpl}F2
-${tmpl}pff = ${tmpl}F1/${tmpl}F2
-EOF
-
-#=======================================================================
-# computing fragmentation index
-#=======================================================================
-
-echo "computing fragmentation index ..."
-# (1) patch pf < 0.4
-# (2) transitional 0.4 < pf < 0.6
-# (3) edge pf > 0.6 and pf - pff < 0
-# (4) perforated pf > 0.6 and pf - pff > 0
-# (5) interior for which pf = 1.0
-# (6) undetermined pf > 0.6 and pf = pff
-
-r.mapcalc file=- << EOF
-${tmpl}pf2 = ${tmpl}pf - ${tmpl}pff
-${tmpl}f1 = if(${tmpl}pf<0.4,1,0)
-${tmpl}f2 = if(${tmpl}pf>=0.4 && ${tmpl}pf<0.6,2,0)
-${tmpl}f3 = if(${tmpl}pf>=0.6 && ${tmpl}pf2<0,3,0)
-${tmpl}f4 = if(${tmpl}pf>0.6 && ${tmpl}pf2>0,4,0)
-${tmpl}f5 = if(${tmpl}pf==1 && ${tmpl}pff==1,5,0)
-${tmpl}f6 = if(${tmpl}pf>0.6 && ${tmpl}pf<1 && ${tmpl}pf2==0,6,0)
-${tmpl}index = ${tmpl}f1+${tmpl}f2+${tmpl}f3+${tmpl}f4+${tmpl}f5+${tmpl}f6
-${tmpl}indexfin2 = ( ${tmpl}A * ${tmpl}index )
-EOF
-
-#create categories
-echo "creating colors and categories ... "
-r.reclass input=${tmpl}indexfin2 output=${tmpl}indexfin3 title="frag index" rules=- <<EOF
-0 = 0 nonforest
-1 = 1 patch
-2 = 2 transitional
-3 = 3 edge
-4 = 4 perforated
-5 = 5 interior
-6 = 6 undef
-EOF
-
-
-# Shrink the region
-g.region save=rforestfrag987654321 --overwrite
-if [ $GIS_FLAG_A -eq 1 ];
-then
- NSRES=`g.region -p -g | grep nsres= | cut -f2 -d"="`
- EWRES=`g.region -p -g | grep ewres= | cut -f2 -d"="`
- NSCOR=$(echo "$SW*$NSRES*0.5" | bc)
- EWCOR=$(echo "$SW*$EWRES*0.5" | bc)
- g.region n=n-$NSCOR s=s+$NSCOR e=e-$EWCOR w=w+$EWCOR
-fi
-
-r.mapcalc "${GIS_OPT_OUTPUT} = ${tmpl}indexfin3"
-r.null ${GIS_OPT_OUTPUT} null=0
-
-#create color codes
-r.colors ${GIS_OPT_OUTPUT} --quiet rules=- << EOF
-0 255:255:0
-1 215:48:39
-2 252:141:89
-3 254:224:139
-4 217:239:139
-5 26:152:80
-6 145:207:96
-EOF
-
-#=======================================================================
-# Clean up (user option: copy pf, pff, pf2 first)
-#=======================================================================
-
-if [ $GIS_FLAG_T -eq 1 ];
-then
- g.copy raster=${tmpl}pf,${GIS_OPT_OUTPUT}_pf
- g.copy raster=${tmpl}pff,${GIS_OPT_OUTPUT}_pff
- g.copy raster=${tmpl}pf2,${GIS_OPT_OUTPUT}_pf2
-fi
-
-echo "Deleting temporary files ...."
-g.remove -f -b --quiet type=raster pattern=${tmpl}*
-
-#=======================================================================
-# computing fragmentation index
-#=======================================================================
-
-if [ $GIS_FLAG_S -eq 1 ];
-then
-
- g.copy raster=${GIS_OPT_OUTPUT},${tmpl}_REPORT
- # Shrink the region if output map was not trimmed
- if [ $GIS_FLAG_A -eq 0 ];
- then
- NSRES=`g.region -p -g | grep nsres= | cut -f2 -d"="`
- EWRES=`g.region -p -g | grep ewres= | cut -f2 -d"="`
- NSCOR=$(echo "$SW*$NSRES*0.5" | bc)
- EWCOR=$(echo "$SW*$EWRES*0.5" | bc)
- g.region save=rforestfrag987654321
- g.region n=n-$NSCOR s=s+$NSCOR e=e-$EWCOR w=w+$EWCOR
- fi
-
- echo "generate map reports ..."
- r.report map=${tmpl}_REPORT units=h,p -n
- echo " "
- echo "------------------------------------------------------"
- echo "Please note in order to avoid edge effects, the region"
- echo "is reduced / shrinked with a number of raster cells"
- echo "equal to 1/2 the moving window size before running"
- echo "r.report, unless the user has already selected the "
- echo "option to trim the output map with number of cells equal"
- echo "to 1/2 the moving window size (default)."
- echo "------------------------------------------------------"
- echo " "
-fi
-
-g.remove -f type=raster name=${tmpl}_REPORT --quiet
-g.region region=rforestfrag987654321
-g.remove --quiet -f type=region name=rforestfrag987654321
-
-
-
-
-
-
-
-
-
-
Modified: grass-addons/grass7/raster/r.forestfrag/r.forestfrag.html
===================================================================
--- grass-addons/grass7/raster/r.forestfrag/r.forestfrag.html 2015-01-24 20:41:54 UTC (rev 64305)
+++ grass-addons/grass7/raster/r.forestfrag/r.forestfrag.html 2015-01-24 23:02:23 UTC (rev 64306)
@@ -1,59 +1,65 @@
<h2>DESCRIPTION</h2>
-<em>r.forestfrag</em> computes the forest fragmentation following the methodology
-proposed by Riitters et. al (2000). See http://www.consecol.org/vol4/iss2/art3/ for
-a detailed explanation.
+<em>r.forestfrag</em> Computes the forest fragmentation following
+the methodology proposed by Riitters et. al (2000). See
+<a href="http://www.consecol.org/vol4/iss2/art3/">this article</a>for a
+detailed explanation.
-<p>It follows a “sliding window” algorithm with overlapping windows.
-The amount of forest and its occurence as adjacent forest pixels within fixed-
-area "moving-windows" surrounding each forest pixel is measured. The
-window size is user-defined. The result is stored at the location of the center pixel.
-Thus, a pixel value in the derived map refers to "between-pixel" fragmentation
-around the corresponding forest location.
+<p>It follows a "sliding window" algorithm with overlapping windows.
+The amount of forest and its occurence as adjacent forest pixels
+within fixed- area "moving-windows" surrounding each forest pixel is
+measured. The window size is user-defined. The result is stored at
+the location of the center pixel. Thus, a pixel value in the derived
+map refers to "between-pixel" fragmentation around the corresponding
+forest location.
-<p> Let Pf be the proportion of pixels in the window that are forested. Define
-Pff (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. Pff (roughly) estimates the conditional probability that, given a
-pixel of forest, its neighbor is also forest. [...] the classification
-model [...] identifies six fragmentation categories:
+<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:
-<ul>
-<li> (1) interior, for which Pf = 1.0</li>
-<li> (2) patch, Pf < 0.4</li>
-<li> (3) transitional, 0.4 < Pf < 0.6</li>
-<li> (4) edge, Pf > 0.6 and Pf - Pff < 0</li>
-<li> (5) perforated, Pf > 0.6 and Pf - Pff > 0</li>
-<li> (6) undetermined, Pf > 0.6 and Pf = Pff</li>
-</ul>
+<ol start="0">
+<li>No forest</li>
+<li>interior, for which Pf = 1.0</li>
+<li>patch, Pf < 0.4</li>
+<li>transitional, 0.4 < Pf < 0.6</li>
+<li>edge, Pf > 0.6 and Pf - Pff < 0</li>
+<li>perforated, Pf > 0.6 and Pf – Pff > 0</li>
+<li>undetermined, Pf > 0.6 and Pf = Pff</li>
+</ol>
-<h2>NOTES</h2>
-<p>The moving window size is user-defined (default=3), but must be given as a odd number.
-if the user gives an even number, the number is silently reduced to the nearest
-odd number < user-defined number.
-<p>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.
-<p>If the user selects to run r.report on the output map, but has selected to leave
-the output untrimmed, r.report will be run on a copy of the output map that is trimmed.
+<h2>NOTES</h2> <p>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.
+<p>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.
+
+<p>The function respects the region. The user has however the option
+to set the region to match the input layer.
+
<h2>SEE ALSO</h2>
-This addon is based on the r.forestfrag.sh
-"http://grasswiki.osgeo.org/wiki/AddOns/GRASS_6#r.forestfrag" addon, adapted to make
-it compatible with GRASS 7.0. It introduces the option to change the moving window
-size and unlike the original addon, it respects the region (there is still the option
-to set the region to the input raster). This addon will not work on GRASS 6.4, please
-use the above mentioned addon.
+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).
-<h2>AUTHORS</h2>
-Emmanuel Sambale: Author original script<br>
-Stefan Sylla: Revised / updated <br>
-Paulo van Breugel: Revised / updated for GRASS 7.0, added the option to determine moving window size.<br>
+<h2>AUTHORS</h2>
+<p><em>Emmanuel Sambale</em>: Author original script
+<p><em>Stefan Sylla</em>: Revised / updated the r.forestfrag.sh for
+GRASS GIS 6
+<p><em>Paulo van Breugel</em>: Rewritten in Python
+for use with GRASS GIS 7.0, and added extra options.
Contact: paulo at ecodiv.org
-<h2>REFERENCES</h2>
-Riitters, K., J. Wickham, R. O'Neill, B. Jones, and E. Smith. 2000. Global-scale patterns of forest
-fragmentation. Conservation Ecology 4(2): 3. [online] URL: http://www.consecol.org/vol4/iss2/art3/
+<h2>REFERENCES</h2>
+<p>Riitters, K., J. Wickham, R. O'Neill, B. Jones,
+and E. Smith. 2000. Global-scale patterns of forest fragmentation.
+Conservation Ecology 4(2): 3. [online] URL:
+http://www.consecol.org/vol4/iss2/art3/
-<p>
-<i>Last changed: $Date$</i>
+<p><i>Last changed: $Date$</i>
Added: grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py
===================================================================
--- grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py (rev 0)
+++ grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py 2015-01-24 23:02:23 UTC (rev 64306)
@@ -0,0 +1,447 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+############################################################################
+#
+# MODULE: r.forestfrag
+#
+# AUTHOR(S): Emmanuel Sambale, Stefan Sylla and Paulo van Breugel
+# Maintainer: Paulo van Breugel (paulo at ecodiv.org)
+#
+# PURPOSE: Creates forest fragmentation index map from a forest-non-forest
+# raster; The index map is based on Riitters, K., J. Wickham,
+# R. O'Neill, B. Jones, and E. Smith. 2000. in: Global-scale
+# patterns of forest fragmentation. Conservation Ecology 4(2): 3.
+# [online] URL:http://www.consecol.org/vol4/iss2/art3/
+#
+#
+# COPYRIGHT: (C) 1997-2015 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: Computes the forest fragmentation index (Riitters et. al 2000).
+#% keyword: raster
+#% keyword: forest
+#% keyword: fragmentation index
+#% keyword: Riitters
+#%End
+
+#%Option
+#% key: input
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of forest raster map (where forest=1, non-forest=0)
+#% required : yes
+#%End
+
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Name output layer
+#% key_desc: name
+#% required: yes
+#%end
+
+#%option
+#% key: window
+#% type: integer
+#% description: Moving window size (odd number)
+#% key_desc: number
+#% answer : 3
+#% required: yes
+#%end
+
+#%flag
+#% key: r
+#% description: Set region to raster?
+#%end
+
+#%flag
+#% key: t
+#% description: keep pf and pff maps
+#%END
+
+#%flag
+#% key: s
+#% description: Run r.report on output map
+#%END
+
+#%flag
+#% key: a
+#% description: trim the output map to avoid border effects?
+#%END
+
+
+# For testing purposes
+#options = {"input":"forest96 at forestfrag", "output":"mosaicindex", "window":"5"}
+#flags = {"r":False, "t":True, "s":True, "a":False}
+
+#----------------------------------------------------------------------------
+# Standard
+#----------------------------------------------------------------------------
+
+# import libraries
+import os
+import sys
+import numpy as np
+import uuid
+import atexit
+import tempfile
+import string
+import grass.script as grass
+
+if not os.environ.has_key("GISBASE"):
+ grass.message("You must be in GRASS GIS to run this program.")
+ sys.exit(1)
+
+# create set to store names of temporary maps to be deleted upon exit
+clean_rast = set()
+
+def cleanup():
+ for rast in clean_rast:
+ grass.run_command("g.remove", flags="f", type="rast", name=rast, quiet=True)
+
+#----------------------------------------------------------------------------
+# Functions
+#----------------------------------------------------------------------------
+
+def CheckLayer(envlay):
+ ffile = grass.find_file(envlay, element = 'cell')
+ if ffile['fullname'] == '':
+ grass.fatal("The layer " + envlay + " does not exist.")
+
+def tmpname(prefix):
+ tmpf = prefix + str(uuid.uuid4())
+ tmpf = string.replace(tmpf, '-', '_')
+ clean_rast.add(tmpf)
+ return tmpf
+
+#----------------------------------------------------------------------------
+# Main
+#----------------------------------------------------------------------------
+
+def main():
+
+ #------------------------------------------------------------------------
+ # Variables
+ #------------------------------------------------------------------------
+
+ ipl = options['input']
+ CheckLayer(ipl)
+ opl = options['output']
+ wz = int(options['window'])
+ if wz % 2 == 0:
+ grass.fatal("Please provide an odd number for the moving window)")
+ flag_r = flags['r']
+ flag_t = flags['t']
+ flag_s = flags['s']
+ flag_a = flags['a']
+
+
+ #set to current input map region (user option, default=current region)
+ if flag_r:
+ grass.message("setting region to input map ...")
+ grass.run_command('g.region', quiet=True, raster=ipl)
+
+ # Check if map values are limited to 1 and 0
+ tmpdir = tempfile.mkdtemp()
+ tmpfile = tmpdir + 'forestfrag1'
+ grass.run_command("r.stats", flags="cn", overwrite=True,
+ quiet=True, input=ipl, output=tmpfile, separator="|")
+ tstf = np.loadtxt(tmpfile, delimiter="|", dtype="int")
+ os.remove(tmpfile)
+ if min(tstf[:,0]) != 0 or max(tstf[:,0]) != 1:
+ grass.fatal("Your input map must be binary, with values 0 and 1")
+
+ #------------------------------------------------------------------------
+ # computing pf values
+ #------------------------------------------------------------------------
+ grass.info("step 1: computing pf values ...\n")
+
+ # 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_')
+ grass.run_command("r.neighbors", quiet=True, input=ipl,
+ output=tmpA2, method="sum", size=wz)
+
+ # 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)
+
+ # Clean up
+ grass.run_command("g.remove", quiet=True,
+ flags="f", type="raster",
+ name=[tmpA2, tmpC3])
+
+ #------------------------------------------------------------------------
+ # 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"
+
+ grass.info("step 2: computing pff values ...\n")
+
+ # calculate number of 'forest-forest' pairs
+ #------------------------------------------------------------------------
+
+ # Compute matrix dimensions
+ SWn= int((wz - 1) / 2)
+
+ # Write mapcalc expression to tmpf - rows
+ tmpfile2 = tempfile.mkstemp()[1]
+ tmpl4 = tmpname('pff1_')
+ text_file = open(tmpfile2, "w")
+ text_file.write(tmpl4 + " = ")
+
+ # Write mapcalc expression to tmpf (rows and columns)
+ rsub=range(SWn, -SWn-1, -1)
+ 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) + "] + ")
+ 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
+ + "[" + str(l-1) + "," + str(k) + "] + ")
+ text_file.write("0")
+ text_file.close()
+
+ grass.run_command("r.mapcalc", file=tmpfile2, quiet=True)
+ os.remove(tmpfile2)
+
+ # number of 'forest patches
+ #------------------------------------------------------------------------
+
+ tmpfile3 = tempfile.mkstemp()[1]
+ tmpl5 = tmpname('pff2_')
+ text_file = open(tmpfile3, "w")
+ text_file.write(tmpl5 + " = ")
+
+
+ # Write mapcalc expression to tmpf - rows and columns
+ 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(k) + "," + str(l) + "]+"
+ + ipl + "[" + 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("0")
+ text_file.close()
+
+ grass.run_command("r.mapcalc", file=tmpfile3, quiet=True)
+ os.remove(tmpfile3)
+
+ # create pff map
+ #------------------------------------------------------------------------
+
+ pff = tmpname('pff3_')
+ grass.mapcalc("$pff = float($tmpl4) / float($tmpl5)", tmpl4=tmpl4,
+ tmpl5=tmpl5, pff=pff,quiet=True)
+
+ #------------------------------------------------------------------------
+ # computing fragmentation index
+ #------------------------------------------------------------------------
+
+ grass.info("step 3: computing fragmentation index ...\n")
+ pf2 = tmpname('pf2_')
+ 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_')
+ grass.run_command("r.series", input=[f1,f2,f3,f4,f5,f6], output=Index,
+ method="sum", quiet=True)
+ indexfin2 = tmpname('indexfin2_')
+ grass.mapcalc("$if2 = $ipl * $Index", if2=indexfin2, ipl=ipl,
+ Index=Index, quiet=True)
+
+ #create categories
+ indexfin3 = tmpname('indexfin3_')
+ tmprul = tempfile.mkstemp()
+ text_file = open(tmprul[1], "w")
+ text_file.write("0 = 0 nonforest\n")
+ text_file.write("1 = 1 patch\n")
+ text_file.write("2 = 2 transitional\n")
+ text_file.write("3 = 3 edge\n")
+ text_file.write("4 = 4 perforated\n")
+ text_file.write("5 = 5 interior\n")
+ text_file.write("6 = 6 undetermined\n")
+ text_file.close()
+ grass.run_command("r.reclass", quiet=True, input=indexfin2, output=indexfin3,
+ title="fragmentation index", rules=tmprul[1])
+ os.remove(tmprul[1])
+
+ # Shrink the region
+ if flag_a:
+ regionoriginal = tmpname('region_')
+ grass.run_command("g.region", save=regionoriginal, quiet=True, overwrite=True)
+ reginfo = grass.parse_command("g.region", flags="gp")
+ NSCOR = SWn * float(reginfo['nsres'])
+ EWCOR = SWn * float(reginfo['ewres'])
+ grass.run_command("g.region",
+ n=float(reginfo['n'])-NSCOR,
+ s=float(reginfo['s'])+NSCOR,
+ e=float(reginfo['e'])-EWCOR,
+ 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)
+ grass.run_command("g.remove", flags="f", quiet=True, type="rast", name=[indexfin3])
+ grass.run_command("g.remove", flags="f", quiet=True, type="rast", name=[indexfin2])
+
+ #create color codes
+ tmpcol = tempfile.mkstemp()
+ text_file = open(tmpcol[1], "w")
+ text_file.write("0 255:255:0\n")
+ text_file.write("1 215:48:39\n")
+ text_file.write("2 252:141:89\n")
+ text_file.write("3 254:224:139\n")
+ text_file.write("4 217:239:139\n")
+ text_file.write("5 26:152:80\n")
+ text_file.write("6 145:207:96\n")
+ text_file.close()
+ grass.run_command("r.colors", quiet=True, map=opl, rules=tmpcol[1])
+ os.remove(tmpcol[1])
+
+ # Function call
+
+ if flag_r:rflag="\n\t-r"
+ else: rflag=""
+ if flag_t: tflag="\n\t-t"
+ else: tflag=""
+ if flag_s: sflag="\n\t-s"
+ else: sflag=""
+ if flag_a: aflag="\n\t-a"
+ else: aflag=""
+ desctxt = "r.forestfrag \n\tinput=" + ipl + "\n\toutput=" + opl + \
+ "\n\twindow=" + str(wz) + rflag + tflag + sflag + aflag
+
+ # Write metadata for main layer
+ tmphist = tempfile.mkstemp()
+ text_file = open(tmphist[1], "w")
+ text_file.write("Forest fragmentation index (6 classes) following Riiters et al. (2000)\n\
+\t(1) patch\n\t(2) transitional\n\t(3) edge\n\t(4) perforated\n\t(5) interior\n\t(6) undetermined\n")
+ text_file.write("\ncreated using:\n")
+ text_file.write(desctxt)
+ text_file.close()
+ grass.run_command("r.support", map=opl, title="Forest fragmentation",
+ units="Fragmentation classes (6)",
+ source1="based on " + ipl,
+ source2="",
+ description="Forst fragmentation index (6 classes)",
+ loadhistory=tmphist[1])
+
+ # Write metadata for intermediate layers
+ if flag_t:
+ # pf layer
+ tmphist = tempfile.mkstemp()
+ text_file = open(tmphist[1], "w")
+ text_file.write("created using:\n")
+ text_file.write(desctxt + "\n")
+ text_file.close()
+ grass.run_command("r.support", map=pf,
+ title="Proportion forested",
+ units="Proportion",
+ source1="based on " + ipl,
+ source2="",
+ description="Proportion of pixels in the moving window that is forested",
+ loadhistory=tmphist[1])
+
+ # pff layer
+ tmphist = tempfile.mkstemp()
+ text_file = open(tmphist[1], "w")
+ text_file.write("created using:\n")
+ text_file.write(desctxt + "\n\n")
+ text_file.write("Proportion of all adjacent (cardinal directions only) pixel pairs that\n")
+ text_file.write("include at least one forest pixel for which both pixels are forested.\n")
+ text_file.write("It thus (roughly) estimates the conditional probability that, given a\n")
+ text_file.write("pixel of forest, its neighbor is also forest")
+ text_file.close()
+ grass.run_command("r.support", map="Mosaic1_pf at forestfrag",
+ title="Conditional probability neighboring cell is forest",
+ units="Proportion",
+ source1="based on " + ipl,
+ source2="",
+ description="Probability neighbor of forest cell is forest",
+ loadhistory=tmphist[1])
+
+
+ #------------------------------------------------------------------------
+ # Report fragmentation index and names of layers created
+ #------------------------------------------------------------------------
+
+ if flag_s:
+ grass.run_command("r.report", map=opl, units=["h","p"],
+ flags="n", page_width=50, quiet=True)
+ grass.info("\n")
+ grass.info("The following layers were created\n")
+ grass.info("The fragmentation index: " + opl +"\n")
+ if flag_t:
+ grass.run_command("g.rename", quiet=True, raster=[pf,opl + "_pf"])
+ grass.run_command("g.rename", quiet=True, raster=[pff,opl + "_pff"])
+ grass.info("The proportion forested (pf): " + opl + "_pf\n")
+ grass.info("The proportion forested pixel pairs: " + opl + "_pff\n")
+ else:
+ grass.run_command("g.remove", flags="f", quiet=True, type="rast", name=[pf,pff])
+
+ #------------------------------------------------------------------------
+ # Clean up
+ #------------------------------------------------------------------------
+ if flag_a:
+ grass.run_command("g.region", region=regionoriginal, quiet=True, overwrite=True)
+ grass.run_command("g.remove", flags="f", type="region", quiet=True, name=regionoriginal)
+ os.removedirs(tmpdir)
+ os.remove(tmphist[1])
+ grass.run_command("g.remove", flags="f", quiet=True, type="rast", name=[f1,f2,f3,f4,f5,f6, pf2])
+ grass.run_command("g.remove", flags="f", quiet=True, type="rast", name=[tmpl4,tmpl5,Index])
+
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ atexit.register(cleanup)
+ sys.exit(main())
+
+
+
+
Property changes on: grass-addons/grass7/raster/r.forestfrag/r.forestfrag.py
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list