[GRASS-dev] Re: r.neighbors modification

GRASS GIS GRASS GIS developer
1

[..]

> I am looking for a way to compute r.neighbors for a ring only and not
> for

the

> whole 5x5, 7x7 etc. moving window.
> It is something like a 7x7 windows but minus the values in the 5x5
> windows.
>
> traditional way:
>
> a a a a a
> a b b b a
> a b x b a
> a b b b a
> a a a a a
>
> ring analysis:
>
> a a a a a
> a a
> a x a
> a a
> a a a a a

I'm guessing that should be:
> a a a a a
> a a
> a x a
> a a
> a a a a a
>
>
> I think there is no post-r.neighbors way to compute statistics in such
> a "ring",

It depends upon the aggregate. Some of them could be computed with a
combination of r.neighbors and r.mapcalc commands.

They could all be done using r.mapcalc, but some of them might require
excessively complex expressions (at a minimum, they would all involve
16 map[r,c] terms). r.mapcalc has n-ary min, max, median and mode; sum
and mean are fairly easy to write (you would probably want to process
the nulls separately to simplify the expressions), variance and stddev
aren't hard but would be somewhat more verbose.

I wrote a short neighbourhood ring analysis for average and
variance analysis (currently only up to 11x11 but I will add rings up to
21x21).
Hamish advice how to improve a script (v.cellstats threat) is partly
incorporated but I am grateful for any further suggestions how to improve this
script.

Martin

Of course, the expressions get longer as the window size increases,
and you would need a separate expression for each size (although you
could probably generate it from a script).

> hence I or somebody else has to modify the source code. I am
> willing to learn some more C programming, but definitely need some help
> to get started.

I'm not sure it's worth modifying r.neighbors for this specific case
(it might be worth you modifying your version, but not us shipping a
modified version). More generally useful would be to allow the
neighbourhood set to be read from a file, similar to r.mfilter (or
even just have a list of offsets as an optional argument).

(attachments)

r.ring (13.3 KB)

2

Martin Wegmann wrote:

I wrote a short neighbourhood ring analysis for average and
variance analysis (currently only up to 11x11 but I will add rings up
to 21x21).
Hamish advice how to improve a script (v.cellstats threat) is partly
incorporated but I am grateful for any further suggestions how to
improve this script.

Hi,

just a few more comments-

this only has to happen once:
#%Module
#% description: ring neighborhood analysis
#%END

you might combine all those flags into one option,
#%option
#% key: size
#% type: integer
#% description: compute var & avg of NxN ring
#% options: 3,5,7,9,11
#% answer: 5
#% required : no
#%END

(defaults and not abreviating in descriptions [these get built into help
pages, translated, etc] are nice)

you call a "cleanup" fn but don't define one.

Hamish

3

Hello Hamish,

thanks for your comments, I added most of them but a few questions remain:

On Friday 24 November 2006 07:38, Hamish wrote:

Martin Wegmann wrote:
> I wrote a short neighbourhood ring analysis for average and
> variance analysis (currently only up to 11x11 but I will add rings up
> to 21x21).
> Hamish advice how to improve a script (v.cellstats threat) is partly
> incorporated but I am grateful for any further suggestions how to
> improve this script.

Hi,

just a few more comments-

this only has to happen once:
#%Module
#% description: ring neighborhood analysis
#%END

done

you might combine all those flags into one option,
#%option
#% key: size
#% type: integer
#% description: compute var & avg of NxN ring
#% options: 3,5,7,9,11
#% answer: 5
#% required : no
#%END

very good - implemented

(defaults and not abreviating in descriptions [these get built into help
pages, translated, etc] are nice)

that means that in
#% description:
it is nicer not to use abbreviations, because they are shown in help pages
etc., right?

you call a "cleanup" fn but don't define one.

I think you mean this part:

# what to do in case of user break:
exitprocedure()
{
echo "User break!" 1>&2
cleanup
exit 1
}
# shell check for user break (signal list: trap -l)
trap "exitprocedure" 2 3 15
## copy from r.univar.sh - end

how do I only call a "cleanup"?

thanks, Martin

Hamish

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

4

On Friday 24 November 2006 12:16, Martin Wegmann wrote:
[...]

> you might combine all those flags into one option,
> #%option
> #% key: size
> #% type: integer
> #% description: compute var & avg of NxN ring
> #% options: 3,5,7,9,11
> #% answer: 5
> #% required : no
> #%END

very good - implemented

is it possible to add a flag when all sizes shall be computed?

I tried a bit with || statements in the if-then part, but had no luck so far.

any ideas?

TIA, Martin

[...]

(attachments)

r.ring (65.4 KB)

5

Martin Wegmann wrote:

> > you might combine all those flags into one option,
> > #%option
> > #% key: size
> > #% type: integer
> > #% description: compute var & avg of NxN ring
> > #% options: 3,5,7,9,11
> > #% answer: 5
> > #% required : no
> > #%END
>
> very good - implemented

is it possible to add a flag when all sizes shall be computed?

I tried a bit with || statements in the if-then part, but had no luck
so far.

any ideas?

from the parser standpoint:

#% options: 3,5,7,9,11
#% multiple: yes

