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@clim’,‘bio2@clim’,‘bio3@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”)