[GRASS-dev] Compute mahalanobis distance using Scipy

[Pietro]

On Mon, Feb 16, 2015 at 10:50 PM, Paulo van Breugel
<p.vanbreugel at gmail.com> wrote:
> On Sat, Feb 14, 2015 at 7:29 PM, Vaclav Petras <wenzeslaus at gmail.com> wrote:
>> On Sat, Feb 14, 2015 at 11:47 AM, Paulo van Breugel
>> <p.vanbreugel at 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