# [GRASS-dev] Compute mahalanobis distance using Scipy

Pietro peter.zamb at gmail.com
Tue Feb 17 02:31:09 PST 2015

```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)
}}}
```