[GRASSLIST:644] Re: Fwd: Fwd: Developing a forest fragmentation index

[maning sambale]

Dear grasslist,

I have made some revisions from my forest fragmentation script.  Thank
you for your contributions.  I have incorporated some of your
suggestions (those that I understood) as well as adding options for
specifying window size, removal of temporary maps and displaying
summary class statistics.

One of the problems I encounter is the processing of very large maps. 
My P4 CPU clock shoots to 100% and then crashes.  I was thinking maybe
of cutting the images to manageable chunks and then process them by
batches.

Next is to overlay watershed basins and create summary statistics for
each class for each basin and import the class statistics to R
readable format.  For further statistical analysis.

Thanks for your ideas!

Maning

r.forest.frag script

#!/bin/sh
############################################################################
#
# MODULE:       r.forestfrag
#
# AUTHOR(S):    Emmanuel Sambale esambale at yahoo.com emmanuel.sambale at gmail.com
#
# PURPOSE:      Creates forest fragmentation index map from a raster
vegetation file
#                The index map was based on 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/
#
# COPYRIGHT:    (C) 2005 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
#%End
#%flag
#%  key: r
#%  description: remove temporary maps
#%END
#%option
#% key: input
#% type: string
#% gisprompt: old,cell,raster
#% description: Name of vegetation map
#% required : yes
#%End
#%option
#% key: window
#% type: integer
#% description: window size default is 5
#% answer : 5
#% required : yes
#%END

if  [ -z "$GISBASE" ] ; then
	echo "You must be in GRASS GIS to run this program."
	exit 1
fi

if [ "$1" != "@ARGS_PARSED@" ] ; then
   exec $GISBASE/bin/g.parser "$0" "$@"
fi

unset LC_ALL
export LC_NUMERIC=C

#set to current input map region
g.region rast=$GIS_OPT_input -p

# get map
r.mapcalc 1990veg2=$GIS_OPT_input

# set null values
r.null map=1990veg2 setnull=3,6

#recode data where nonforest = 0 & forest = 1
r.mapcalc "A=(if(1990veg2 == 4,1))+(if(1990veg2 == 1||1990veg2 ==
2||1990veg2 == 5,0))"

echo "computing pf values ..."
# count number of non-forest pixels value=1
r.mapcalc "B=if(1990veg2 == 1||1990veg2 == 2||1990veg2 == 4||1990veg2 == 5,1)"

# C number of forest pixels in a 5x5 window
# D number of non-forest pixels in a 5x5 window

r.neighbors input=A output=C method=sum size="$GIS_OPT_window"
r.neighbors input=B output=D method=sum size="$GIS_OPT_window"

# count number of forest pixels value=1

# this one is better create pf map
r.mapcalc << EOF
E = 1.0 * C
F = 1.0 * D
pf = (E/F)
EOF

echo "computing pff values ..."
## pixels with both forest from the cardinal directions
# 0--x--0
# |  |  |
# x--x--x
# |  |  |
# 0--x--0

r.mapcalc "F4 = (A * A [1,0]) * (A * A [-1,0]) * (A * A [0,-1]) * (A * A [0,1])"

#count both forest pixels
r.neighbors input=F4 output=F5 method=sum size="$GIS_OPT_window"

# create pff map
r.mapcalc << EOF
F6 = 1.0 * F5
pff = (F6/E)
EOF

echo "computing fragmentation index ..."
# (1) interior, for which Pf = 1.0; (2) patch, Pf < 0.4; (3)
transitional, 0.4 < Pf < 0.6; (4) edge,
# Pf > 0.6 and Pf - Pff > 0; (5) perforated, Pf > 0.6 and Pf – Pff <
0, and (6) undetermined,
# Pf > 0.6 and Pf = Pffr.mapcalc "Pff2 = pf - pff"

r.mapcalc "Pff3 = (if (pff2 > 0,1)) + (if (pff2 < 0,2)) + (if (pff2 == 0,3))"

# recode pff to 0-.4, .4-.6, .6-.99, 1
r.mapcalc "Pff4 = (if (pff < 0.4,1)) + (if (pff >= 0.4 && pff < .6,2))
+ (if (pff >= 0.6 && pff < .99,3)) + (if (pff == 1,4))"


r.mapcalc "F11 = pff4 == 1"
r.mapcalc "F21 = pff4 == 2"
r.mapcalc "F31 = (pff4 == 3) && (pff3 == 1)"
r.mapcalc "F41 = (pff4 == 3) && (pff3 == 2)"
r.mapcalc "F51 = pf == 1"
r.mapcalc "F61 = (pff4 == 3) && (pff3 == 3)"

r.mapcalc "index = if (if (F11==1,1)) + (if (F21==1,2)) + (if
(F31==1,3)) + (if (F41==1,4)) + (if (F51==1,5)) + (if (F61==1,6))"

r.mapcalc "indexfin2= ( A * index )"

#create color codes

echo "creating color codes and categories ..."
r.colors indexfin2 color=rules << EOF
0 yellow
1 215:48:39
2 252:141:89
3 254:224:139
4 217:239:139
5 26:152:80
6 145:207:96
7 145:207:97
8 145:207:98
EOF

#create categories
cat recl.txt
cat recl.txt | r.reclass indexfin2 out=indexfin3 title="frag index"

r.mapcalc indexfin4=indexfin3

if [ $GIS_FLAG_r -eq 1 ] ; then
  echo "Temporary files deleted ...."
  g.remove rast=A,B,C,D,E,F,F11,F21,F31,F4,F41,F5,F51,F6,F61,Pff2,Pff3,Pff4
  g.remove rast=indexfin3
  g.remove rast=indexfin2
  g.remove rast=index
fi

#create color codes
r.colors indexfin4 color=rules << EOF
0 yellow
1 215:48:39
2 252:141:89
3 254:224:139
4 217:239:139
5 26:152:80
6 145:207:96
7 145:207:97
8 145:207:98
EOF

echo "generate map reports ..."
r.report map=indexfin4 units=h,p null=* nsteps=255 -e -i

reclass table
recl.txt
0 = 0 nonforest
1 = 1 patch
2 = 2 transitional
3 = 3 edge
4 = 4 perforated
5 = 5 interior
6 thru 8 = 6 undet