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 <gl...@gclements.plus.com> _______________________________________________ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev