Bhattacharyya distance, Euclidean distance, and Spectral angle calculator for multivariate samples

distances.py

""" - Bhattacharyya distance, Euclidean distance, and Spectral angle calculator for multivariate samples - 

It assumes that there are the same amount of non-valid values in every band.

Input files are multiband raster layers

Alfonso Diz-Lois <dizlois@gmail.com>"""



def Calculate_Distances(file1,file2,nondatavalue=-9999):
    from osgeo import gdal
    import numpy as np
    
    
    #Open files and check the number of bands...
    ga=gdal.Open(file1)
    gb=gdal.Open(file2)
    Bands=ga.RasterCount
    if not Bands==gb.RasterCount:
        print "Error found: different number of bands in %s and %s"%(file1,file2)
        sys.exit()
    signatures=["a","b"]
    #Put all the bands together in a (bands*features) array
    for z in signatures:
        exec "band1=g%s.GetRasterBand(1).ReadAsArray()"%(z)
        exec "dist_%s=np.zeros((Bands,len(band1[np.where(band1!=%s)])))"%(z,nondatavalue)
        for i in xrange(Bands):
            exec "%s = g%s.GetRasterBand(i+1).ReadAsArray()" % (("Band"+str(i+1)),z)
            exec ("%s =np.ravel(Band"+str(i+1)+")") %("list"+str(i+1)+str(z))
            exec ("%s =[ x for x in %s if x!=%s ]") %("list"+str(i+1)+str(z),"list"+str(i+1)+str(z),nondatavalue)
            exec "dist_%s[i,:]=%s"%(z,"list"+str(i+1)+str(z))
    
    #Now, we need to get the mean and covariance matrix for each distribution...
    
    mean1=dist_a.mean( axis=1 )
    mean2=dist_b.mean( axis=1 )
    cov1 = np.cov (dist_a)
    cov2 = np.cov (dist_b)

    # We calculate the Bhattacharyya distance as follow...
    
    cov_mean=(cov1+cov2)/2
    inv_cov_mean=np.linalg.inv ( cov_mean )
    s = np.subtract(mean1,mean2)
    det1=np.linalg.det(cov1)
    det2=np.linalg.det(cov2)
    detmean=np.linalg.det(cov_mean)    
    bat_distance=(0.125)*np.dot(np.dot(s.T,inv_cov_mean),s)+(0.5)*np.log(detmean/np.sqrt(det1*det2))

    print "Bhattacharyya Distance= ", bat_distance
    
    #Euclidean distance...
    
    euc_distance=np.linalg.norm(mean1-mean2)
    
    print "Euclidean Distance= ", euc_distance
    
    #And Spectral Angle...
    
    xyterm=0.
    x_quad=0.
    y_quad=0.
    for i in xrange(len(mean1)):
        xyterm+=mean1[i]*mean2[i]
        x_quad+=(mean1[i])**2
        y_quad+=(mean2[i])**2
    Spectral_angle=np.rad2deg((np.arccos(xyterm/((np.sqrt(x_quad))*(np.sqrt(y_quad))))))
    print "Spectral angle= ", Spectral_angle
    return (bat_distance,euc_distance,Spectral_angle)