Opened 12 years ago

Last modified 3 years ago

#857 new enhancement

add simple neighborhood functions to r.mapcalc

Reported by: dickeya Owned by: grass-dev@…
Priority: normal Milestone: 7.6.2
Component: Raster Version: svn-trunk
Keywords: r.mapcalc Cc:
CPU: All Platform: All

Description

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)

Change History (10)

comment:1 in reply to:  description Changed 12 years ago by 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.
  1. 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?.

comment:2 Changed 9 years ago by neteler

Milestone: 6.4.07.0.0
Version: unspecifiedsvn-trunk

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

comment:3 Changed 6 years ago by martinl

Milestone: 7.0.07.0.5

comment:4 Changed 5 years ago by martinl

Milestone: 7.0.57.3.0

comment:5 Changed 5 years ago by martinl

Milestone: 7.3.07.4.0

Milestone renamed

comment:6 Changed 4 years ago by neteler

Milestone: 7.4.07.4.1

Ticket retargeted after milestone closed

comment:7 Changed 4 years ago by neteler

Milestone: 7.4.17.4.2

comment:8 Changed 3 years ago by martinl

Milestone: 7.4.27.6.0

All enhancement tickets should be assigned to 7.6 milestone.

comment:9 Changed 3 years ago by martinl

Milestone: 7.6.07.6.1

Ticket retargeted after milestone closed

comment:10 Changed 3 years ago by martinl

Milestone: 7.6.17.6.2

Ticket retargeted after milestone closed

Note: See TracTickets for help on using tickets.