[GRASS5] Re: [GRASSLIST:1549] Use of neighbourhood operator in r.mapcalc

Roy Sanderson wrote:

I have a raster map with over 200 points, each with a different category
value, and I need to buffer around each point by one pixel, but retaining
the original category value in the buffered cells. Since r.grow and
r.buffer assign new categories to the buffered cells, I've tried to use the
neighbourhood operators in r.mapcalc. Using a script, containing the line

r.mapcalc buffmap = if(oldmap, oldmap, \
if(oldmap[0,-1], oldmap[0,-1] \
  if(oldmap[0,1], oldmap[0,1] \
   if(oldmap[-1,0], oldmap[-1,0] \
    if(oldmap[1,0], oldmap[1,0])))))

which is very similar to the example script in Shapiro & Westervelt's
r.mapcalc guide. The puzzle is that whilst this works fine in Grass 4.3,
in Grass5.0pre5 the buffering is incomplete, with the western buffer pixel
on cell too far to the west. Have I overlooked something?

Prior to 5.0.2, the horizontal modifier for CELL (integer) maps was
effectively doubled (the shift was applied twice). This bug was fixed
in 5.0.2. Note that the bug was introduced when r.mapcalc was
completely re-written, somewhere around the time of the -pre3 release,
so this isn't applicable to 4.3.

A number of significant bugs have been fixed (in r.mapcalc and in
other programs) between -pre5 and the current 5.0.2 release. I
strongly recommend upgrading to either 5.0.2 or 5.0.3-RC5.

Aside: this operation seems to be sufficiently popular that we should
consider adding an option to either r.buffer or r.grow to support it.
Using r.mapcalc quickly becomes impractical if you want a buffer
larger than one cell.

--
Glynn Clements <glynn.clements@virgin.net>

> r.mapcalc buffmap = if(oldmap, oldmap, \
> if(oldmap[0,-1], oldmap[0,-1] \
> if(oldmap[0,1], oldmap[0,1] \
> if(oldmap[-1,0], oldmap[-1,0] \
> if(oldmap[1,0], oldmap[1,0])))))

Aside: this operation seems to be sufficiently popular that we should
consider adding an option to either r.buffer or r.grow to support it.
Using r.mapcalc quickly becomes impractical if you want a buffer
larger than one cell.

[enter r.grow2]

Great, now I can get rid of my crappy script.

My use:
I've used this with a series of "region codes" at various site locations
which 'grow' outwards in a loop until the entire non-null region is
filled.

General method:
s.to.rast taking category number for raster value.
r.mapcalc base map so area of interest is 0 and no-go is NULL.
r.patch sites,basemap
r.grow_endless_loop
^C, check, r.grow_endless_loop [or {while find(==0)} might work..]

This allows me to 'grow' around corners; s.surf.rst+r.reclass has no
concept of impenetrable barriers, eg data from two parallel lakes
pollute each other.

One moderately important caveat:
For correct results I found I had to alternate the sign of the [-1,0]'s
(search order) after each iteration in order for it to make it around
certain corners and/or fill in some pixels which only have data on one
adjoining side. I cycled the search order but maybe randomly is more
appropriate.

thanks, this will be useful.
Hamish

Hamish wrote:

One moderately important caveat:
For correct results I found I had to alternate the sign of the [-1,0]'s
(search order) after each iteration in order for it to make it around
certain corners and/or fill in some pixels which only have data on one
adjoining side. I cycled the search order but maybe randomly is more
appropriate.

I don't understand what you are talking about here. The only
difference that the ordering will make is that it determines *which*
non-NULL neighbour fills in a given NULL cell.

Do you have some non-NULL cells which are being plotted in the same
colour as NULLs? That could result in confusion.

--
Glynn Clements <glynn.clements@virgin.net>

> One moderately important caveat:
> For correct results I found I had to alternate the sign of the
> [-1,0]'s(search order) after each iteration in order for it to make
> it around certain corners and/or fill in some pixels which only have
> data on one adjoining side. I cycled the search order but maybe
> randomly is more appropriate.

I don't understand what you are talking about here. The only
difference that the ordering will make is that it determines *which*
non-NULL neighbour fills in a given NULL cell.

Ok, this was a problem with the r.mapcalc solution:

r.mapcalc buffmap = if(oldmap, oldmap, \
if(oldmap[0,-1], oldmap[0,-1] \
  if(oldmap[0,1], oldmap[0,1] \
   if(oldmap[-1,0], oldmap[-1,0] \
    if(oldmap[1,0], oldmap[1,0])))))

