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

[Glynn Clements]

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 at virgin.net>