On Mon, Feb 16, 2015 at 10:50 PM, Paulo van Breugel
<p.vanbreu...@gmail.com> wrote:
> On Sat, Feb 14, 2015 at 7:29 PM, Vaclav Petras <wenzesl...@gmail.com> wrote:
>> On Sat, Feb 14, 2015 at 11:47 AM, Paulo van Breugel
>> <p.vanbreu...@gmail.com> wrote:
>>> For a quick solution, what about using r.tile to split the input data in
>>> tiles and compute the mahalanobis distance per tile.
>>
>> See PyGRASS GridModule class which will do tiling (and a lot of other
>> things) for you.
>>
>> http://grass.osgeo.org/grass70/manuals/libpython/pygrass.modules.grid.html
>
> The Grid Module does indeed what I want, except that I want to run a python
> function (the mahalanobis distance computation) on the tiles, while the
> GridModule requires a raster GRASS r.* command. Any way around that? Writing
> an intermediate addon/function that I can then call with the GridModule?
mhh the GridModule it is heavily using the Module interface, therefore
I think you should provide a module and then use the GridModule.
Therefore I wrote a function that do exactly what you are looking for:
{{{
import numpy as np
from grass.pygrass.gis.region import Region
from grass.pygrass.raster import RasterRow
from grass.pygrass.raster.buffer import Buffer
from grass.pygrass.utils import split_in_chunk
def mahalanobis_distance(array, param):
# compute mahalanobis distance, here just a dummy computation
dist = array / param
return dist
def tiled_function(raster_input, raster_output, func,
out_mtype='DCELL', overwrite=False, nrows=100,
**kwargs):
reg = Region()
with RasterRow(raster_input, mode='r') as rinp:
with RasterRow(raster_output, mode='w', mtype=out_mtype) as rout:
buf = Buffer((reg.cols, ), mtype=out_mtype) # take a row buffer
for chunk in split_in_chunk(rinp, nrows):
ichunk = np.array(chunk)
for row in func(ichunk, **kwargs):
buf[:] = row[:]
rout.put_row(buf)
tiled_function('elevation', 'distance', mahalanobis_distance,
overwrite=True, param=2.)
}}}
I'm not sure where should I put this function, maybe in pygrass.raster
or in pygrass.utils (but the last one it is a bit generic...).
All the best
Pietro
_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev
