[GRASS-dev] Compute mahalanobis distance using Scipy

[Pietro]

Hi,

On Tue, Feb 17, 2015 at 10:44 AM, Paulo van Breugel
<p.vanbreugel at gmail.com> wrote:
> Can raster_input be multiple rasters?

yes, I've define a new split_rasters function, that should does what
you want, note, the code is not tested, just to provide a rough idea.

{{{
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_distances(X):
    m = np.mean(X, axis=0)
    L = np.cov(m)
    cX = X - m[np.newaxis, :]
    tmp = solve_triangular(L, cX.T, lower=True).T
    tmp **= 2

def split_rasters(rasts, nrows=100):
    generators = [split_in_chunk(rst, nrows) for rst in rasts]
    for chunk in zip(*generators):
        yield chunk


def tiled_function(raster_inputs, raster_output, func,
                   out_mtype='DCELL', overwrite=False, nrows=100,
                   **kwargs):
    reg = Region()
    rasts = [RasterRow(rst) for rst in raster_inputs]
    # open input raster maps
    for rst in rasts:
        rst.open('r')
    # open the output map
    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_rasters(rasts, nrows):
            ichunk = np.array(chunk)
            for row in func(ichunk, **kwargs):
                buf[:] = row[:]
                rout.put_row(buf)


tiled_function('elevation', 'distance', mahalanobis_distances,
               overwrite=True)
}}}