[GRASSLIST:4893] Re: conditional growing of a region

[Roger Bivand]

> Dear all,
>
> I have been trying to regionalise a DEM into particular landforms. In
> this instance, I am after identifying terraces. My original attempt was
> with a classification type method whereby I identified threshold in
> variables such as slope gradients and flatness. While this is one  way
> of dealing with this question, I am a little disappointed by the
> results because they don't respect the visual geomorphological
> interpretation of the landscape. It is striking that terraces are chunks
>  of land with a particular aspect and slope value (shallow slope
> obviously), yet thresholding the slope and aspect ignore the
> topological adjacency of pixels. Therefore, I'd like to build an
> algorithm that would work like the "magic wand" in photoshop. That  is a
> method that looks at the characteristics of a seed point and test  the
> membership of adjacent pixels according to some decision rule.
>
> In practical terms, how does one go about comparing each pixel of  a
> buffer zone to a seed point. I found the r.grow module which will  do to
> expand the seed point by one pixel in all direction, but once  this ring
> of neighbours is determined, how does one test them one  by one? Or is
> there more nifty trick to deal with this?

You can get quite a long way with the indices used by r.mapcalc:

layer[-1,-1] is offset left and down (check the manual - directions are
not always intuitive!!) by one cell from the [0,0] cell.

You could also look at r.mfilter - by designing a filter matrix, you can
"package" the comparison style you would like without having to repeat
r.mapcalc scripts. I would try with r.mapcalc first to get a feel for
things, then convert your best filter to an r.mfilter matrix. For
instance, Chrisman's textbook has a filter example on p. 166 (1st edition)
of a slope filter:

TITLE 3x3 Sobel x filter
MATRIX 3
-1 0 1
-2 0 2
-1 0 1
DIVISOR 8
TYPE P

TITLE 3x3 Sobel y filter
MATRIX 3
-1 -2 -1
0 0 0
1 2 1
DIVISOR 8
TYPE P

with pseudo-code to filter:

slope.in.x = Filter elevation
      using "Sobel x filter";
slope.in.y = Filter elevation
      using "Sobel y filter";
slope = sqrt((slope.in.x ^ 2)
      + (slope.in.y ^ 2));

An r.mfilter matrix for taking the average of N, S. E, and W cells - one
step - is:

TITLE 3x3 rook's move
MATRIX 3
0 1 0
1 0 1
0 1 0
DIVISOR 0
TYPE P

r.mfilter gives lots of flexibility, matrices can be bigger, and you can
fine-tune the filtering to pick out the kinds of things you need -
provided they don't change too much across the region (aspect might if
valleys curve).

Roger

>
> Thanks for your help,
>
> Thomas
>
> Thomas Dewez
> PhD Student
> Dept Geography & Earth Sciences
> Brunel University (West London)
> Uxbridge UB8 3PH, Mddx
> United Kingdom
>
> Phone: +44-(0)1895-203215
> Fax:     +44-(0)1895-203217
>
> e-mail: thomas.dewez at brunel.ac.uk


-- 
Roger Bivand
NHH, Breiviksveien 40, N-5045 Bergen, Norway
(travelling but still accessible)