<div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Feb 17, 2015 at 7:22 AM, Pietro <span dir="ltr"><<a href="mailto:peter.zamb@gmail.com" target="_blank">peter.zamb@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">On Mon, Feb 16, 2015 at 10:50 PM, Paulo van Breugel<br>
<span><<a href="mailto:p.vanbreugel@gmail.com" target="_blank">p.vanbreugel@gmail.com</a>> wrote:<br>
> On Sat, Feb 14, 2015 at 7:29 PM, Vaclav Petras <<a href="mailto:wenzeslaus@gmail.com" target="_blank">wenzeslaus@gmail.com</a>> wrote:<br>
>> On Sat, Feb 14, 2015 at 11:47 AM, Paulo van Breugel<br>
>> <<a href="mailto:p.vanbreugel@gmail.com" target="_blank">p.vanbreugel@gmail.com</a>> wrote:<br>
>>> For a quick solution, what about using r.tile to split the input data in<br>
>>> tiles and compute the mahalanobis distance per tile.<br>
>><br>
>> See PyGRASS GridModule class which will do tiling (and a lot of other<br>
>> things) for you.<br>
>><br>
>> <a href="http://grass.osgeo.org/grass70/manuals/libpython/pygrass.modules.grid.html" target="_blank">http://grass.osgeo.org/grass70/manuals/libpython/pygrass.modules.grid.html</a><br>
><br>
> The Grid Module does indeed what I want, except that I want to run a python<br>
> function (the mahalanobis distance computation) on the tiles, while the<br>
> GridModule requires a raster GRASS r.* command. Any way around that? Writing<br>
> an intermediate addon/function that I can then call with the GridModule?<br>
<br>
</span>mhh the GridModule it is heavily using the Module interface, therefore<br>
I think you should provide a module and then use the GridModule.<br>
<br>
Therefore I wrote a function that do exactly what you are looking for:<br>
<br>
{{{<br>
import numpy as np<br>
<br>
from grass.pygrass.gis.region import Region<br>
from grass.pygrass.raster import RasterRow<br>
from grass.pygrass.raster.buffer import Buffer<br>
from grass.pygrass.utils import split_in_chunk<br>
<br>
<br>
def mahalanobis_distance(array, param):<br>
    # compute mahalanobis distance, here just a dummy computation<br>
    dist = array / param<br>
    return dist<br>
<br>
<br>
def tiled_function(raster_input, raster_output, func,<br>
                   out_mtype='DCELL', overwrite=False, nrows=100,<br>
                   **kwargs):<br>
    reg = Region()<br>
    with RasterRow(raster_input, mode='r') as rinp:<br>
        with RasterRow(raster_output, mode='w', mtype=out_mtype) as rout:<br>
            buf = Buffer((reg.cols, ), mtype=out_mtype)  # take a row buffer<br>
            for chunk in split_in_chunk(rinp, nrows):<br>
                ichunk = np.array(chunk)<br>
                for row in func(ichunk, **kwargs):<br>
                    buf[:] = row[:]<br>
                    rout.put_row(buf)<br>
<br>
<br>
tiled_function('elevation', 'distance', mahalanobis_distance,<br>
               overwrite=True, param=2.)<br>
}}}<br>
<br></blockquote><div><br></div><div>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).<br><br>{{{<br>

<p style="margin:0px;text-indent:0px">def mahalanobis_distances(X, m, L):</p>
<p style="margin:0px;text-indent:0px">       cX = X - m[np.newaxis, :]</p>
<p style="margin:0px;text-indent:0px">       tmp = solve_triangular(L, cX.T, lower=True).T</p>
<p style="margin:0px;text-indent:0px">       tmp **= 2</p>}}}<br></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
I'm not sure where should I put this function, maybe in pygrass.raster<br>
or in pygrass.utils (but the last one it is a bit generic...).<br>
<br>
All the best<br>
<span><font color="#888888"><br>
Pietro<br>
</font></span></blockquote></div><br></div></div>