[GRASS-dev] Compute mahalanobis distance using Scipy
Paulo van Breugel
p.vanbreugel at gmail.com
Thu Jan 29 09:30:16 PST 2015
I would like to compute a raster layer with for each raster cell the
mahalanobis distance to the centre of the environmental space formed by all
reference data points (raster cells). In R this can be done as explained
here [1].
. 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)
ref = ['bio1 at clim','bio2 at clim','bio3 at clim']
import numpy as np
from scipy.spatial import distance
import scipy.linalg as linalg
import grass.script as grass
import grass.script.array as garray
# Create covariance table (could do this in python instead)
text_file = open("tmpfile", "w")
text_file.write(grass.read_command("r.covar", quiet=True, map=ref))
text_file.close()
covar = np.loadtxt(tmpcov, skiprows=1)
os.remove(tmpcov)
# Import data
dat_ref = None
stat_mean = None
layer = garray.array()
for i, map in enumerate(ref):
layer.read(map, null=np.nan)
s = len(ref)
r, c = layer.shape
if dat_ref is None:
dat_ref = np.empty((s, r, c), dtype=np.double)
if stat_mean is None:
stat_mean = np.empty((s), dtype=np.double)
dat_ref[i,:,:] = layer
stat_mean[i] = np.nanmean(layer)
del layer
# Compute mahalanobis distance
r, c = dat_ref[1,:,:].shape
stat_mah = garray.array()
for i in xrange(r):
for j in xrange(c):
cell_ref = dat_ref[...,i,j]
stat_mah[i,j] = distance.mahalanobis(cell_ref, stat_mean
, linalg.inv(covar))
stat_mah.write("mahalanobisMap")
[1]
https://pvanb.wordpress.com/2014/05/13/a-new-method-and-tool-exdet-to-evaluate-novelty-environmental-conditions/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-dev/attachments/20150129/7719c6b1/attachment.html>
More information about the grass-dev
mailing list