Bingham statistics

[Antoine-Cate]

Only Fisher stats are available in stereonet_math.py, and these are not ideal to calculate the mean value of a series of planes or lines.
A new function (bingham_stats) that would provide the orientation of the three axes and the Bingham best fit plane would be beneficial. I need to calculate these for a project, so I hope to be able to propose a solution to this issue within the next two months.

[joferkington]

Have you had a look at various functionality in https://github.com/joferkington/mplstereonet/blob/master/mplstereonet/analysis.py ?

For orientations, the axes are identical to what you'd get by estimating the full Bingham distribution. (You actually have to use the same method to estimate the axes when fitting a Bingham distribution -- it's still eigenvectors of the covariance matrix of orientations to get the axes.) You don't get the additional shape parameters of the distribution, of course, but the fitting methods there may very well do what you need.

That having been said, I'd welcome a bingham_stats function, though it should go in analysis instead of stereonet_math, as fisher_stats should have originally.

[stims99]

Ugly piece of code but this does bingham stats:
def eigvector(w, v, list):
esort = np.argsort(w)

for c in esort:
    
    t = round(np.degrees(math.atan(v[0,c]/v[1,c]) + math.pi if v[1,c] < 0 else math.atan(v[0,c]/v[1,c])),2)
    t = t if t > 0 else t + 360
    p = round(np.degrees(math.asin(v[2,c])),2)
    
    
    if p < 0:
        dip = 90 + p
        dipdir = t
    else:
        dip = 90 - p
        dipdir = t + 180 if t < 180 else t - 180

    list.append([dip,dipdir])

def dir_cos(value):
lm = np.inner(value.Jx, value.Jy)
ln = np.inner(value.Jx, value.Jz)
mn = np.inner(value.Jy, value.Jz)
l2 = np.inner(value.Jx, value.Jx)
m2 = np.inner(value.Jy, value.Jy)
n2 = np.inner(value.Jz, value.Jz)

return np.array([l2,lm,ln, lm,m2,mn, ln, mn, n2]).reshape(3,3)

#initialize empty dictionaries
dict_mean = dict()
#dict_ori is a dictionary for each joint set
#loop through each joint set and compute bingham stats

for key, value in dict_ori.items():
value = value.reset_index(drop = True)
#directional cosine matric
dcm = sp.dir_cos(value)

#hermitian matrix to return eigenvalues
w, v = np.linalg.eigh(dcm, UPLO = 'L')
#compute the unit eigenvalues, makes converting to dip/dip azimuth easier
w_unit = np.around(w/np.sum(w),decimals =3)

#make sure eigenvectors sorted small to large (e3, e2, e1)
esort = np.argsort(w)

#initialize an empty list
ori_list = []

#calculate the eigenvector dip and dip azimuth (in degrees). rotate/flip if plotting in upper hemisphere
sp.eigvector(w, v, ori_list)

ori_list.append([w_unit.tolist(), len(value)])
ori_flat = [item for sublist in ori_list for item in sublist]

dict_mean[key] = ori_flat

ori_cols = ['e3-Dip', 'e3-Dip Az.', 'e2-Dip', 'e2-Dip Az.', 'e1-Dip', 'e1-Dip Az.', 'Eigenvalues', 'Count']

df_mean = pd.DataFrame.from_dict(dict_mean, orient = 'index', columns = ori_cols)
df_mean.sort_index()