[GRASS-dev] Compute mahalanobis distance using Scipy

[Paulo van Breugel]

On Tue, Feb 17, 2015 at 7:22 AM, Pietro <peter.zamb at gmail.com> wrote:

> 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.)
> }}}
>
>
Can raster_input be multiple rasters? The input of the mahalanobis distance
requires as input an ndarray with values of multiple rasters, i.e. X below.
m is an array with average values per raster, L is the covariance matrix of
the input rasters (both m and L can of course be created using standard
grass functions).

{{{

def mahalanobis_distances(X, m, L):

    cX = X - m[np.newaxis, :]

    tmp = solve_triangular(L, cX.T, lower=True).T

    tmp **= 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-dev/attachments/20150217/b0b6a30d/attachment.html>