<div dir="ltr">

<p style="margin:0px;text-indent:0px">I would like to compute a raster layer with for each raster cell the mahalanobis distance to the centre of the environmental space<span style="color:rgb(153,51,0)"><em></em></span>  formed by all reference data points (raster cells). In R this can be done as explained here [1]. <br></p><p style="margin:0px;text-indent:0px">. I would like to do this using python only (no dependency on R). I came up with the following, which works, but is very slow. I guess this is because it loops over every raster cell to compute the mahalanobis distance? Any idea how this can be done faster (I am very new to python, so bear with me if I am making stupid mistakes)<br></p><br><p style="margin:0px;text-indent:0px"><br></p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace"></span></p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">ref = ['bio1@clim','bio2</span><span style="font-family:monospace,monospace"><span style="font-family:monospace,monospace">@clim</span>','bio3</span><span style="font-family:monospace,monospace"><span style="font-family:monospace,monospace">@clim</span>']<br></span></p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace"></span><span style="font-family:monospace,monospace"></span></p><p style="margin:0px;text-indent:0px"></p><p style="margin:0px;text-indent:0px">

</p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">import numpy as np</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">from scipy.spatial import distance</span></p><p style="margin:0px;text-indent:0px">

</p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">import scipy.linalg as linalg<br></span></p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">import grass.script as grass</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">import grass.script.array as garray</span></p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace"><br></span></p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace"># Create covariance table (could do this in python instead)<br></span></p><span style="font-family:monospace,monospace">
</span><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    text_file = open("tmpfile", "w")</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    text_file.write(grass.read_command("r.covar", quiet=True, map=ref))</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    text_file.close()</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    covar = np.loadtxt(tmpcov, skiprows=1)</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    </span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    os.remove(tmpcov)</span></p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace"> </span></p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace"># Import data<br></span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    dat_ref = None</span></p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">stat_mean = None<br></span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    layer = garray.array()</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace"><br></span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    for i, map in enumerate(ref):</span></p><span style="font-family:monospace,monospace"></span><span style="font-family:monospace,monospace">    layer.read(map, null=np.nan)</span><br><span style="font-family:monospace,monospace"></span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    s = len(ref)</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    r, c = layer.shape</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    if dat_ref is None:</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">        dat_ref = np.empty((s, r, c), dtype=np.double)</span></p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    if stat_mean is None:</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">        stat_mean = np.empty((s), dtype=np.double)</span></p><span style="font-family:monospace,monospace">


</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    dat_ref[i,:,:] = layer</span></p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    stat_mean[i] = np.nanmean(layer)</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px">        </p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">del layer</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace"><br></span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    # Compute mahalanobis distance</span></p><span style="font-family:monospace,monospace">

</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    r, c = dat_ref[1,:,:].shape</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    stat_mah = garray.array()</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">    for i in xrange(r):</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">            for j in xrange(c):</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">                    cell_ref = dat_ref[...,i,j]</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">                    stat_mah[i,j] = distance.mahalanobis(cell_ref, stat_mean</span></p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">, linalg.inv(covar))</span></p><span style="font-family:monospace,monospace">
</span><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">stat_mah.write("mahalanobisMap")</span></p><span style="font-family:monospace,monospace"></span><p style="margin:0px;text-indent:0px"><br></p><p style="margin:0px;text-indent:0px"><br></p><p style="margin:0px;text-indent:0px">[1] <a href="https://pvanb.wordpress.com/2014/05/13/a-new-method-and-tool-exdet-to-evaluate-novelty-environmental-conditions/">https://pvanb.wordpress.com/2014/05/13/a-new-method-and-tool-exdet-to-evaluate-novelty-environmental-conditions/</a></p>
<p style="margin:0px;text-indent:0px">    </p></div>