If you had:
(* = NULL, 0 is cell to be filled, 1 is a value)

111
*0*
***

The 0 would be set to NULL (if(oldmap[0,-1]); I reset the 0 base-map at
each iteration to remove these NULLs; and continuing the loop, the cell
would never get a value. Only changing the search order so oldmap[-1,0]
came first would fill the cell correctly.

That was one problem, the other was when two categories met, the
0-valued cell would be consistently filled by one of the categories over
the other based on their position relative to the target cell, due to
the search order always testing for say the left hand cell first.

I think the attached PNG is showing a combination of this and the NULL
problem mentioned above, but it has been 6 months and I think I fixed
the worst of it by hand, so it isn't the most informative example,
sorry.
[purple started growing from off screen SW, green from off screen NE,
they meet in the middle]

Summary: I think a constant search order leads to undesirable artifacts.

I haven't studied r.grow2 to comment on how it would act.

Do you have some non-NULL cells which are being plotted in the same
colour as NULLs? That could result in confusion.

No.

regards,
Hamish

(attachments)

grow_biased.png

Hamish wrote:

> > One moderately important caveat:
> > For correct results I found I had to alternate the sign of the
> > [-1,0]'s(search order) after each iteration in order for it to make
> > it around certain corners and/or fill in some pixels which only have
> > data on one adjoining side. I cycled the search order but maybe
> > randomly is more appropriate.
>
> I don't understand what you are talking about here. The only
> difference that the ordering will make is that it determines *which*
> non-NULL neighbour fills in a given NULL cell.

Ok, this was a problem with the r.mapcalc solution:

> r.mapcalc buffmap = if(oldmap, oldmap, \
> if(oldmap[0,-1], oldmap[0,-1] \
> if(oldmap[0,1], oldmap[0,1] \
> if(oldmap[-1,0], oldmap[-1,0] \
> if(oldmap[1,0], oldmap[1,0])))))

If you had:
(* = NULL, 0 is cell to be filled, 1 is a value)

111
*0*
***

The 0 would be set to NULL (if(oldmap[0,-1]); I reset the 0 base-map at
each iteration to remove these NULLs; and continuing the loop, the cell
would never get a value. Only changing the search order so oldmap[-1,0]
came first would fill the cell correctly.

The r.mapcalc example above (once you replace the missing commas) is
clearly designed for the case where "empty" cells are represented by
zero rather than NULL. If you want it to work reliably with NULLs for
empty cells, you would have to use not(isnull(...)), i.e.

  r.mapcalc buffmap = \
   if(not(isnull(oldmap)), oldmap, \
    if(not(isnull(oldmap[0,-1])), oldmap[0,-1], \
     if(not(isnull(oldmap[0, 1])), oldmap[0, 1], \
      if(not(isnull(oldmap[-1,0])), oldmap[-1,0], \
       if(not(isnull(oldmap[ 1,0])), oldmap[ 1,0], \
        null())))))

An expression of the form if(NULL,...) always evaluates to NULL
regardless of the other arguments.

That was one problem, the other was when two categories met, the
0-valued cell would be consistently filled by one of the categories over
the other based on their position relative to the target cell, due to
the search order always testing for say the left hand cell first.

I think the attached PNG is showing a combination of this and the NULL
problem mentioned above, but it has been 6 months and I think I fixed
the worst of it by hand, so it isn't the most informative example,
sorry.
[purple started growing from off screen SW, green from off screen NE,
they meet in the middle]

Summary: I think a constant search order leads to undesirable artifacts.

The handling of cells for which there is no unique nearest non-NULL
neighbour is bound to be arbitrary. Ultimately, the issue is how you
wish to assign categories to the "border" cells.

A fixed search order will always use the category from the cells which
lie in the first direction tries (i.e. North in the above example). A
varying search order will tend to "share" the categories; that may be
an advantage (e.g. if you are subsequently computing area totals), or
a disadvantage (e.g. due to semi-random artifacts in the shape of the
boundary).

I haven't studied r.grow2 to comment on how it would act.

The latest version of r.grow2 allows a radius to be specified, in
which case it will behave more like r.buffer. This is bound to be
preferable to an iterative approach, as any "approximations" (i.e.
choosing an arbitrary neighbour when there are multiple candidates)
are only performed once, not once per iteration.

FWIW, I will be adding options to use either the Manhattan or maximum
metrics as well as the Euclidean metric (repeated iterations of r.grow
effectively use the Manhattan metric, expanding single cells to
"diamonds").

--
Glynn Clements <glynn.clements@virgin.net>