[GRASSLIST:418] Re: Fwd: Fwd: Developing a forest fragmentation index
Jachym Cepicky
jachym.cepicky at centrum.cz
Tue Mar 28 08:49:31 EST 2006
Hallo,
I have only few comments to your script (look after "!!!" string in the
script):
#!/bin/sh
# Create pf map; Pf = (# forest pixels) / (# non-forest pixels)
# A=forest pixels
# B=non-forest pixels
# category vegation code
#1 Agriculture
#2 Brushland
#3 Clouds
#4 Forest
#5 Grassland
#6 River
# !!!
#%Module
#% description: r.forest
#%End
#%option
#% key: input
#% type: string
#% description: raster input map
#% required : yes
#%end
#%option
#% key: output
#% type: string
#% description: raster output map
#% required : yes
#%end
if [ "$1" != "@ARGS_PARSED@" ] ; then
exec $GISBASE/etc/bin/cmd/g.parser "$0" "$@"
fi
# !!! get map
g.copy rast=$GIS_OPT_input,$GIS_OPT_output
# !!! set null values
r.null map=$GIS_OPT_output setnull=3,6
#recode data
# !!! count number of forest pixels value=1
r.mapcalc "A=(if(1990veg == 4,1))+(if(1990veg == 1||1990veg ==
2||1990veg == 5,0))"
# count number of non-forest pixels value=1
r.mapcalc "B=if(1990veg == 1||1990veg == 2||1990veg == 4||1990veg == 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=5
r.neighbors input=B output=D method=sum size=5
## pixels with both forest from the cardinal directions
#create filters
# !!! could this not be done with help of r.mapcalc neighbord
# operations?
r.mfilter input=A output=Nf filter=filtrN.txt
r.mfilter input=A output=Sf filter=filtrS.txt
r.mfilter input=A output=Wf filter=filtrW.txt
r.mfilter input=A output=Ef filter=filtrE.txt
# !!! recode filtrd to 2=1 1=0 and recode
r.mapcalc "F3 = eval ( \
Nf2 = if(Nf == 1||Nf == 0,0))+(if(Nf == 2,1)), \
Sf2=(if(Sf == 1||Sf == 0,0))+(if(Sf == 2,1)), \
Ef2=(if(Ef == 1||Ef == 0,0))+(if(Ef == 2,1)), \
Wf2=(if(Wf == 1||Wf == 0,0))+(if(Wf == 2,1)), \
Nf2 + Sf2 + Ef2 + Wf2))"
# !!! count both forest pixels
r.neighbors input=F4 output=F5 method=sum size=5
#frag model
#code Category Pf Pff
# 1 Patch < 40%
# 2 Transitional >= 40% but < 60%
# 3 Edge > 60% Pf - Pff > 0 1
# 4 Perforated > 60% Pf ??? Pff < 0 2
# 5 Interior 1
# 6 Undetermined > 0.6 Pf = Pff 3
r.mapcalc "Pff3 = eval ( \
Pff2 = (1.0 * C) /(1.0 * D) - (1.0 * F5)/E, \
(if (pff2 > 0,1)) + (if (pff2 < 0,2)) + (if (pff2 == 0,3)) \
)"
#0-.4 1
#.4 - .6 2
#.6 - 0.99 3
#1 4
# !!!
r.mapcalc "Pff4 = (if ((F6/E) < 0.4,1)) + (if ((F6/E) >= 0.4 && (F6/E) < .6,2))
+ (if ((F6/E) >= 0.6 && (F6/E) < .99,3)) / (if ((F6/E) == 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 "indexfin = if (index == 0 ||index == 7 ||index == 8,0) + if
(index == 1,1) + if /
(index == 2,2) + if (index == 3,3) + if (index == 4,4) + if (index ==
5,5) + if (index == 6,6)"
g.remove rast=A,B,C,D,E,Ef,Ef2,F,F3,F4,F5,F6,Nf,Nf2,Sf,Wf,Wf2,F11,F21,F31,F41,F51,F61,Pff2,Pff3,Pff4,Sf2,index
I made only few changes, so you could have an idea, how to rewrite your
script a bit, so it does not perform so many calculations - or so it
does not to create so many temporary maps. I did not modify the end of
the script.
Hamish allready wrote you, that the order of priorities is
1. It works
2. Clarity
--
3. Efficiency
Your script worked and now you can try to make it work better. When you
finish these optimalisations, try to answer following question: Is this
script fast enought for my purposes? If our answer is "no", then you
have to rewrite it (or at least some parts of it) to C. Before you do
so, you should get really working prototype.
>
> Cheers,
> Maning
>
> On 3/22/06, maning sambale <emmanuel.sambale at gmail.com> wrote:
> > Thank you for the responses. Will dig into to the manuals you recommended.
> > Will post later when I made some progress. Right now i'm into the
> > generating the Pf which is simply r.neighbors and r.mapcalc commands .
> >
> > What is bugging me is generating Pff, which counts from the moving
> > window the forest pixel pairs in 4 cardinal directions. (see
> > attached)
> >
> > Pff = (# of pixel pairs with both pixel forest) / (# of pixel pairs
> > with at least one pixel forest)
> >
> > Cheers,
> > Maning
> >
> > PS. how do I install r.li tarball in grass in cygwin?
> >
> > On 3/21/06, Jachym Cepicky <jachym.cepicky at centrum.cz> wrote:
> > > Hallo,
> > > I have some experience with neighborhood operations on rastermaps, so I
> > > could help you to make the module together.
> > >
> > > Perhaps, you could be able make some prototype with r.mapcalc module?
> > > (see http://grass.itc.it/gdp/raster/mapcalc-algebra.pdf and
> > > http://grass.itc.it/gdp/raster/mapcalc.pdf)
> > >
> > > Jachym
> > >
> > > maning sambale wrote:
> > >
> > > >Dear list,
> > > >
> > > >I hope somebody can help me on this. I am a MS student studying
> > > >forest loss and fragmentation here in the Philippines. I am interested
> > > >in implementing Mr. Kurt Ritters forest fragmentation model
> > > >(http://www.ecologyandsociety.org/vol4/iss2/art3/) and James Hurd's
> > > >adaptation to watersehd units
> > > >(http://resac.uconn.edu/publications/tech_papers/pdf_paper/Forest_Fragmentation_ASPRS2002.pdf)
> > > >in my study area and relate the fragmentation to socio-economic
> > > >drivers.
> > > >
> > > >Mr. Ritters gave me the C source code for his model, however, I am not
> > > >programmer so I don't know how to implement these. Right now, I am
> > > >doing my image analysis (i.* modules) and classification using GRASS
> > > >and import the output to Arcview were EPA has an Attila extension
> > > >(http://www.epa.gov/esd/land-sci/attila/index.htm) for classifying
> > > >them into the forest fragmenation index. What I want to do is create
> > > >a GRASS script for the forest fragmentation index.
> > > >
> > > >I have been using GRASS for just a couple of months and I am not very
> > > >familiar with most of the commands. I am willing to collaborate with
> > > >anybody from this list doing related research on forest fragmentation.
> > > > As I have mentioned above, I am not a programmer thus my limitation
> > > >in understanding C code. However, I am willing to do GRASS related
> > > >contributions related to my study.
> > > >
> > > >The method is described below:
> > > >
> > > >The basis for the forest fragmentation index was the fragmentation
> > > >model developed by Ritters et al. (2000) which was used and expanded
> > > >by Hurd et. al (2002) on Landsat TM data. The objective of the model
> > > >is to provide additional information beyond the commonly used measures
> > > >of forest proportions (i. e. percent of forest cover) by providing
> > > >information on the spatial configuration of fragmentation and
> > > >connectivity of forest within the study
> > > >region.
> > > >
> > > >The amount of forest and its occurrence is measured as adjacent forest
> > > >pixels within a 5 x 5 windows surrounding each forest pixel. This
> > > >information will be used to classify the window by the type of
> > > >fragmentation. The result was 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
> > > >(Riitters et al. 2000). The model generates two values that
> > > >characterize a forest pixel located at the center of a sliding window
> > > >of fixed size, Pf and Pff:
> > > >Pf is the proportion of pixels in the window that are forested. Pff
> > > >(strictly) is 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.
> > > >
> > > >Pf = ( # of forest pixels) / (# of all non water pixels)
> > > >Pff = (# of pixel pairs with both pixel forest) / (# of pixel pairs
> > > >with at least one pixel forest)
> > > >
> > > >The classification model identifies six fragmentation categories: (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
> > > >
> > > >Any tips in implementing the generation of forest fragmentation index
> > > >either as a GRASS script or a step-by-step GRASS commands would be
> > > >very helpful (r.mapcalc syntax?).
> > > >
> > > >I have also attached Ritters's C code as well as my previous
> > > >correspondence to him. Also attached are images from Hurd's study.
> > > >
> > > >Thank you very and looking forward for any response.
> > > >
> > > >Cheers,
> > > >
> > > >Maning Sambale
> > > >
> > > >
> > > >---------- Forwarded message ----------
> > > >From: Kurt H Riitters <kriitters at fs.fed.us>
> > > >Date: Oct 23, 2004 2:46 AM
> > > >Subject: Re: Developing a forest fragmentation index
> > > >To: emmanuel sambale <emmanuel.sambale at gmail.com>
> > > >Cc: Kurt H Riitters <kriitters at fs.fed.us>
> > > >
> > > >
> > > >Emanuel,
> > > > Attached is the main C program. As I must say to everyone, I am not able
> > > >to provide technical support for using the program, and I do not guarantee
> > > >that it is bug-free. You can share this with anyone you want to, but if
> > > >you share it, then you must answer any questions about it.
> > > >
> > > > Probably the interesting part will be the algorithm for the moving
> > > >window, and the definition of the functions (like Pf and Pff). If it were
> > > >me, I'd think about adapting those parts to new programs of your design.
> > > >
> > > > The input / output format is assumed to be an old format that I adapted
> > > >from somewhere else...and it is not standard and is not used by anyone
> > > >else... so you will want to re-write the input / output routines to read
> > > >files in other formats.
> > > >
> > > > I do not have any written documentation, but there are a lot of comments
> > > >in the code that will help understand the program flow.
> > > >
> > > > Lastly, since this code is a research tool, there are many places
> > > >containing references to code or subroutines that are no longer included.
> > > >
> > > > Oh, I should say... this program uses a moving window to calculate the
> > > >X and Y values that are in the model that appears in that 2000 paper. One
> > > >run of the program makes a map of Pf; a second run makes a map of Pff.
> > > >Once you have the output maps of the Pf and Pff values, it is necessary to
> > > >put them together to perform the classification into categories like
> > > >'perforated' 'edge' etc. Arc or envi would be convenient for that.
> > > >
> > > > Best regards,
> > > >
> > > >Kurt Riitters
> > > >RWU-4803
> > > >Forest Health Monitoring
> > > >US Forest Service
> > > >3041 Cornwallis Road
> > > >Research Triangle Park NC 27709
> > > >919-549-4015
> > > >kriitters at fs.fed.us
> > > >
> > > >
> > > >
> > > >>>Emmanuel,
> > > >>> Thank you for the message.
> > > >>> Of course you may use the methods; they are public knowledge.
> > > >>>
> > > >>> I wrote C programs to process the landcover maps. I am happy to send
> > > >>>them to you, however, if you are not a programmer then it will be very
> > > >>>difficult to adapt the programs to your problems. Even another
> > > >>>
> > > >>>
> > > >programmer
> > > >
> > > >
> > > >>>may find it difficult since I am not a very good programmer either.
> > > >>>
> > > >>> I do not have any arc or envi scripts since I use the C programs...
> > > >>>
> > > >>>
> > > >But
> > > >
> > > >
> > > >>>it is a very simple moving window (arc: focal function) approach that
> > > >>>should be easy to implement in another language... For example, arc
> > > >>>
> > > >>>
> > > >should
> > > >
> > > >
> > > >>>be able to do the "Pf" part of the model with the Grid: focalsum
> > > >>>
> > >
> > >
> >
> >
> > --
> > |---------|----------------------------------------------------------|
> > | __.-._ |"Ohhh. Great warrior. Wars not make one great." -Yoda |
> > | '-._"7' |"Freedom is still the most radical idea of all" -N.Branden|
> > | /'.-c |Linux registered user #402901, http://counter.li.org/ |
> > | | /T |Philippine Biodiversity Web Map |
> > | _)_/LI |www.geocities.com/esambale/philbiodivmap/philbirds.html |
> > |---------|----------------------------------------------------------|
> >
> >
>
>
> --
> |---------|----------------------------------------------------------|
> | __.-._ |"Ohhh. Great warrior. Wars not make one great." -Yoda |
> | '-._"7' |"Freedom is still the most radical idea of all" -N.Branden|
> | /'.-c |Linux registered user #402901, http://counter.li.org/ |
> | | /T |Philippine Biodiversity Web Map |
> | _)_/LI |www.geocities.com/esambale/philbiodivmap/philbirds.html |
> |---------|----------------------------------------------------------|
>
>
> --
> |---------|----------------------------------------------------------|
> | __.-._ |"Ohhh. Great warrior. Wars not make one great." -Yoda |
> | '-._"7' |"Freedom is still the most radical idea of all" -N.Branden|
> | /'.-c |Linux registered user #402901, http://counter.li.org/ |
> | | /T |Philippine Biodiversity Web Map |
> | _)_/LI |www.geocities.com/esambale/philbiodivmap/philbirds.html |
> |---------|----------------------------------------------------------|
>
--
Jachym Cepicky
e-mail: jachym.cepicky at centrum.cz
URL: http://les-ejk.cz
More information about the grass-user
mailing list