"""
Testtext for documentation
contains :py:func:`get_best_LC`
"""
from collections import defaultdict
import numpy as np
from ppcpy.misc.helper import mean_stable
from ppcpy.misc.helper import default_to_regular
elastic2raman = {355: 387, 532: 607}
[docs]
def lc_for_cldFreeGrps(data_cube, retrieval:str) -> list:
"""Estimate the lidar constant from the optical profiles.
Parameters
----------
data_cube : object
Main PicassoProc object.
retrieval : str
Retrieval type. 'klett' or 'raman'.
Returns
-------
LCs : list
Lidar constant for retrieval type per channel per cloud free period.
Notes
-----
- For NR, done directly form the optical profiles, whereas in the matlab version, the ``LC*olAttri387.sigRatio`` is taken.
.. TODO:: Change back to Picasso version to check if lidar calibration constants get more similar.
.. TODO:: Check if LC's are normalized with respect to the mean of the profiles.
"""
print("retival", retrieval)
height = data_cube.retrievals_highres['range']
hres = data_cube.rawdata_dict['measurement_height_resolution']['var_data']
config_dict = data_cube.polly_config_dict
heightFullOverlap = [np.array(config_dict['heightFullOverlap']) for i in data_cube.clFreeGrps]
print('LCMeanWindow', config_dict['LCMeanWindow'],
'LCMeanMinIndx', config_dict['LCMeanMinIndx'],
'LCMeanMaxIndx', config_dict['LCMeanMaxIndx'])
LCs = defaultdict(list)
for i, cldFree in enumerate(data_cube.clFreeGrps):
cldFreeTime = np.array(data_cube.retrievals_highres['time'])[cldFree]
cldFree = cldFree[0], cldFree[1] + 1
profiles = data_cube.retrievals_profile[retrieval][i]
for channel in profiles:
wv, t, tel = channel.split('_')
# Telescope type dependent configurations:
if tel == 'NR':
key_smooth = f'smoothWin_{retrieval}_NR_'
else:
key_smooth = f'smoothWin_{retrieval}_'
hFullOverlap = heightFullOverlap[i][data_cube.gf(wv, t, tel)][0]
hBaseInd = np.argmax(
height >= (hFullOverlap + config_dict[f'{key_smooth}{wv}'] / 2 * hres))
# Elastic signals:
sig = profiles[channel]['signal']
signal = np.nanmean(np.squeeze(
data_cube.retrievals_highres[f'sig{sig}'][slice(*cldFree), :, data_cube.gf(wv, t, tel)]), axis=0)
molBsc = data_cube.mol_profiles[f'mBsc_{wv}'][i, :]
molExt = data_cube.mol_profiles[f'mExt_{wv}'][i, :]
# Check for avaiabel retrievals:
if not ('aerExt' in profiles[channel] and 'aerBsc' in profiles[channel]):
print(f'No availabel retrievals, skipping {channel} {cldFree}')
continue
aerExt = profiles[channel]['aerExt'].copy()
aerExt[:hBaseInd] = aerExt[hBaseInd]
aerBsc = profiles[channel]['aerBsc']
aerOD = np.nancumsum(aerExt * np.concatenate(([height[0]], np.diff(height))))
molOD = np.nancumsum(molExt * np.concatenate(([height[0]], np.diff(height))))
trans = np.exp(-2 * (aerOD + molOD))
bsc = molBsc + aerBsc
LC = signal * height**2 / bsc / trans
LC_stable, _, LCStd = mean_stable(
x=LC,
win=config_dict['LCMeanWindow'],
minBin=config_dict['LCMeanMinIndx'],
maxBin=config_dict['LCMeanMaxIndx']
)
print(f"cldFreGrp {i}, Channel {wv} {t} {tel}, LC_stable {LC_stable}, LCStd {LCStd}")
if LC_stable is None:
print(f'Can not find a stable LC value, skipping {channel} {cldFree}')
continue
LCs[channel].append({
'LC': LC_stable, 'LCStd': LC_stable * LCStd,
'time_start': int(cldFreeTime[0]), 'time_end': int(cldFreeTime[1])
})
if retrieval == 'raman' and int(wv) in elastic2raman.keys():
wv_r = elastic2raman[int(wv)]
# Inelastic signals:
signal_r = np.nanmean(np.squeeze(
data_cube.retrievals_highres[f'sig{sig}'][slice(*cldFree), :, data_cube.gf(wv_r, t, tel)]), axis=0)
#molBsc_r = data_cube.mol_profiles[f'mBsc_{wv_r}'][i, :]
molExt_r = data_cube.mol_profiles[f'mExt_{wv_r}'][i, :]
aerExt_r = aerExt * (int(wv)/int(wv_r))**config_dict['angstrexp']
aerOD_r = np.cumsum(aerExt_r * np.concatenate(([height[0]], np.diff(height))))
molOD_r = np.cumsum(molExt_r * np.concatenate(([height[0]], np.diff(height))))
trans_r = np.exp(- (aerOD + molOD + aerOD_r + molOD_r))
bsc = molBsc
LC_r = signal_r * height**2 / bsc / trans_r
LC_r_stable, _, LCStd_r = mean_stable(
x=LC_r,
win=config_dict['LCMeanWindow'],
minBin=config_dict['LCMeanMinIndx'],
maxBin=config_dict['LCMeanMaxIndx']
)
print(f"cldFreGrp {i}, Channel {wv_r} {t} {tel}, LC_stable {LC_r_stable}, LCStd {LCStd_r}")
if LC_r_stable is None:
print(f'Can not find a stable LC value, skipping {channel} {cldFree}')
continue
LCs[f"{wv_r}_{t}_{tel}"].append({
'LC': LC_stable, 'LCStd': LC_stable * LCStd,
'time_start': int(cldFreeTime[0]), 'time_end': int(cldFreeTime[1])
})
return default_to_regular(LCs)