Old forest fragmentation scripts

forest_fraction_riitters.sh

#!/bin/sh -x
############################################################################
#
# MODULE:       r.forestfrag
#
# AUTHOR(S):    Emmanuel Sambale [hidden email] [hidden email]
#
# 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_ALLc
export LC_NUMERIC=C

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

# get map (assuming 'landclass96 is a raster with forest/non-forest categories)
r.mapcalc 1990veg2=$GIS_OPT_input

# set null values (water areas will not be included)
#r.null map=1990veg2 setnull=6

#recode data: nonforest = 0, forest = 1
r.mapcalc "A=(if(1990veg2==1,1,0))"

echo "computing pf values ..."
# create grid with all values = 1 (needed to calculate pf):
r.mapcalc "B=if(1990veg2,1,1)"

# C: number of forest pixels within a 3x3 moving window
# D: create grid with value=number of pixels within 3x3 moving window

r.neighbors input=A output=C method=sum size=3 --o
r.neighbors input=B output=D method=sum size=3 --o

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

## generate grid where value=1 if forest pixel is surrounded by forest pixels in cardinal directions:
# 0--x--0
# |  |  |
# x--x--x
# |  |  |
# 0--x--0

# 3x3 moving window
r.mapcalc "F4 = (A * A [1,0]) * (A * A [-1,0]) * (A * A [0,-1]) * (A * A [0,1])"

# create grid with pixel-value = number of pixels completely surrounded by forest within 3x3 moving window:
r.neighbors input=F4 output=F5 method=sum size=3 --o

# 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 = Pff

r.mapcalc "pf2 = pf - pff"

r.mapcalc "f1=if(pf<0.4,1,0)" # patch
r.mapcalc "f2=if(pf>=0.4 && pf<0.6,2,0)" # transitional
r.mapcalc "f3=if(pf>=0.6 && pf2<0,3,0)" # edge
r.mapcalc "f4=if(pf>0.6 && pf2>0,4,0)" # perforate
r.mapcalc "f5=if(pf==1 && pff==1,5,0)" # interior
r.mapcalc "f6=if(pf>0.6 && pf<1 && pf2==0,6,0)" # undetermined

r.mapcalc "index=f1+f2+f3+f4+f5+f6"

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

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
EOF

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

r.mapcalc indexfin4=indexfin3

echo "Temporary files deleted ...."
  g.remove rast=A,B,C,D,E,F,f1,f2,f3,f4,f5,f6,pf2
  g.remove rast=indexfin3
  g.remove rast=indexfin2
  g.remove rast=index

#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
EOF

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

## create reclass table "recl.txt":
#cat > recl2.txt
#0 = 0 nonforest
#1 = 1 patch
#2 = 2 transitional
#3 = 3 edge
#4 = 4 perforated
#5 = 5 interior
# 6 thru 8 = 6 undef

r.forestfrag

#!/bin/sh -x

############################################################################
#
# MODULE:       r.forestfrag
#
# AUTHOR(S):    Emmanuel Sambale esambale@yahoo.com emmanuel.sambale@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.
#
#############################################################################

# raster of landcover values of at least 3 values (1=forest, 2=non-forest land, 3=water)
COVER=forest_bin
WINDOW=5
LAND=landmask
############## nothing to change below

if  [ -z "$GISBASE" ] ; then
    echo "You must be in GRASS GIS to run this program."
    exit 1
fi
 
#unset LC_ALL
#export LC_NUMERIC=C


#Redone script
#set to current input map region
g.region rast=$COVER -pm

#recode landcover into binary forest/nonforest
#recode data where nonforest = 0 & forest = 1
r.mapcalc "A=(if('$COVER' == 1,1,0))"

# count number of non-forest pixels value=1
r.mapcalc "B=(if('$LAND' == 1,1,null()))"


echo "computing pf values ..."
 
# C number of forest pixels in a 5x5 window
# D number of land pixels in a 5x5 window

r.neighbors input=A output=C method=sum size="$WINDOW" --o
r.neighbors input=B output=D method=sum size="$WINDOW" --o


r.mapcalc "pf = (C*1.0/D)"

# !!!the above command doesn't give percentage values when used at full zoom/region

echo "computing pff values ..."
## pixels with both forest from the cardinal directions
# 0--x--0
# |  |  |
# x--x--x
# |  |  |
# 0--x--0
 
r.mapcalc "E = (A * A [1,0]) * (A * A [-1,0]) * (A * A [0,-1]) * (A * A [0,1])"
 
#count both forest pixels
r.neighbors input=E output=F method=sum size="$WINDOW" --o

r.mapcalc "pff = (F*1.0/C)"

#echo "computing fragmentation index ..."


r.mapcalc "pff2 = pf - pff"

#index 
r.mapcalc "f1 = if (pf > 0.9,1,null())"
r.mapcalc "f2 = if (pf > 0.6 && pf < 0.9 && pff2 > 0,2,null())" 
r.mapcalc "f3 = if (pf > 0.6 && pf < 0.9 && pff2 < 0.1,3,null())"
r.mapcalc "f4 = if (pf > 0.4 && pf < 0.6,4,null())"
r.mapcalc "f5 = if (pf > 0.0 && pf < 0.39,5,null())"
r.mapcalc "f6 = if (pf > 0.6 && pff2 == 0.00,6,null())"

# mop up all unclassifief forest to undertermined forest
r.mapcalc "f7 =(if('$COVER' == 1,6,null()))"

r.mask in=A maskcats=1 -o
r.patch in=f1,f2,f3,f4,f5,f6,f7 out=index --o
g.remove rast=MASK

#create categories
cat recl.txt
cat recl.txt | r.reclass index out=index_fin title="frag index" --o

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


d.mon x0
r.colors map=pf color=grey
d.rast.leg pf

d.mon x1
r.colors map=pff color=grey
d.rast.leg pff

d.mon x2
d.rast.leg index_fin

d.mon x3
d.rast.leg forest_bin

r.report index_fin units=h,c,p

recl.txt

0 = 0 nonforest
1 = 1 interior
2 = 2 perforated
3 = 3 edge
4 = 4 transitional
5 = 5 patch
6 thru 8 = 6 undetermined