[GRASS-dev] Re: [GRASS GIS] #857: add simple neighborhood functions to r.mapcalc

[GRASS GIS]

#857: add simple neighborhood functions to r.mapcalc
--------------------------+-------------------------------------------------
  Reporter:  dickeya      |       Owner:  grass-dev at lists.osgeo.org
      Type:  enhancement  |      Status:  new                      
  Priority:  normal       |   Milestone:  6.4.0                    
 Component:  Raster       |     Version:  unspecified              
Resolution:               |    Keywords:  r.mapcalc                
  Platform:  All          |         Cpu:  All                      
--------------------------+-------------------------------------------------
Comment (by glynn):

 Replying to [ticket:857 dickeya]:

 > I find myself using a lot of neighbor functions in r.mapcalc, but find
 the syntax to be tedious.  Anything larger than a 3x3 neighborhood and
 writing out "raster[-1,-1] + raster[-1,0] + raster[-1,-1] + ..." is quite
 a task.

 For a sum, use `r.mfilter[.fp]` or `r.neighbors ... weight=...`

 > It's especially hard to do an average when nodata cells are present.

 Convert nulls first in a separate pass.

 > It would be great to have some built in neighborhood functions that
 would handle the math, taking into account nodata cells, and working on a
 rectangular or circular area.  Something like:
 >
 {{{
 neighborhood(raster, width, height, shape, function)
 }}}

 There area a couple of issues with this:

 1. The functions which make up the bulk of r.mapcalc all operate on 1-D
 arrays of size G_window_cols(). So we can't "leverage" the existing
 functions for this purpose.

 2. The above neighborhood() function cannot be implemented using the
 existing function interface; it would need to be a new language construct,
 recognised by the parser, the evaluator, and the I/O code.

 If neighborhood() was a normal function, the above call would result in
 its first argument being the current row of the specified raster, which is
 no good for a neighborhood operation (and all of its other arguments would
 also be row-sized arrays of CELL/FCELL/DCELL values; a numeric literal is
 passed as a row-sized array filled with that value).

 I suspect that it might be more fruitful to add new aggregates to
 r.neighbors.

 That still leaves a gap between what can be implemented with a combination
 of r.mapcalc + r.neighbors + r.mapcalc and what could be implemented with
 a more general language. E.g. there isn't any combination which will allow
 you to calculate:
 {{{
 output[r,c] = sum<i,j>(input[r+i,c+j] <op> kernel[i,j])
 }}}
 for any <op> except multiplication.

 Implementing that specific case in r.mapcalc (or r.neighbors) would be a
 lot of work, but not entirely out of the question. My main concern is that
 it would be the thin end of the wedge; either we draw a line in the sand,
 or we end up with Perl.

 IMHO, if you can't do what you want with r.mapcalc + r.neighbors, you'll
 need a system with enough (virtual) memory to use Matlab/Octave/NumPy.

-- 
Ticket URL: <http://trac.osgeo.org/grass/ticket/857#comment:1>
GRASS GIS <http://grass.osgeo.org>