import numpy as np
from libpyhat.derived.utils import compute_b_a, reflectance
[docs]
def bd_func1(data, wvls, kernel):
wvl_left = wvls[0]
wvl_center = wvls[1]
wvl_right = wvls[2]
b, a = compute_b_a([wvl_left, wvl_center, wvl_right])
r_left = reflectance(data, wvl_left, kernel=kernel[0])
r_center = reflectance(data, wvl_center, kernel=kernel[1])
r_right = reflectance(data, wvl_right, kernel=kernel[2])
return 1.0 - (r_center / ((a * r_left) + (b * r_right)))
[docs]
def bd_func2(data, wvls, kernel):
wvl_left = wvls[0]
wvl_center = wvls[1]
wvl_right = wvls[2]
b, a = compute_b_a([wvl_left, wvl_center, wvl_right])
r_left = reflectance(data, wvl_left, kernel=kernel[0])
r_center = reflectance(data, wvl_center, kernel=kernel[1])
r_right = reflectance(data, wvl_right, kernel=kernel[2])
return 1.0 - (r_center / ((a * r_left) + (b * r_right)))
[docs]
def min_2_bands(data, wvls1, kernel1, wvls2, kernel2):
band_wvls1 = data.closest_wvl(wvls1)
bd1 = bd_func2(data, band_wvls1, kernel1)
band_wvls2 = data.closest_wvl(wvls2, kernel2)
bd2 = bd_func2(data, band_wvls2)
return np.min([bd1, bd2])
[docs]
def sh_func(data, wvls, kernel):
wvl_left = wvls[0]
wvl_center = wvls[1]
wvl_right = wvls[2]
b, a = compute_b_a([wvl_left, wvl_center, wvl_right])
r_left = reflectance(data, wvl_left, kernel=kernel[0])
r_center = reflectance(data, wvl_center, kernel=kernel[1])
r_right = reflectance(data, wvl_right, kernel=kernel[2])
return 1.0 - (((a * r_left) + (b * r_right)) / r_center)
[docs]
def rpeak1_func(data, vnir_wvls, threshold=0.000001):
vnir_spectra = data.df[data.spect_label][vnir_wvls]
wavelength = np.array(vnir_spectra.columns.values, dtype=float)
rpeaks = []
print(
"Calculating rpeak1. This involves a polynomial fit to each "
"spectrum, and may take a while on large cubes..."
)
for i in vnir_spectra.index:
try:
fit_result = np.polynomial.Polynomial.fit(
wavelength, np.array(vnir_spectra.loc[i]), 5
)
deriv = fit_result.deriv()
rpeaks.append(
vnir_spectra.loc[i][np.min(wavelength[deriv(wavelength) < threshold])]
)
except:
rpeaks.append(0)
return np.array(rpeaks, dtype=float)
[docs]
def crism_sumutil_bdicont(data, low_wvl1, high_wvl1, low_wvl2, high_wvl2):
# Adapted from CRISM CAT crism_sumutil_bdicont, assuming hyperspectral data
wvls_1um = data.wvls[(data.wvls > low_wvl1) * (data.wvls < high_wvl1)]
pcnt75_wvls = []
pcnt75 = np.percentile(
data.df[data.spect_label][wvls_1um], 75, interpolation="nearest", axis=1
)
data_1um = data.df[data.spect_label][wvls_1um]
print(
"Iterating over each spectrum to find wvl of the 75th percentile "
"reflectance between " + str(low_wvl1) + " and " + str(high_wvl1)
)
print("This can take a few minutes!")
for i in range(data.df.shape[0]): # is there a faster way to do this?
# Probably.
pcnt75_wvls.append(wvls_1um[data_1um.iloc[i, :] == pcnt75[i]][0])
wvls_2p4 = data.wvls[(data.wvls > low_wvl2) * (data.wvls < high_wvl2)]
median_2p4 = np.median(data.df[data.spect_label][wvls_2p4], axis=1)
median_2p4_wvl = np.median(wvls_2p4)
cont_slope = (median_2p4 - pcnt75) / (median_2p4_wvl - pcnt75_wvls)
return cont_slope, median_2p4, median_2p4_wvl, pcnt75, pcnt75_wvls