then you can have e.g. a 2 cell wide ring with

size=5,7

[There are two options here: multiple individual sizes or combination
sizes. The rest talks about multiple sizes within the same map as that's
more interesting to me.]

[Will r.mfilter do the same as the output _avg map?]

from the g.parser help page:
NOTES
An option can be instructed to allow multiple inputs by adding the
following line:

#% multiple : yes

While this will only directly change the Usage section of the help
screen, the option's environmental string may be easily parsed from
within a script. For example, individual comma separated identities for
an option named "input" can be parsed with the following Bash shell
code:

IFS=,
for opt in $GIS_OPT_INPUT ; do
    ... "$opt"
done

I am not sure how to do the r.mapcalc part then, maybe >> append each
mapcalc ring size rules to a file, then run r.mapcalc redirecting the
conglomorate file to stdin,

r.mapcalc < $TMP.rulesfile

it would take too much to fully paste all combinations in to the script
without some sort of dynamic rule building.

if you want multiple individual maps, set DO_5=0, DO_7=1 in the loop
and then test for that, e.g.:

DO_3=0
DO_5=0
DO_7=0
# ...

IFS=,
for size in $GIS_OPT_SIZE ; do
  if [ $size -eq 3 ] ; then
    DO_3=1
  fi
  if [ $size -eq 5 ] ; then
    DO_5=1
  fi
  if [ $size -eq 7 ] ; then
    DO_7=1
  fi
# ...

done

(or something like that)

Hamish

6

On Tuesday 28 November 2006 00:38, Hamish wrote:

Martin Wegmann wrote:
> > > you might combine all those flags into one option,
> > > #%option
> > > #% key: size
> > > #% type: integer
> > > #% description: compute var & avg of NxN ring
> > > #% options: 3,5,7,9,11
> > > #% answer: 5
> > > #% required : no
> > > #%END
> >
> > very good - implemented
>
> is it possible to add a flag when all sizes shall be computed?
>
> I tried a bit with || statements in the if-then part, but had no luck
> so far.
>
> any ideas?

from the parser standpoint:

#% options: 3,5,7,9,11
#% multiple: yes

then you can have e.g. a 2 cell wide ring with

size=5,7

[There are two options here: multiple individual sizes or combination
sizes. The rest talks about multiple sizes within the same map as that's
more interesting to me.]

[Will r.mfilter do the same as the output _avg map?]

just quick test and reply:

r.covar gives:
1.000000 0.969710
0.969710 1.000000

same algorithm but with moving window shows a few differences (see jpg).
(divide by 1000 to get correlation values between -1 and 1)

Martin

from the g.parser help page:
NOTES
An option can be instructed to allow multiple inputs by adding the
following line:

#% multiple : yes

While this will only directly change the Usage section of the help
screen, the option's environmental string may be easily parsed from
within a script. For example, individual comma separated identities for
an option named "input" can be parsed with the following Bash shell
code:

IFS=,
for opt in $GIS_OPT_INPUT ; do
    ... "$opt"
done

I am not sure how to do the r.mapcalc part then, maybe >> append each
mapcalc ring size rules to a file, then run r.mapcalc redirecting the
conglomorate file to stdin,

r.mapcalc < $TMP.rulesfile

it would take too much to fully paste all combinations in to the script
without some sort of dynamic rule building.

if you want multiple individual maps, set DO_5=0, DO_7=1 in the loop
and then test for that, e.g.:

DO_3=0
DO_5=0
DO_7=0
# ...

IFS=,
for size in $GIS_OPT_SIZE ; do
  if [ $size -eq 3 ] ; then
    DO_3=1
  fi
  if [ $size -eq 5 ] ; then
    DO_5=1
  fi
  if [ $size -eq 7 ] ; then
    DO_7=1
  fi
# ...

done

(or something like that)

Hamish

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

(attachments)

mfilter_ring_correlation.jpg

7

Martin Wegmann wrote:

I wrote a short neighbourhood ring analysis for average and
variance analysis (currently only up to 11x11 but I will add rings up to
21x21).
Hamish advice how to improve a script (v.cellstats threat) is partly
incorporated but I am grateful for any further suggestions how to improve this
script.

Average can be computed using r.mfilter; except, that it needs to be
modified to use floating-point.

Variance could be computed using two passes of r.mfilter and two of
r.mapcalc, although that algorithm is less stable in the case of very
small variance relative to the mean (cf historical r.univar.sh
problems).

Apart from that, I suggest generating the r.mapcalc expression using
awk, e.g.:

  # expects $size to be the radius of the ring
  awk -vsize=$size -vinf=$GIS_OPT_input -voutf=$GIS_OPT_output \
  'BEGIN {
  print outf" = (0\\"
  for (i=-size; i<size; i++) {
  print "+"inf"["(-size)","(-i)"]\\" ;
  print "+"inf"["( size)","( i)"]\\" ;
  print "+"inf"["( i)","(-size)"]\\" ;
  print "+"inf"["(-i)","( size)"]\\" ;
  }
  print ")/(8*size)"
  }' < /dev/null | r.mapcalc

--
Glynn Clements <glynn@gclements.plus.com>