[GRASS-dev] grass7 - python script - moving window: how to simplify/accelerate r.mapcalc?

[Glynn Clements]

Helmut Kudrnovsky wrote:

> in a python script within a mapcalc-expression I have moving window with 6
> cells in all directions and a if(x,a,b) for each cell in
> 
> r.mapcalc "elevation_percentile_step2 = (100.0 / 48.0) * \
> (if(eudem_osttirol[3,3] < eudem_osttirol , 1, 0 ) + if(eudem_osttirol[2,2] < eudem_osttirol , 1, 0 ) \
> + if(eudem_osttirol[1,1] < eudem_osttirol , 1, 0 ) + if(eudem_osttirol[-3,3] < eudem_osttirol , 1, 0 ) \
> + if(eudem_osttirol[-2,2] < eudem_osttirol , 1, 0 ) + if(eudem_osttirol[-1,1] < eudem_osttirol , 1, 0 ) \

[snip]

> + if(eudem_osttirol[-4,-4] < eudem_osttirol , 1, 0 ) + if(eudem_osttirol[-5,-5] < eudem_osttirol , 1, 0 ) \
> + if(eudem_osttirol[-6,-6] < eudem_osttirol , 1, 0 ) + if(eudem_osttirol[-4,0] < eudem_osttirol , 1, 0 ) \
> + if(eudem_osttirol[-5,0] < eudem_osttirol , 1, 0 ) + if(eudem_osttirol[-6,0] < eudem_osttirol , 1, 0 ))"
> 
> any idea/hint to simplify and/or accelerate such a python script with such a
> large moving window?

The most straightforward change is to eliminate the if() functions. 
Any boolean value (e.g. the result of a comparison) is an integer,
with 1 for true, 0 for false. So "if(a<b,1,0)" is equivalent to just
"(a<b)". This should also provide a minor increase in efficiency.

It may be desirable to generate the expression dynamically, e.g.

	offsets = [d
	    for j in xrange(1,6+1)
	    for i in [j,-j]
	    for d in [(i,0),(0,i),(i,i),(i,-i)]]
	terms = ["(eudem_osttirol[%d,%d] < eudem_osttirol)" % d
	         for d in offsets]
	expr = "elevation_percentile_step2 = (100.0 / 48.0) * \\\n(%s)" % " + \\\n".join(terms)

This makes it easier to change e.g. the size of the window or the name
of the map.

But if you want any significant gain in efficiency, the obvious
solution is to modify r.neighbors to add a new aggregate (similar to
interspersion, but using less-than rather than not-equal). The weight=
option can be used to specify the window shape (for an aggregate which
lacks a weighted version, the weights are converted to a mask, where
cells with a non-zero weight are used and those with a zero weight
ignored).

If the region is small, you may get away with using grass.script.array
and NumPy. E.g. (untested)

	import numpy
	import grass.script.array
	data = grass.script.array.array()
	data.read("eudem_osttirol")
	offsets = [d
	    for j in xrange(1,6+1)
	    for i in [j,-j]
	    for d in [(i,0),(0,i),(i,i),(i,-i)]]
	rows,cols = data.shape
	center = data[6:,6:][:rows-13,:cols-13]
	count = numpy.zeros((rows-13,cols-13),dtype=int)
	for dy,dx in offsets:
	    cmp = data[6-dy:,6-dx:][:rows-13,:cols-13] < center
	    count += cmp
	perc = grass.script.array.array()
	perc[...] = 0
	perc[6:-7,6:-7] = count * 100.0 / 48
	perc.write("elevation_percentile_step2")

Realistically, this requires that the raw data is significantly
smaller than the available memory.

This type of problem comes up often enough that we really want an
aggregate (method) which evaluates an arbitrary expression with the
cell values and weights as inputs.

The language would need to support iterating over arrays (which rules
out r.mapcalc), reasonably easy to embed, and preferably efficient for
the case where you're evaluating the same expression many times with
different inputs (i.e. it shouldn't parse/compile/etc the expression
each time).

-- 
Glynn Clements <glynn at gclements.plus.com>