On Tue, Feb 17, 2015 at 7:22 AM, Pietro <peter.z...@gmail.com> wrote:
> On Mon, Feb 16, 2015 at 10:50 PM, Paulo van Breugel > <p.vanbreu...@gmail.com> wrote: > > On Sat, Feb 14, 2015 at 7:29 PM, Vaclav Petras <wenzesl...@gmail.com> > wrote: > >> On Sat, Feb 14, 2015 at 11:47 AM, Paulo van Breugel > >> <p.vanbreu...@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 >
_______________________________________________ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev