add simple neighborhood functions to r.mapcalc

[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. It's especially hard to do an average when nodata cells are present.

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)

[glynn]

Replying to 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.

[neteler]

Perhaps the new Python interface in GRASS 7 is a useful start here.

[martinl]

Milestone renamed

[neteler]

Ticket retargeted after milestone closed

[martinl]

All enhancement tickets should be assigned to 7.6 milestone.

[martinl]

Ticket retargeted after milestone closed

[martinl]

Ticket retargeted after milestone closed