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

Tue Apr 11 03:35:50 EDT 2006

Hallo,

On Mon, Apr 10, 2006 at 06:10:33PM +0800, maning sambale wrote:
> 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.
> 

If you are using mostly r.mapcalc, it should not crashe - in which step
does this happen?

some operations should be performed faster, if you would use g.copy and
r.reclass instead of r.mapcalc -- look after "!!!" string below

> 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

Jachym

> 
#!/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 -> g.copy
g.copy  rast=$GIS_OPT_input,1990veg2

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

#recode data where nonforest = 0 & forest = 1
# !!! r.mapcalc -> r.reclass
#r.mapcalc "A=(if(1990veg2 == 4,1))+(if(1990veg2 == 1||1990veg2 ==
2||1990veg2 == 5,0))"
r.reclass in=1990veg2 out=A << EOF
4 = 1
2 = 5
2 = 5
end
EOF

echo "computing pf values ..."
# count number of non-forest pixels value=1
# r.mapcalc -> r.reclass !!!
#r.mapcalc "B=if(1990veg2 == 1||1990veg2 == 2||1990veg2 == 4||1990veg2 == 5,1)"
r.reclass in=1990veg2 out=B << EOF
1 = 1
2 = 1
4 = 1
5 = 1
end
EOF

# 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 -> r.reclass
# 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.reclass in=pff4 out=F11 << EOF
1 = 1
end
EOF

r.reclass in=pff4 out=F21 << EOF
2 = 1
end
EOF

r.mapcalc "F31 = (pff4 == 3) && (pff3 == 1)"
r.mapcalc "F41 = (pff4 == 3) && (pff3 == 2)"

r.reclass in=pf out=F51 << EOF
1 = 1
end
EOF

r.mapcalc "F61 = (pff4 == 3) && (pff3 == 3)"

# !!! would this work too?
# 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 "index = F11 + F21*2 + F31*3 + F41*4 + F51*5 + F61*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 -> g.copy !!!
g.copy indexfin3,indexfin4

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

-- 
Jachym Cepicky                                  
e-mail: jachym.cepicky at centrum.cz           
URL: http://les-ejk.cz                      
GPG: http://les-ejk.cz/gnupg_public_key/    
-----------------------------------------   
OFFICE:                                     
Department of Geoinformation Technologies   
LDF MZLU v Brně                             
Zemědělská 3                                
613 00 Brno                                 
e-mail: xcepicky at node.mendelu.cz            
URL:    http://mapserver.mendelu.cz
Tel.:   +420 545 134 514