From 723b6f0d3d828761f67e994ef803d06fe1aa213e Mon Sep 17 00:00:00 2001 From: HavardStridBuholdt Date: Tue, 11 Aug 2026 17:52:41 +0200 Subject: [PATCH] Refine calibration constant selection. --- ppcpy/calibration/lidarconstant.py | 179 ++++-- ppcpy/calibration/polarization.py | 570 +++++++++++------- ppcpy/calibration/select.py | 59 +- ppcpy/config/json2nc-mapper_NR_att_bsc.json | 4 +- ppcpy/config/json2nc-mapper_OC_att_bsc.json | 6 +- ppcpy/config/json2nc-mapper_att_bsc.json | 6 +- ppcpy/config/json2nc-mapper_profiles.json | 24 +- .../config/json2nc-mapper_quasi_results.json | 10 +- .../json2nc-mapper_quasi_results_V2.json | 10 +- ppcpy/config/json2nc-mapper_vol_depol.json | 6 +- ppcpy/interface/picassoProc.py | 284 +++++++-- ppcpy/io/sql_interaction.py | 47 +- ppcpy/io/write2nc.py | 4 +- ppcpy/qc/transCor.py | 5 +- ppcpy/retrievals/depolarization.py | 8 +- ppcpy/retrievals/highres.py | 7 +- ppcpy/retrievals/quasi.py | 2 +- 17 files changed, 851 insertions(+), 380 deletions(-) diff --git a/ppcpy/calibration/lidarconstant.py b/ppcpy/calibration/lidarconstant.py index 8ccb2d2..34d8ded 100644 --- a/ppcpy/calibration/lidarconstant.py +++ b/ppcpy/calibration/lidarconstant.py @@ -12,8 +12,81 @@ elastic2raman:dict = {355: 387, 532: 607} -def lc_for_cldFreeGrps(data_cube, retrieval:str, collect_debug:bool=False) -> list: - """Estimate the lidar constant from the optical profiles. + +def loadDefaults(data_cube, **defaults) -> dict: + """Prepare default Lidar calibration values. + + Parameters + ---------- + data_cube : object + Main PicassoProc object. + LC : list, optional + Default Lidar constant value per channel. + LCStd : list, optional + Default Lidar constant error per channel. + + Returns + ------- + defaultDict : dict + Default Lidar calibration result per channel. + + Each channel contains a list with one single sub-dict with entries: + + ``LC`` : float + Default Lidar calibration constant. + + ``LCStd`` : float + Default uncertainty of lidar calibration constant. + + ``method`` : str + Name of retrieval method. + + Notes + ----- + Default values are by standard taken from their config variable but can be + overwritten if passed as an input to this function. + The order of ``LC`` and ``LCStd`` must match the channel order in + ``data_cube.retrievals_highres['channel']``. + + .. TODO:: Consider allowing default values for a single channel to passed as input. + + **History** + + - 2026-08-07: First edition by Buholdt + + + Example + ------- + >> loadDefaults(data_cube, + LC=[1e13, 1, 1, 1, 4e14, 1, ...], + LCStd=[1e-3, 1, 1, 1, 2e-2, 1, ...] + ) + """ + + default_values = data_cube.polly_config_dict | defaults + default_LC = np.asarray(default_values['LC']) + default_LCStd = np.asarray(default_values['LCStd']) + + defaultDict = {} + channels = [ + (355, 'FR'), (532, 'FR'), (1064, 'FR'), + (387, 'FR'), (607, 'FR'), + (355, 'NR'), (532, 'NR'), + (387, 'NR'), (607, 'NR'), + ] + + for (wv, tel) in channels: + defaultDict[f"{wv}_total_{tel}"] = [{ + 'LC': float(np.squeeze(default_LC[data_cube.gf(wv, 'total', tel)])), + 'LCStd': float(np.squeeze(default_LCStd[data_cube.gf(wv, 'total', tel)])), + 'method': 'default' + }] + + return defaultDict + + +def lc_for_cldFreeGrps(data_cube, retrieval:str, collect_debug:bool=False) -> dict: + """Estimate the lidar calibration constant from the optical profiles. Parameters ---------- @@ -26,24 +99,48 @@ def lc_for_cldFreeGrps(data_cube, retrieval:str, collect_debug:bool=False) -> li Returns ------- - LCs : list - Lidar constant for retrieval type per channel per cloud free period. + LCs : dict + Lidar calibration results for ``retrieval`` retrieved optical profiles per channel. + + Each channel contains a list of sub-dicts with entries: + + ``LC`` : float + Lidar calibration constant. + + ``LCStd`` : float + Uncertainty of lidar calibration constant. + + ``time_start``, ``time_end`` : int + Start and stop times for successful calibration. + + ``method`` : str + Name of retrieval method. + + The number of elements in each list depends on the number of successful retrievals. Notes ----- - - For NR, done directly form the optical profiles, whereas in the matlab version, the ``LC*olAttri387.sigRatio`` is taken. - - Through the config variable 'flagUseRetrievedExt4LCCalc', the extinction used to calculate the LCs can be specified. - if 'flagUseRetrievedExt4LCCalc' is True the retrieved extinction will be used otherwise the extinction approximated by - the backscatter times the assumed lidar constant will be used. - - Missing Rotational Raman and Aeronet LC retrieval. + For NR channels, the LC is calculated directly form the optical profiles, whereas in the matlab version, + it is estimated by multiplying the respective FR LC with ``olAttri387.sigRatio``. + + The function uses the following configuration flags: + + - ``flagUseRetrievedExt4LCCalc``: If enabled the retrieved extinction when calculating the LCs. + If disabled the extinction will be estimated by the retrieved + backscatter times the assumed LR. + + + The options for Rotational Raman and Aeronet LC retrievals are currently missing. .. TODO:: Check if LC's are normalized with respect to the mean of the profiles. + .. TODO:: Add option for Aeronet and rotational Raman retrieved LC. **History** xxxx-xx-xx: First edition by ... 2026-03-18: Changed beta_mol for inelastic wavelengths and added the 'flagUseRetrievedExt4LCCalc' variable. + """ logging.info(f'LC retrieval: {retrieval} method') @@ -61,7 +158,7 @@ def lc_for_cldFreeGrps(data_cube, retrieval:str, collect_debug:bool=False) -> li for channel in profiles: wv, t, tel = channel.split('_') - # Telescope type dependent configurations: + ## Telescope type dependent configurations if tel == 'NR': key_smooth = f'smoothWin_{retrieval}_NR_' key_LR = 'LR_NR_' @@ -73,29 +170,29 @@ def lc_for_cldFreeGrps(data_cube, retrieval:str, collect_debug:bool=False) -> li hBaseInd = np.argmax( height >= (hFullOverlap + config_dict[f'{key_smooth}{wv}'] / 2 * hres)) - # Elastic signal: + ## Elastic signal 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, :].copy() molExt = data_cube.mol_profiles[f'mExt_{wv}'][i, :].copy() - # Check for avaiabel retrievals: + ## Check for available retrievals if not ('aerExt' in profiles[channel] and 'aerBsc' in profiles[channel]): - logging.warning(f'No availabel retrievals, skipping {channel} {cldFree}') + logging.warning(f'No available retrievals, skipping {channel} {cldFree}') continue - # Backscatter and extinction retrievals: + ## Backscatter and extinction retrievals aerBsc = profiles[channel]['aerBsc'].copy() if config_dict['flagUseRetrievedExt4LCCalc'] & ~config_dict['flagPicassoComparison']: - logging.info('Using Retrieved Exticntion') + logging.info("Using Retrieved Extinction") aerExt = profiles[channel]['aerExt'].copy() else: - logging.info('Using approximated Extinction') + logging.info("Using approximated Extinction") aerBsc[aerBsc <= 0] = np.nan aerExt = aerBsc * config_dict[f'{key_LR}{wv}'] - # Interpolate extinction to ground + ## Interpolate extinction to ground aerExt[:hBaseInd + 1] = aerExt[hBaseInd] ## Optical depth (OD) @@ -115,22 +212,22 @@ def lc_for_cldFreeGrps(data_cube, retrieval:str, collect_debug:bool=False) -> li minBin=config_dict['LCMeanMinIndx'], maxBin=config_dict['LCMeanMaxIndx'] ) - logging.info(f'cldFreGrp {i}, Channel {wv} {t} {tel}, LC_stable {LC_stable}, LCStd {LCStd}') + logging.info(f"cldFreGrp {i}, Channel {wv} {t} {tel}, LC_stable {LC_stable}, LCStd {LCStd}") if LC_stable is None: - logging.warning(f'Can not find a stable LC value, skipping {wv} nm {t} {tel} channel for cloud free period {cldFree}') + logging.warning(f"Can not find a stable LC value, skipping {wv} nm {t} {tel} channel for cloud free period {cldFree}") continue - + + ## save LC result + LCs[channel].append({ + 'LC': LC_stable, 'LCStd': LC_stable * LCStd, + 'time_start': int(cldFreeTime[0]), 'time_end': int(cldFreeTime[1]), + 'method': retrieval + }) + + ## Collect debug info if collect_debug: - LCs[channel].append({ - 'LC': LC_stable, 'LCStd': LC_stable * LCStd, 'LC_profile': LC, - 'time_start': int(cldFreeTime[0]), 'time_end': int(cldFreeTime[1]) - }) - else: - LCs[channel].append({ - 'LC': LC_stable, 'LCStd': LC_stable * LCStd, - 'time_start': int(cldFreeTime[0]), 'time_end': int(cldFreeTime[1]) - }) + LCs[channel][-1]['LC_profile'] = LC # ----------------------------------------------------------------------------------- # LC for raman / inelastic channels @@ -155,7 +252,7 @@ def lc_for_cldFreeGrps(data_cube, retrieval:str, collect_debug:bool=False) -> li if config_dict['flagPicassoComparison']: bsc_r = molBsc - ## Lidar clibration constant + ## Lidar calibration constant LC_r = (signal_r * height**2) / (bsc_r * trans_r) LC_r[LC_r <= 0] = np.nan LC_r_stable, _, LCStd_r = mean_stable( @@ -164,22 +261,22 @@ def lc_for_cldFreeGrps(data_cube, retrieval:str, collect_debug:bool=False) -> li minBin=config_dict['LCMeanMinIndx'], maxBin=config_dict['LCMeanMaxIndx'] ) - logging.info(f'cldFreGrp {i}, Channel {wv_r} {t} {tel}, LC_stable {LC_r_stable}, LCStd {LCStd_r}') + logging.info(f"cldFreGrp {i}, Channel {wv_r} {t} {tel}, LC_stable {LC_r_stable}, LCStd {LCStd_r}") if LC_r_stable is None: - logging.warning(f'Can not find a stable LC value, skipping {wv_r} nm {t} {tel} channel for cloud free period {cldFree}') + logging.warning(f"Can not find a stable LC value, skipping {wv_r} nm {t} {tel} channel for cloud free period {cldFree}") continue + ## Save LC result + LCs[f"{wv_r}_{t}_{tel}"].append({ + 'LC': LC_r_stable, 'LCStd': LC_r_stable * LCStd_r, + 'time_start': int(cldFreeTime[0]), 'time_end': int(cldFreeTime[1]), + 'method': retrieval + }) + + ## Collect debug info if collect_debug: - LCs[f"{wv_r}_{t}_{tel}"].append({ - 'LC': LC_r_stable, 'LCStd': LC_r_stable * LCStd_r, 'LC_profile': LC_r, - 'time_start': int(cldFreeTime[0]), 'time_end': int(cldFreeTime[1]) - }) - else: - LCs[f"{wv_r}_{t}_{tel}"].append({ - 'LC': LC_r_stable, 'LCStd': LC_r_stable * LCStd_r, - 'time_start': int(cldFreeTime[0]), 'time_end': int(cldFreeTime[1]) - }) + LCs[f"{wv_r}_{t}_{tel}"][-1]['LC_profile'] = LC_r return default_to_regular(LCs) diff --git a/ppcpy/calibration/polarization.py b/ppcpy/calibration/polarization.py index daf7729..b247f08 100644 --- a/ppcpy/calibration/polarization.py +++ b/ppcpy/calibration/polarization.py @@ -40,6 +40,73 @@ def smooth_signal(signal:np.ndarray, window_len:int) -> np.ndarray: return uniform_filter(signal, window_len) +def loadDefaults(data_cube, **defaults) -> dict: + """Prepare default Depol calibration values. + + Parameters + ---------- + data_cube : object + Main PicassoProc object. + polCaliEta355 : float, optional + Default Depol calibration constant at 355 nm. + polCaliEtaStd355 : float, optional + Default Depol calibration constant error at 355 nm. + polCaliEta532 : float, optional + Default Depol calibration constant at 532 nm. + polCaliEtaStd532 : float, optional + Default Depol calibration constant error at 532 nm. + polCaliEta1064 : float, optional + Default Depol calibration constant at 1064 nm. + polCaliEtaStd1064 : float, optional + Default Depol calibration constant error at 1064 nm. + + Returns + ------- + defaultDict : dict + Default polarization calibration result per wavelength. + + Each wavelength contains a list with one single sub-dict with entries: + + ``eta`` : float + Default Depol calibration constant. + + ``eta_std`` : float + Defaults uncertainty of Depol calibration constant. + + ``method`` : str + Name of retrieval method. + + Notes + ----- + Default values are by standard taken from their config variable but can be + overwritten if passed as an input to this function. + + **History** + + - 2026-08-07: First edition by Buholdt + + + Example + ------- + >> loadDefaults(data_cube, + polCaliEta355=46.7, + polCaliEtaStd355=2.7e-3 + ) + """ + + default_values = data_cube.polly_config_dict | defaults + defaultDict = {} + + for wv in [355, 532, 1064]: + defaultDict[f'{wv}_FR'] = [{ + 'eta': float(default_values[f'polCaliEta{wv}']), + 'eta_std': float(default_values[f'polCaliEtaStd{wv}']), + 'method': 'default' + }] + + return defaultDict + + def loadGHK(data_cube): """Prepare the GHK parameters, especially if given in TR convert them into GHK. @@ -47,22 +114,43 @@ def loadGHK(data_cube): ---------- data_cube : object Main PicassoProc object. + + Yields + ------ + data_cube.polly_config_dict : dict + Updated to parameters: + TR --> Removed + G --> filled and converted to array + H --> filled and converted to array + K --> filled and converted to array + voldepol_error_355 --> converted to array + voldepol_error_532 --> converted to array + voldepol_error_1064 --> converted to array + + Notes + ----- + .. TODO:: + - Write a proper docstring. + + ** History ** + + - xx-xx-xxxx: First edition by ... + """ - print('starting loadGHK') - #print('flag_532_total', flag_532_total_FR) - #print('flag_532_cross', flag_532_cross_FR) - - print('data_cube keys ', data_cube.__dict__.keys()) - G = np.array(data_cube.polly_config_dict['G']).astype(float) - H = np.array(data_cube.polly_config_dict['H']).astype(float) - K = np.array(data_cube.polly_config_dict['K']).astype(float) - #print(TR[flag_532_total_FR]) - #print(TR[flag_532_cross_FR]) - #if data_cube.polly_config_dict['H'][0] == -999: + logging.info("Starting loadGHK") + # print('flag_532_total', flag_532_total_FR) + # print('flag_532_cross', flag_532_cross_FR) + + G = np.asarray(data_cube.polly_config_dict['G'], dtype=float) + H = np.asarray(data_cube.polly_config_dict['H'], dtype=float) + K = np.asarray(data_cube.polly_config_dict['K'], dtype=float) + # print(TR[flag_532_total_FR]) + # print(TR[flag_532_cross_FR]) + # if data_cube.polly_config_dict['H'][0] == -999: if (np.all(np.isclose(H, -999))): - TR = np.array(data_cube.polly_config_dict['TR']).astype(float) - print('H is empty -> calculate parameters') + TR = np.asarray(data_cube.polly_config_dict['TR'], dtype=float) + logging.info('H is empty -> calculate parameters') K[data_cube.flag_355_total_FR] = 1.0 K[data_cube.flag_532_total_FR] = 1.0 @@ -83,49 +171,34 @@ def loadGHK(data_cube): H[data_cube.flag_1064_cross_FR] = onemx_onepx(TR[data_cube.flag_1064_cross_FR]) if np.any(data_cube.flag_532_total_NR): - print("GHK for 532NR") + logging.info("GHK for 532 NR") K[data_cube.flag_532_total_NR] = 1.0 G[data_cube.flag_532_total_NR] = 1.0 H[data_cube.flag_532_total_NR] = onemx_onepx(TR[data_cube.flag_532_total_NR]) - print("H", H[data_cube.flag_532_total_NR]) + logging.info(f"G: {G[data_cube.flag_532_total_NR]}, H {H[data_cube.flag_532_total_NR]}, K: {K[data_cube.flag_532_total_NR]}.") if np.any(data_cube.flag_532_cross_DFOV): - print("GHK for 532DFOV") + logging.info("GHK for 532 DFOV") K[data_cube.flag_532_cross_DFOV] = 1.0 G[data_cube.flag_532_cross_DFOV] = 1.0 H[data_cube.flag_532_cross_DFOV] = onemx_onepx(TR[data_cube.flag_532_cross_DFOV]) - print("H", H[data_cube.flag_532_cross_DFOV]) - print('TR', TR) + logging.info(f"G: {G[data_cube.flag_532_cross_DFOV]}, H: {H[data_cube.flag_532_cross_DFOV]}, K: {K[data_cube.flag_532_cross_DFOV]}.") + logging.info(f"TR: {TR}") else: - print("Using GHK from config file") - #print('TR', TR) - #print('TR from H', (1-H)/(1+H)) - print('G', G) - print('H', H) - print('K', K) + logging.info("Using GHK from config file") + # print('TR', TR) + # print('TR from H', (1-H)/(1+H)) + logging.info(f"G: {G}, H: {H}, K: {K}") data_cube.polly_config_dict.pop('TR', None) # remove the TR from the config to avoid inconsistencies - data_cube.polly_config_dict['G'] = np.array(G) - data_cube.polly_config_dict['H'] = np.array(H) - data_cube.polly_config_dict['K'] = np.array(K) - data_cube.polly_config_dict['voldepol_error_355'] = np.array(data_cube.polly_config_dict['voldepol_error_355']) - data_cube.polly_config_dict['voldepol_error_532'] = np.array(data_cube.polly_config_dict['voldepol_error_532']) - data_cube.polly_config_dict['voldepol_error_1064'] = np.array(data_cube.polly_config_dict['voldepol_error_1064']) + data_cube.polly_config_dict['G'] = np.asarray(G) + data_cube.polly_config_dict['H'] = np.asarray(H) + data_cube.polly_config_dict['K'] = np.asarray(K) + data_cube.polly_config_dict['voldepol_error_355'] = np.asarray(data_cube.polly_config_dict['voldepol_error_355']) + data_cube.polly_config_dict['voldepol_error_532'] = np.asarray(data_cube.polly_config_dict['voldepol_error_532']) + data_cube.polly_config_dict['voldepol_error_1064'] = np.asarray(data_cube.polly_config_dict['voldepol_error_1064']) -""" -.. TODO:: can this be removed? -"TR": [0.898, 1086, 1, 1, 1.45, 778.8, 1, 1, 1, 1, 1, 1, 1], -"G": [-999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999 ], -"H": [-999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999 ], -"K": [-999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999 ], -"voldepol_error_355": [0.003, 0, 0], -"voldepol_error_532": [0.004, 0, 0], -"voldepol_error_1064": [0.005, 0, 0], -""" - - - -def calibrateGHK(data_cube) -> dict: - """Estimate the polarization calibration from the delta 90 Method [1]_ +def calibrateGHK(data_cube, collect_debug:bool=False) -> dict: + """Estimate the polarization calibration from the Delta-90° method [1]_. Parameters ---------- @@ -135,22 +208,38 @@ def calibrateGHK(data_cube) -> dict: Returns ------- pol_cali : dict - polarization factors from delta 90 for each wavelength containing - sub-dicts with 'eta', 'eta_std', 'time_start', 'time_end', 'status' + polarization calibration results from Delta-90° method for each wavelength. + + Each wavelength contains a list of sub-dicts. One per + successful retrieval period, with entries: + + ``eta`` : float + Depol calibration constant. + + ``eta_std`` : float + Uncertainty of Depol calibration constant. + + ``time_start``, ``time_end`` : int + Start and stop times for successful calibration. + ``method`` : str + Name of retrieval method. + + The number of element in each list depends on the number of successful retrievals. + + References + ---------- + .. [1] Freudenthaler, V. About the effects of polarising optics on lidar signals and the Delta90 calibration. + Atmos. Meas. Tech., 9, 4181–4255 (2016). + + Notes ----- - - **History** - Function is called here https://github.com/PollyNET/Pollynet_Processing_Chain/blob/5f5e4d0fd3dcebe7f87220cf802fcd6f414fe235/lib/interface/picassoProcV3.m#L548 The two most relevant functions here are https://github.com/PollyNET/Pollynet_Processing_Chain/blob/dev/lib/calibration/pollyPolCaliGHK.m which also calls https://github.com/PollyNET/Pollynet_Processing_Chain/blob/dev/lib/calibration/depolCaliGHK.m - - References - ---------- - .. [1] Freudenthaler 2016 + **History** """ @@ -158,64 +247,63 @@ def calibrateGHK(data_cube) -> dict: tel = 'FR' # currently only implemented in the far range receiver for wv in [355, 532, 1064]: - if np.any(data_cube.gf(wv, 'total', tel)) and np.any(data_cube.gf(wv, 'cross', tel)): - logging.info(f'and even a {wv} channel') + logging.info(f"Channels: {wv} total {tel} | {wv} cross {tel}") + if not np.any(data_cube.gf(wv, 'total', tel)) or not np.any(data_cube.gf(wv, 'cross', tel)): + logging.warning(f"Total or cross signal missing at {wv} {tel}. Skipping calibration for this channel.") + continue - sigBGCor_total = np.squeeze(data_cube.retrievals_highres['sigBGCor'][:, :, data_cube.gf(wv, 'total', tel)]) - bg_total = np.squeeze(data_cube.retrievals_highres['BG'][:, data_cube.gf(wv, 'total', tel)]) - sigBGCor_cross = np.squeeze(data_cube.retrievals_highres['sigBGCor'][:, :, data_cube.gf(wv, 'cross', tel)]) - bg_cross = np.squeeze(data_cube.retrievals_highres['BG'][:, data_cube.gf(wv, 'cross', tel)]) - - pol_cali[f"{wv}_{tel}"] = depol_cali_ghk( - signal_t=sigBGCor_total, bg_t=bg_total, - signal_x=sigBGCor_cross, bg_x=bg_cross, - time=data_cube.retrievals_highres['time'], - pol_cali_pang_start_time=data_cube.retrievals_highres['depol_cal_ang_p_time_start'], - pol_cali_pang_stop_time=data_cube.retrievals_highres['depol_cal_ang_p_time_end'], - pol_cali_nang_start_time=data_cube.retrievals_highres['depol_cal_ang_n_time_start'], - pol_cali_nang_stop_time=data_cube.retrievals_highres['depol_cal_ang_n_time_end'], - # K should be 0d? - K=np.squeeze(data_cube.polly_config_dict['K'][data_cube.gf(wv, 'total', tel)]), - cali_h_indx_range=[data_cube.polly_config_dict[f'depol_cal_minbin_{wv}'], - data_cube.polly_config_dict[f'depol_cal_maxbin_{wv}']], - SNRmin=data_cube.polly_config_dict[f'depol_cal_SNRmin_{wv}'], - sig_max=data_cube.polly_config_dict[f'depol_cal_sigMax_{wv}'], - rel_std_dplus=data_cube.polly_config_dict[f'rel_std_dplus_{wv}'], - rel_std_dminus=data_cube.polly_config_dict[f'rel_std_dminus_{wv}'], - segment_len=data_cube.polly_config_dict[f'depol_cal_segmentLen_{wv}'], - smooth_win=data_cube.polly_config_dict[f'depol_cal_smoothWin_{wv}'], - collect_debug=False - ) - print(pol_cali[f'{wv}_{tel}']) - logging.info(f"pol_cali_{wv} {pol_cali[f'{wv}_{tel}']}") - else: - logging.warning(f'calibrateGHK no {wv} channel') - - # TODO handling of default and database calibrations + sigBGCor_total = np.squeeze(data_cube.retrievals_highres['sigBGCor'][:, :, data_cube.gf(wv, 'total', tel)]) + bg_total = np.squeeze(data_cube.retrievals_highres['BG'][:, data_cube.gf(wv, 'total', tel)]) + sigBGCor_cross = np.squeeze(data_cube.retrievals_highres['sigBGCor'][:, :, data_cube.gf(wv, 'cross', tel)]) + bg_cross = np.squeeze(data_cube.retrievals_highres['BG'][:, data_cube.gf(wv, 'cross', tel)]) + + pol_cali[f"{wv}_{tel}"] = depol_cali_ghk( + signal_t=sigBGCor_total, bg_t=bg_total, + signal_x=sigBGCor_cross, bg_x=bg_cross, + time=data_cube.retrievals_highres['time'], + pol_cali_pang_start_time=data_cube.retrievals_highres['depol_cal_ang_p_time_start'], + pol_cali_pang_stop_time=data_cube.retrievals_highres['depol_cal_ang_p_time_end'], + pol_cali_nang_start_time=data_cube.retrievals_highres['depol_cal_ang_n_time_start'], + pol_cali_nang_stop_time=data_cube.retrievals_highres['depol_cal_ang_n_time_end'], + # K should be 0d? + K=np.squeeze(data_cube.polly_config_dict['K'][data_cube.gf(wv, 'total', tel)]), + cali_h_indx_range=[data_cube.polly_config_dict[f'depol_cal_minbin_{wv}'], + data_cube.polly_config_dict[f'depol_cal_maxbin_{wv}']], + SNRmin=data_cube.polly_config_dict[f'depol_cal_SNRmin_{wv}'], + sig_max=data_cube.polly_config_dict[f'depol_cal_sigMax_{wv}'], + rel_std_dplus=data_cube.polly_config_dict[f'rel_std_dplus_{wv}'], + rel_std_dminus=data_cube.polly_config_dict[f'rel_std_dminus_{wv}'], + segment_len=data_cube.polly_config_dict[f'depol_cal_segmentLen_{wv}'], + smooth_win=data_cube.polly_config_dict[f'depol_cal_smoothWin_{wv}'], + collect_debug=collect_debug, + flagPicassoComparison=data_cube.polly_config_dict['flagPicassoComparison'] + ) + logging.info(f"Calibration results at {wv} {tel}: {pol_cali[f'{wv}_{tel}']}") return pol_cali -def depol_cali_ghk(signal_t:np.ndarray, bg_t:np.ndarray, signal_x:np.ndarray, bg_x:np.ndarray, +def depol_cali_ghk(signal_t:np.ndarray, bg_t:np.ndarray, signal_x:np.ndarray, bg_x:np.ndarray, time:np.ndarray, pol_cali_pang_start_time:np.ndarray, pol_cali_pang_stop_time:np.ndarray, pol_cali_nang_start_time:np.ndarray, pol_cali_nang_stop_time:np.ndarray, K:float, cali_h_indx_range:list|tuple, SNRmin:list, sig_max:list, rel_std_dplus:float, rel_std_dminus:float, - segment_len:int, smooth_win:int, collect_debug:bool=False) -> dict: + segment_len:int, smooth_win:int, collect_debug:bool=False, + flagPicassoComparison:bool=False) -> list: """Polarization calibration for PollyXT lidar system. Parameters ---------- signal_t : ndarray Background-removed photon count signal at the total channel. - Shape: (n_bins, n_profiles) + Shape: (n_profiles, n_bins) bg_t : ndarray - Background at the total channel. Shape: (n_bins, n_profiles) + Background at the total channel. Shape: (n_profiles) signal_x : ndarray Background-removed photon count signal at the cross channel. - Shape: (n_bins, n_profiles) + Shape: (n_profiles, n_bins) bg_x : ndarray - Background at the cross channel. Shape: (n_bins, n_profiles) + Background at the cross channel. Shape: (n_profiles) time : ndarray Datetime array representing the measurement time of each profile. pol_cali_pang_start_time, pol_cali_pang_stop_time : ndarray @@ -236,38 +324,53 @@ def depol_cali_ghk(signal_t:np.ndarray, bg_t:np.ndarray, signal_x:np.ndarray, bg Segment length for testing the variability of calibration results. smooth_win : int Width of the sliding window for smoothing the signal. - collect_debug : bool, default=False - store and return the intermediate results + collect_debug : bool, optional + Store and return the intermediate results. Default is False. Returns ------- - pol_cali_eta : list - Eta values from polarization calibration. - pol_cali_eta_std : list - Uncertainty of eta values from calibration. - pol_cali_start_time, pol_cali_stop_time : list - Start and stop times for successful calibration. - cali_status : int - 1 if calibration is successful, 0 otherwise. - global_attri : dict, optional - Information about the depolarization calibration. + results : list of dicts + Containing successful retrieved depol calibration constant and retrieval information. + + Each dict contains the following entries: + + ``eta`` : float + Eta values from polarization calibration. + + ``eta_std`` : float + Uncertainty of eta values from calibration. + + ``time_start``, ``time_end`` : int + Start and stop times for successful calibration. + + ``method`` : str + Name of retrieval method. + + The number of element in the list depends on the number of successful retrievals. + + Notes + ----- + .. TODO:: Why is ``pol_cali_nang_start_time`` returned as ``time_start`` of the calibration + period instead of ``pol_cali_start_time``??? + + **History** + + - xxxx-xx-xx: First edition by ... + - 2026-07-30: Made output datatype consistent and changed from `nanmean` to + `nansum` in signal aggregation to be consistent with SNR requirement. + - 2026-08-06: Overhalled function structure and input output consistency. + """ - # Initialize outputs and intermediate storage - pol_cali_eta, pol_cali_eta_std = [], [] - mean_dplus, mean_dminus, std_dplus, std_dminus = [], [], [], [] - pol_cali_start_time, pol_cali_stop_time = [], [] - if collect_debug: - global_attri = defaultdict(list) # the beauty of a proper programming language + + ## Initialize output + results = [] if signal_t.size == 0 or signal_x.size == 0: - logging.warning("Warning: No data for polarization calibration.") - #return pol_cali_eta, pol_cali_eta_std, pol_cali_start_time, pol_cali_stop_time, 0, global_attri - return {'status': 0} + logging.warning("No signal for calibration.") + return results # the iteration of days can be omitted if unixtimestamps are used - print('pol_cali_nang_start_time', pol_cali_nang_start_time) - - time = np.array(time) + time = np.asarray(time) for i_depol_cal in range(len(pol_cali_nang_start_time)): indx_45p = np.where( (time >= pol_cali_pang_start_time[i_depol_cal]) & @@ -276,41 +379,51 @@ def depol_cali_ghk(signal_t:np.ndarray, bg_t:np.ndarray, signal_x:np.ndarray, bg indx_45m = np.where( (time >= pol_cali_nang_start_time[i_depol_cal]) & (time <= pol_cali_nang_stop_time[i_depol_cal]))[0] + if len(indx_45p) < 4 or len(indx_45m) < 4: - logging.warning(f'calibrateGHK array to short {len(indx_45p)}{len(indx_45m)} in period {i_depol_cal}') - break - this_cali_start_time = min(pol_cali_pang_start_time[i_depol_cal], - pol_cali_nang_start_time[i_depol_cal]) - this_cali_stop_time = max(pol_cali_pang_stop_time[i_depol_cal], - pol_cali_nang_stop_time[i_depol_cal]) - # Exclude the first and last profiles + logging.warning(f"Not enough calibration profiles in clibration period {i_depol_cal}. Skipping this period.") + continue + + ## Get start and end time of calibration period + pol_cali_start_time = min( + pol_cali_pang_start_time[i_depol_cal], + pol_cali_nang_start_time[i_depol_cal] + ) + pol_cali_stop_time = max( + pol_cali_pang_stop_time[i_depol_cal], + pol_cali_nang_stop_time[i_depol_cal] + ) + + ## Exclude the first and last profiles indx_45m = indx_45m[1:-1] indx_45p = indx_45p[1:-1] - # matlab -> python swap from signal_t[:, indx_45p] to signal_t[indx_45p,:] - # to be a profile - sig_t_p = np.nanmean(signal_t[indx_45p, :], axis=0) - bg_t_p = np.nanmean(bg_t[indx_45p], axis=0) + ## Calculating SNR (only sum should be used when aggregating signals for SNR calculations!) + func = np.nansum + if flagPicassoComparison: + func = np.nanmean + + sig_t_p = func(signal_t[indx_45p, :], axis=0) + bg_t_p = func(bg_t[indx_45p], axis=0) snr_t_p = calc_snr(sig_t_p, bg_t_p) indx_bad_t_p = (snr_t_p <= SNRmin[0]) | (sig_t_p >= sig_max[0]) - sig_t_m = np.nanmean(signal_t[indx_45m, :], axis=0) - bg_t_m = np.nanmean(bg_t[indx_45m], axis=0) + sig_t_m = func(signal_t[indx_45m, :], axis=0) + bg_t_m = func(bg_t[indx_45m], axis=0) snr_t_m = calc_snr(sig_t_m, bg_t_m) indx_bad_t_m = (snr_t_m <= SNRmin[1]) | (sig_t_m >= sig_max[1]) - sig_x_p = np.nanmean(signal_x[indx_45p, :], axis=0) - bg_x_p = np.nanmean(bg_x[indx_45p], axis=0) + sig_x_p = func(signal_x[indx_45p, :], axis=0) + bg_x_p = func(bg_x[indx_45p], axis=0) snr_x_p = calc_snr(sig_x_p, bg_x_p) indx_bad_x_p = (snr_x_p <= SNRmin[2]) | (sig_x_p >= sig_max[2]) - sig_x_m = np.nanmean(signal_x[indx_45m, :], axis=0) - bg_x_m = np.nanmean(bg_x[indx_45m], axis=0) + sig_x_m = func(signal_x[indx_45m, :], axis=0) + bg_x_m = func(bg_x[indx_45m], axis=0) snr_x_m = calc_snr(sig_x_m, bg_x_m) indx_bad_x_m = (snr_x_m <= SNRmin[3]) | (sig_x_m >= sig_max[3]) - # Calculate dplus and dminus - #print('smooth_win', smooth_win) + ## Calculate dplus and dminus dplus = smooth_signal(sig_x_p, smooth_win) / smooth_signal(sig_t_p, smooth_win) dminus = smooth_signal(sig_x_m, smooth_win) / smooth_signal(sig_t_m, smooth_win) dplus = np.where(np.isfinite(dplus), dplus, np.nan) @@ -318,122 +431,145 @@ def depol_cali_ghk(signal_t:np.ndarray, bg_t:np.ndarray, signal_x:np.ndarray, bg dplus[indx_bad_t_p | indx_bad_x_p] = np.nan dminus[indx_bad_t_m | indx_bad_x_m] = np.nan - # Subset the calibration range - dplus = dplus[cali_h_indx_range[0]:cali_h_indx_range[1]] - dminus = dminus[cali_h_indx_range[0]:cali_h_indx_range[1]] + ## Subset the calibration range + dplus = dplus[cali_h_indx_range[0]:cali_h_indx_range[1]+1] + dminus = dminus[cali_h_indx_range[0]:cali_h_indx_range[1]+1] if np.all(np.isnan(dplus)) or np.all(np.isnan(dminus)): - logging.warning(f'calibrateGHK all values in dplus or dminus masked in period {i_depol_cal}') - print(f'calibrateGHK all values in dplus or dminus masked in period {i_depol_cal}, len(dplus) {len(dplus)}') - print(' snr_t_p ', np.sum((snr_t_p <= SNRmin[0])[cali_h_indx_range[0]:cali_h_indx_range[1]])) - print(' sig_t_p ', np.sum((sig_t_p >= sig_max[0])[cali_h_indx_range[0]:cali_h_indx_range[1]])) - print('> indx_bad_t_p in height interval', np.sum(indx_bad_t_p[cali_h_indx_range[0]:cali_h_indx_range[1]])) - print(' snr_t_m ', np.sum((snr_t_m <= SNRmin[1])[cali_h_indx_range[0]:cali_h_indx_range[1]])) - print(' sig_t_m ', np.sum((sig_t_m >= sig_max[1])[cali_h_indx_range[0]:cali_h_indx_range[1]])) - print('> indx_bad_t_m in height interval', np.sum(indx_bad_t_m[cali_h_indx_range[0]:cali_h_indx_range[1]])) - print(' snr_x_p ', np.sum((snr_x_p <= SNRmin[2])[cali_h_indx_range[0]:cali_h_indx_range[1]])) - print(' sig_x_p ', np.sum((sig_x_p >= sig_max[2])[cali_h_indx_range[0]:cali_h_indx_range[1]])) - print('> indx_bad_x_p in height interval', np.sum(indx_bad_x_p[cali_h_indx_range[0]:cali_h_indx_range[1]])) - print(' snr_x_m ', np.sum((snr_x_m <= SNRmin[3])[cali_h_indx_range[0]:cali_h_indx_range[1]])) - print(' sig_x_m ', np.sum((sig_x_m >= sig_max[3])[cali_h_indx_range[0]:cali_h_indx_range[1]])) - print('> indx_bad_x_m in height interval', np.sum(indx_bad_x_m[cali_h_indx_range[0]:cali_h_indx_range[1]])) + logging.warning(f"No valid plus or minus 45° calibration found in calibration period {i_depol_cal}. Skipping this period.") + + ## Debug info + logging.debug(f"CalibrateGHK all values in dplus or dminus masked in period {i_depol_cal}, len(dplus) {len(dplus)}") + logging.debug(f" snr_t_p {np.sum((snr_t_p <= SNRmin[0])[cali_h_indx_range[0]:cali_h_indx_range[1]+1])}") + logging.debug(f" sig_t_p {np.sum((sig_t_p >= sig_max[0])[cali_h_indx_range[0]:cali_h_indx_range[1]+1])}") + logging.debug(f"> indx_bad_t_p in height interval {np.sum(indx_bad_t_p[cali_h_indx_range[0]:cali_h_indx_range[1]+1])}") + logging.debug(f" snr_t_m {np.sum((snr_t_m <= SNRmin[1])[cali_h_indx_range[0]:cali_h_indx_range[1]+1])}") + logging.debug(f" sig_t_m {np.sum((sig_t_m >= sig_max[1])[cali_h_indx_range[0]:cali_h_indx_range[1]+1])}") + logging.debug(f"> indx_bad_t_m in height interval {np.sum(indx_bad_t_m[cali_h_indx_range[0]:cali_h_indx_range[1]+1])}") + logging.debug(f" snr_x_p {np.sum((snr_x_p <= SNRmin[2])[cali_h_indx_range[0]:cali_h_indx_range[1]+1])}") + logging.debug(f" sig_x_p {np.sum((sig_x_p >= sig_max[2])[cali_h_indx_range[0]:cali_h_indx_range[1]+1])}") + logging.debug(f"> indx_bad_x_p in height interval {np.sum(indx_bad_x_p[cali_h_indx_range[0]:cali_h_indx_range[1]+1])}") + logging.debug(f" snr_x_m {np.sum((snr_x_m <= SNRmin[3])[cali_h_indx_range[0]:cali_h_indx_range[1]+1])}") + logging.debug(f" sig_x_m {np.sum((sig_x_m >= sig_max[3])[cali_h_indx_range[0]:cali_h_indx_range[1]+1])}") + logging.debug(f"> indx_bad_x_m in height interval {np.sum(indx_bad_x_m[cali_h_indx_range[0]:cali_h_indx_range[1]+1])}") continue - # Analyze segments for stability - #print('before analyze segments', len(dplus), segment_len) - seg = analyze_segments(dplus, dminus, segment_len, rel_std_dplus, rel_std_dminus) - if seg.shape[0] == 0: - logging.warning(f'calibrateGHK no stable segment found in period {i_depol_cal}') + ## Analyze segments for stability + mean_seg_dpluses, std_seg_dpluses, mean_seg_dminuses, std_seg_dminuses = analyze_segments( + dplus=dplus, + dminus=dminus, + segment_len=segment_len, + rel_std_dplus=rel_std_dplus, + rel_std_dminus=rel_std_dminus + ) + if mean_seg_dpluses.size == 0: + logging.warning(f"No stable calibration segment found in calibration period {i_depol_cal}. Skipping this period.") continue - # translate manually - # min(sqrt((std_dplus_tmp./mean_dplus_tmp).^2 + (std_dminus_tmp./mean_dminus_tmp).^2)); - indx_best_seg = np.argmin(np.sqrt((seg[:, 1]/seg[:, 0])**2 + (seg[:, 3]/seg[:, 2])**2)) - # the best segment searching was flawed by the AI translate - best_segment = seg[indx_best_seg] - mean_dplus.append(best_segment[0]) - std_dplus.append(best_segment[1]) - mean_dminus.append(best_segment[2]) - std_dminus.append(best_segment[3]) - pol_cali_start_time.append(this_cali_start_time) - pol_cali_stop_time.append(this_cali_stop_time) - + ## Find optimal segment + best_seg_idx = np.argmin( + np.sqrt((std_seg_dpluses/mean_seg_dpluses)**2 + (std_seg_dminuses/mean_seg_dminuses)**2) + ) + mean_dplus = mean_seg_dminuses[best_seg_idx] + std_dplus = std_seg_dpluses[best_seg_idx] + mean_dminus = mean_seg_dminuses[best_seg_idx] + std_dminus = std_seg_dminuses[best_seg_idx] + + ## Polarization calibration constant + pol_cali_eta = float(1 / K * np.sqrt(mean_dplus * mean_dminus)) + pol_cali_eta_std = float(0.5 * (mean_dplus * std_dminus + mean_dminus * std_dplus) / np.sqrt(mean_dplus * mean_dminus)) + + ## save calibration result + results.append({ + 'eta': pol_cali_eta, 'eta_std': pol_cali_eta_std, + # 'time_start': pol_cali_start_time, + 'time_start': pol_cali_nang_start_time[i_depol_cal], # TODO: Why are we returning pol_cali_nang_start_time insted of pol_cali_start_time? + 'time_end': pol_cali_stop_time, + 'method': 'D90' + }) + + ## collect debug info if collect_debug: - global_attri['sig_t_p'].append(sig_t_p) - global_attri['sig_t_m'].append(sig_t_m) - global_attri['sig_x_p'].append(sig_x_p) - global_attri['sig_x_m'].append(sig_x_m) - global_attri['cali_h_indx_range'].append(cali_h_indx_range) - global_attri['indx_45p'].append(indx_45p) - global_attri['indx_45m'].append(indx_45m) - global_attri['dplus'].append(dplus) - global_attri['dminus'].append(dminus) - global_attri['segment_len'].append(segment_len) - global_attri['indx_best_seg'].append(indx_best_seg) - global_attri['segment_results'].append(seg) - global_attri['K'].append(K) - global_attri['cali_time'].append(np.mean([this_cali_start_time, this_cali_stop_time])) - - if not mean_dplus or not mean_dminus: - logging.warning("Plus or minus 45° calibration is missing.") - #return pol_cali_eta, pol_cali_eta_std, pol_cali_start_time, pol_cali_stop_time, 0, global_attri - return {'status': 0} - - pol_cali_eta = [float(1 / K * np.sqrt(dp * dm)) for dp, dm in zip(mean_dplus, mean_dminus)] - pol_cali_eta_std = [float(0.5 * (dp * std_dm + dm * std_dp) / np.sqrt(dp * dm)) for - dp, std_dp, dm, std_dm in zip(mean_dplus, std_dplus, mean_dminus, std_dminus)] - - results = [ - {'eta': e[0], 'eta_std': e[1], 'time_start': e[2], 'time_end': e[3], 'status': 1} - for e in zip(pol_cali_eta, pol_cali_eta_std, pol_cali_nang_start_time, pol_cali_stop_time)] + results[-1]['sig_t_p'] = sig_t_p + results[-1]['sig_t_m'] = sig_t_m + results[-1]['sig_x_p'] = sig_x_p + results[-1]['sig_x_m'] = sig_x_m + results[-1]['cali_h_indx_range'] = cali_h_indx_range + results[-1]['indx_45p'] = indx_45p + results[-1]['indx_45m'] = indx_45m + results[-1]['dplus'] = dplus + results[-1]['dminus'] = dminus + results[-1]['segment_len'] = segment_len + results[-1]['indx_best_seg'] = best_seg_idx + results[-1]['mean_dplus_seg'] = mean_seg_dpluses + results[-1]['std_dplus_seg'] = std_seg_dpluses + results[-1]['mean_dminus_seg'] = mean_seg_dminuses + results[-1]['std_dminus_seg'] = std_seg_dminuses + results[-1]['K'] = K + results[-1]['cali_time'] = np.mean([ + pol_cali_start_time[i_depol_cal], + pol_cali_stop_time[i_depol_cal]]) - if collect_debug: - results['global_attri'] = dict(global_attri) return results def analyze_segments(dplus:np.ndarray, dminus:np.ndarray, segment_len:int, - rel_std_dplus:float, rel_std_dminus:float) -> np.ndarray: - """... + rel_std_dplus:float, rel_std_dminus:float) -> tuple: + """Analyze calibration segment. Parameters ---------- dplus : ndarray - ... + Plus 45° calibration... dminus : ndarray - ... + Minus 45° calibration... segment_len : int Segment length for testing the variability of calibration results. rel_std_dplus, rel_std_dminus : float - Maximum relative uncertainty of dplus and dminus allowed. + Maximum relative uncertainty of `dplus` and `dminus` allowed. Returns ------- - ndarray - ... - - Notes - ----- - .. TODO:: Finish docstring. + mean_dpluses : ndarray + Mean plus 45° calibration per segment. + std_dpluses : ndarray + Standard deviation of plus 45° calibration per segment. + mean_dminuses : ndarray + Mean minus 45° calibration per segment. + std_dminuses : ndarray + Standard deviation of minus 45° calibration per segment. """ - results = [] + mean_dpluses = [] + std_dpluses = [] + mean_dminuses = [] + std_dminuses = [] + for i in range(len(dplus) - segment_len): - #print(i, i+segment_len) seg_dplus = dplus[i:i + segment_len] seg_dminus = dminus[i:i + segment_len] + if np.sum(~np.isnan(seg_dplus)) <= segment_len / 4 or np.sum(~np.isnan(seg_dminus)) <= segment_len / 4: continue + mean_dp = np.nanmean(seg_dplus) std_dp = np.nanstd(seg_dplus) mean_dm = np.nanmean(seg_dminus) std_dm = np.nanstd(seg_dminus) - #print('mean_dp', mean_dp, 'std_dp', std_dp, '-> ', std_dp / mean_dp, rel_std_dplus) - #print('mean_dm', mean_dm, 'std_dm', std_dm, '-> ', std_dm / mean_dm, rel_std_dminus) - + if std_dp / mean_dp <= rel_std_dplus and std_dm / mean_dm <= rel_std_dminus: - results.append([mean_dp, std_dp, mean_dm, std_dm]) - return np.array(results) + mean_dpluses.append(mean_dp) + std_dpluses.append(std_dp) + mean_dminuses.append(mean_dm) + std_dminuses.append(std_dm) + + # Convert to ndarray + mean_dpluses = np.asarray(mean_dpluses) + std_dpluses = np.asarray(std_dpluses) + mean_dminuses = np.asarray(mean_dminuses) + std_dminuses = np.asarray(std_dminuses) + + return mean_dpluses, std_dpluses, mean_dminuses, std_dminuses """ [data.polCaliEta532, data.polCaliEtaStd532, data.polCaliTime, data.polCali532Attri] = diff --git a/ppcpy/calibration/select.py b/ppcpy/calibration/select.py index c513d1d..4ae5e03 100644 --- a/ppcpy/calibration/select.py +++ b/ppcpy/calibration/select.py @@ -5,7 +5,7 @@ import matplotlib -def single_best(d:dict, name_val:str, name_min:str, relative:bool=False) -> dict: +def single_best(d:dict, name_val:str, name_min:str, name_method:str, relative:bool=False) -> dict: """Select the best calibration constant @@ -36,39 +36,53 @@ def single_best(d:dict, name_val:str, name_min:str, relative:bool=False) -> dict **History** - 2026-02-16: Added additional checks to hinder negative LCs to be chosen. - - 2026-03-27: generalized to also hold for depolarization calibration + - 2026-03-27: generalized to also hold for depolarization calibration. + """ best = {} for k, l in d.items(): - val = np.array([e[name_val] for e in l if e[name_val] >= 0]) - min = np.array([e[name_min] for e in l if e[name_val] >= 0]) + val = np.array([e.get(name_val) for e in l if e.get(name_val, np.nan) >= 0]) + min = np.array([e.get(name_min) for e in l if e.get(name_val, np.nan) >= 0]) + method = np.array([e.get(name_method, 'unknown') for e in l if e.get(name_val, np.nan) >= 0]) + + if len(val) == 0 or len(min) == 0: + continue if relative: - best[k] = val[np.argmin(min / val)] + idx = np.argmin(min / val) else: - best[k] = val[np.argmin(min)] + idx = np.argmin(min) + + best[k] = {name_val:val[idx], name_min:min[idx], name_method:method[idx]} return best -def plot_cals(d, param, used=None): - """plot the calibration constants +def plot_cals(d:dict, param:str, used:dict=None): + """Plot the calibration constants. + + Produces a scatter plot of the stored calibration constants (CC). + CCs retrieved form the measurement are marked as filled circles, + while those loaded from the database are marked as hollow circles. + The chosen optimal CC and its default value are marked by horizontal + gray dotted and dashed lines respectively. Parameters ---------- d : dict - the dict as in data_cube + Dict storing all CCs. ``pol_cali`` or ``LC``. param : str - the parameter to extract - used : dict - the LCused or etaused (will produce a dashed line in the plot) - + Name of the parameter to extract. + used : dict, optional + Dict storing the used CCs. LCused or etaused. + Will produce a horizontal dashed line in the plot if added. + Default is None. Examples -------- - >>> plot_cals(data_cube.pol_cali, 'eta', used=data_cube.etaused) + >>> plot_cals(data_cube.LC, 'LC', used=data_cube.LCused) """ @@ -81,22 +95,29 @@ def plot_cals(d, param, used=None): for c in channels: guess_yscale = [] - fig, ax = plt.subplots(figsize=(8,4)) + fig, ax = plt.subplots(figsize=(8, 4)) + if used and c in used.keys(): + ax.axhline(used[c][param], color='dimgrey', ls=':', label="used") + for k, v in d.items(): if not c in v: continue + if "default" in k: + ax.axhline(d[k][c][0][param], color='dimgray', alpha=0.6, ls='--', label="default") + continue + if 'db' in k: marker = 'o' fillstyle = 'none' else: marker = '.' - fillstyle = 'full' + fillstyle = 'full' time_mean = np.array([np.mean([e['time_start'], e['time_end']]) for e in v[c]]).astype('datetime64[s]') eta = [e[param] for e in v[c]] - if used and c in used.keys(): - ax.axhline(used[c], color='dimgrey', ls=':') + # if used and c in used.keys(): + # ax.axhline(used[c][param], color='dimgrey', ls=':', label="used") ax.plot(time_mean, eta, marker, label=k, fillstyle=fillstyle) guess_yscale.append(np.min(eta)*0.95) @@ -104,7 +125,7 @@ def plot_cals(d, param, used=None): guess_yscale.append(np.mean(eta)*1.2) guess_yscale.append(np.mean(eta)*0.8) - ax.set_ylim(np.max([0,np.min(guess_yscale)]), np.max(guess_yscale)) + ax.set_ylim(np.max([0, np.min(guess_yscale)]), np.max(guess_yscale)) ax.set_ylabel(param) ax.legend() ax.set_title(c) diff --git a/ppcpy/config/json2nc-mapper_NR_att_bsc.json b/ppcpy/config/json2nc-mapper_NR_att_bsc.json index e24db75..18c2b5f 100644 --- a/ppcpy/config/json2nc-mapper_NR_att_bsc.json +++ b/ppcpy/config/json2nc-mapper_NR_att_bsc.json @@ -105,7 +105,7 @@ "plot_range": "0., 1.5e-05", "plot_scale": "linear", "retrieving_info": { - "Lidar_calibration_constant_used": {"value": "__LCused[355_total_NR]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[355_total_NR][LC]", "unit": "", "method": "__LCused[355_total_NR][method]"} }, "comment": "This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter." } @@ -126,7 +126,7 @@ "plot_range": "0., 5.e-06", "plot_scale": "linear", "retrieving_info": { - "Lidar_calibration_constant_used": {"value": "__LCused[532_total_NR]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[532_total_NR][LC]", "unit": "", "method": "__LCused[532_total_NR][method]"} }, "comment": "This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter." } diff --git a/ppcpy/config/json2nc-mapper_OC_att_bsc.json b/ppcpy/config/json2nc-mapper_OC_att_bsc.json index db3dfc9..aff27a7 100644 --- a/ppcpy/config/json2nc-mapper_OC_att_bsc.json +++ b/ppcpy/config/json2nc-mapper_OC_att_bsc.json @@ -105,7 +105,7 @@ "plot_range": "0.0, 2.0e-05", "plot_scale": "linear", "retrieving_info": { - "Lidar_calibration_constant_used": {"value": "__LCused[355_total_OC]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[355_total_FR][LC]", "unit": "", "method": "__LCused[355_total_FR][method]"} }, "comment": "This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter." } @@ -126,7 +126,7 @@ "plot_range": "0.0, 2.0e-05", "plot_scale": "linear", "retrieving_info": { - "Lidar_calibration_constant_used": {"value": "__LCused[532_total_OC]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[532_total_FR][LC]", "unit": "", "method": "__LCused[532_total_FR][method]"} }, "comment": "This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter." } @@ -147,7 +147,7 @@ "plot_range": "0.0, 1.5.e-05", "plot_scale": "linear", "retrieving_info": { - "Lidar_calibration_constant_used": {"value": "__LCused[1064_total_OC]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[1064_total_FR][LC]", "unit": "", "method": "__LCused[1064_total_FR][method]"} }, "comment": "This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter." } diff --git a/ppcpy/config/json2nc-mapper_att_bsc.json b/ppcpy/config/json2nc-mapper_att_bsc.json index d156f95..67d6069 100644 --- a/ppcpy/config/json2nc-mapper_att_bsc.json +++ b/ppcpy/config/json2nc-mapper_att_bsc.json @@ -105,7 +105,7 @@ "plot_range": "0.0, 2.0e-05", "plot_scale": "linear", "retrieving_info": { - "Lidar_calibration_constant_used": {"value": "__LCused[355_total_FR]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[355_total_FR][LC]", "unit": "", "method": "__LCused[355_total_FR][method]"} }, "comment": "This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter." } @@ -126,7 +126,7 @@ "plot_range": "0.0, 2.0e-05", "plot_scale": "linear", "retrieving_info": { - "Lidar_calibration_constant_used": {"value": "__LCused[532_total_FR]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[532_total_FR][LC]", "unit": "", "method": "__LCused[532_total_FR][method]"} }, "comment": "This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter." } @@ -147,7 +147,7 @@ "plot_range": "0.0, 1.5e-05", "plot_scale": "linear", "retrieving_info": { - "Lidar_calibration_constant_used": {"value": "__LCused[1064_total_FR]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[1064_total_FR][LC]", "unit": "", "method": "__LCused[1064_total_FR][method]"} }, "comment": "This parameter is calculated with taking into account of the effects of lidar constants. Therefore, it reflects the concentration of aerosol and molecule backscatter." } diff --git a/ppcpy/config/json2nc-mapper_profiles.json b/ppcpy/config/json2nc-mapper_profiles.json index 3bb2144..72a75c6 100644 --- a/ppcpy/config/json2nc-mapper_profiles.json +++ b/ppcpy/config/json2nc-mapper_profiles.json @@ -330,7 +330,7 @@ "source": "__device", "retrieving_info": { "Smoothing window": {"value": "__polly_config_dict[smoothWin_klett_355]", "unit": "[m]"}, - "eta":{"value": "__etaused[355_FR]", "unit": ""} + "eta":{"value": "__etaused[355_FR][eta]", "unit": "", "method": "__etaused[355_FR][method]"} }, "comment": "Depolarization channel was calibrated with +- 45 \\degree method. You can find more information in Freudenthaler, V., et al. (2009). \"Depolarization ratio profiling at several wavelengths in pure Saharan dust during SAMUM 2006.\" Tellus B 61(1): 165-179." } @@ -365,7 +365,7 @@ "source": "__device", "retrieving_info": { "Smoothing window": {"value": "__polly_config_dict[smoothWin_klett_532]", "unit": "[m]"}, - "eta":{"value": "__etaused[532_FR]", "unit": ""} + "eta":{"value": "__etaused[532_FR][eta]", "unit": "", "method": "__etaused[532_FR][method]"} }, "comment": "Depolarization channel was calibrated with +- 45 \\degree method. You can find more information in Freudenthaler, V., et al. (2009). \"Depolarization ratio profiling at several wavelengths in pure Saharan dust during SAMUM 2006.\" Tellus B 61(1): 165-179." } @@ -400,7 +400,7 @@ "source": "__device", "retrieving_info": { "Smoothing window": {"value": "__polly_config_dict[smoothWin_klett_1064]", "unit": "[m]"}, - "eta":{"value": "__etaused[1064_FR]", "unit": ""} + "eta":{"value": "__etaused[1064_FR][eta]", "unit": "", "method": "__etaused[1064_FR][method]"} }, "comment": "Depolarization channel was calibrated with +- 45 \\degree method. You can find more information in Freudenthaler, V., et al. (2009). \"Depolarization ratio profiling at several wavelengths in pure Saharan dust during SAMUM 2006.\" Tellus B 61(1): 165-179." } @@ -435,7 +435,7 @@ "source": "__device", "retrieving_info": { "Smoothing window": {"value": "__polly_config_dict[smoothWin_raman_355]", "unit": "[m]"}, - "eta":{"value": "__etaused[355_FR]", "unit": ""} + "eta":{"value": "__etaused[355_FR][eta]", "unit": "", "method": "__etaused[355_FR][method]"} }, "comment": "Depolarization channel was calibrated with +- 45 \\degree method. You can find more information in Freudenthaler, V., et al. (2009). \"Depolarization ratio profiling at several wavelengths in pure Saharan dust during SAMUM 2006.\" Tellus B 61(1): 165-179." } @@ -470,7 +470,7 @@ "source": "__device", "retrieving_info": { "Smoothing window": {"value": "__polly_config_dict[smoothWin_raman_532]", "unit": "[m]"}, - "eta":{"value": "__etaused[532_FR]", "unit": ""} + "eta":{"value": "__etaused[532_FR][eta]", "unit": "", "method": "__etaused[532_FR][method]"} }, "comment": "Depolarization channel was calibrated with +- 45 \\degree method. You can find more information in Freudenthaler, V., et al. (2009). \"Depolarization ratio profiling at several wavelengths in pure Saharan dust during SAMUM 2006.\" Tellus B 61(1): 165-179." } @@ -505,7 +505,7 @@ "source": "__device", "retrieving_info": { "Smoothing window": {"value": "__polly_config_dict[smoothWin_raman_1064]", "unit": "[m]"}, - "eta":{"value": "__etaused[1064_FR]", "unit": ""} + "eta":{"value": "__etaused[1064_FR][eta]", "unit": "", "method": "__etaused[1064_FR][method]"} }, "comment": "Depolarization channel was calibrated with +- 45 \\degree method. You can find more information in Freudenthaler, V., et al. (2009). \"Depolarization ratio profiling at several wavelengths in pure Saharan dust during SAMUM 2006.\" Tellus B 61(1): 165-179." } @@ -540,7 +540,7 @@ "source": "__device", "retrieving_info": { "Smoothing window": {"value": "__polly_config_dict[smoothWin_klett_355]", "unit": "[m]"}, - "eta":{"value": "__etaused[355_FR]", "unit": ""} + "eta":{"value": "__etaused[355_FR][eta]", "unit": "", "method": "__etaused[355_FR][method]"} }, "comment": "The aerosol backscatter profile was retrieved by Raman method. The uncertainty of particle depolarization ratio will be very large at aerosol-free altitude. Please take care!" } @@ -575,7 +575,7 @@ "source": "__device", "retrieving_info": { "Smoothing window": {"value": "__polly_config_dict[smoothWin_klett_532]", "unit": "[m]"}, - "eta":{"value": "__etaused[532_FR]", "unit": ""} + "eta":{"value": "__etaused[532_FR][eta]", "unit": "", "method": "__etaused[532_FR][method]"} }, "comment": "The aerosol backscatter profile was retrieved by Raman method. The uncertainty of particle depolarization ratio will be very large at aerosol-free altitude. Please take care!" } @@ -610,7 +610,7 @@ "source": "__device", "retrieving_info": { "Smoothing window": {"value": "__polly_config_dict[smoothWin_klett_1064]", "unit": "[m]"}, - "eta":{"value": "__etaused[355_FR]", "unit": ""} + "eta":{"value": "__etaused[355_FR][eta]", "unit": "", "method": "__etaused[1064_FR][method]"} }, "comment": "The aerosol backscatter profile was retrieved by Raman method. The uncertainty of particle depolarization ratio will be very large at aerosol-free altitude. Please take care!" } @@ -645,7 +645,7 @@ "source": "__device", "retrieving_info": { "Smoothing window": {"value": "__polly_config_dict[smoothWin_raman_355]", "unit": "[m]"}, - "eta":{"value": "__etaused[355_FR]", "unit": ""} + "eta":{"value": "__etaused[355_FR][eta]", "unit": "", "method": "__etaused[355_FR][method]"} }, "comment": "The aerosol backscatter profile was retrieved by Raman method. The uncertainty of particle depolarization ratio will be very large at aerosol-free altitude. Please take care!" } @@ -680,7 +680,7 @@ "source": "__device", "retrieving_info": { "Smoothing window": {"value": "__polly_config_dict[smoothWin_raman_532]", "unit": "[m]"}, - "eta":{"value": "__etaused[532_FR]", "unit": ""} + "eta":{"value": "__etaused[532_FR][eta]", "unit": "", "method": "__etaused[532_FR][method]"} }, "comment": "The aerosol backscatter profile was retrieved by Raman method. The uncertainty of particle depolarization ratio will be very large at aerosol-free altitude. Please take care!" } @@ -715,7 +715,7 @@ "source": "__device", "retrieving_info": { "Smoothing window": {"value": "__polly_config_dict[smoothWin_raman_1064]", "unit": "[m]"}, - "eta":{"value": "__etaused[1064_FR]", "unit": ""} + "eta":{"value": "__etaused[1064_FR][eta]", "unit": "", "method": "__etaused[1064_FR][method]"} }, "comment": "The aerosol backscatter profile was retrieved by Raman method. The uncertainty of particle depolarization ratio will be very large at aerosol-free altitude. Please take care!" } diff --git a/ppcpy/config/json2nc-mapper_quasi_results.json b/ppcpy/config/json2nc-mapper_quasi_results.json index 283f7c1..5ac45d4 100644 --- a/ppcpy/config/json2nc-mapper_quasi_results.json +++ b/ppcpy/config/json2nc-mapper_quasi_results.json @@ -103,7 +103,7 @@ "standard_name": "quasi_bsc_355", "retrieving_info": { "Fixed lidar ratio": {"value": "__polly_config_dict[LR355]", "unit": "[Sr]"}, - "Lidar_calibration_constant_used": {"value": "__LCused[355_total_FR]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[355_total_FR][LC]", "unit": "", "method": "__LCused[355_total_FR][method]"} }, "comment": "This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!" } @@ -122,7 +122,7 @@ "standard_name": "quasi_bsc_532", "retrieving_info": { "Fixed lidar ratio": {"value": "__polly_config_dict[LR532]", "unit": "[Sr]"}, - "Lidar_calibration_constant_used": {"value": "__LCused[532_total_FR]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[532_total_FR][LC]", "unit": "", "method": "__LCused[532_total_FR][method]"} }, "comment": "This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!" } @@ -141,7 +141,7 @@ "standard_name": "quasi_bsc_1064", "retrieving_info": { "Fixed lidar ratio": {"value": "__polly_config_dict[LR1064]", "unit": "[Sr]"}, - "Lidar_calibration_constant_used": {"value": "__LCused[1064_total_FR]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[1064_total_FR][LC]", "unit": "", "method": "__LCused[1064_total_FR][method]"} }, "comment": "This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!" } @@ -214,7 +214,7 @@ "standard_name": "quasi_voldepol_532", "retrieving_info": { "Fixed lidar ratio": {"value": "__polly_config_dict[LR532]", "unit": "[Sr]"}, - "eta": {"value": "__etaused[532_FR]", "unit": ""} + "eta": {"value": "__etaused[532_FR][eta]", "unit": "", "method": "__etaused[532_FR][method]"} }, "comment": "This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin." } @@ -233,7 +233,7 @@ "standard_name": "quasi_pardepol_532", "retrieving_info": { "Fixed lidar ratio": {"value": "__polly_config_dict[LR532]", "unit": "[Sr]"}, - "eta": {"value": "__etaused[532_FR]", "unit": ""} + "eta": {"value": "__etaused[532_FR][eta]", "unit": "", "method": "__etaused[532_FR][method]"} }, "comment": "This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin." } diff --git a/ppcpy/config/json2nc-mapper_quasi_results_V2.json b/ppcpy/config/json2nc-mapper_quasi_results_V2.json index 283f7c1..7632688 100644 --- a/ppcpy/config/json2nc-mapper_quasi_results_V2.json +++ b/ppcpy/config/json2nc-mapper_quasi_results_V2.json @@ -103,7 +103,7 @@ "standard_name": "quasi_bsc_355", "retrieving_info": { "Fixed lidar ratio": {"value": "__polly_config_dict[LR355]", "unit": "[Sr]"}, - "Lidar_calibration_constant_used": {"value": "__LCused[355_total_FR]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[355_total_FR][LC]", "unit": "", "method": "__LCused[355_total_FR][method]"} }, "comment": "This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!" } @@ -122,7 +122,7 @@ "standard_name": "quasi_bsc_532", "retrieving_info": { "Fixed lidar ratio": {"value": "__polly_config_dict[LR532]", "unit": "[Sr]"}, - "Lidar_calibration_constant_used": {"value": "__LCused[532_total_FR]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[532_total_FR][LC]", "unit": "", "method": "__LCused[632_total_FR][method]"} }, "comment": "This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!" } @@ -141,7 +141,7 @@ "standard_name": "quasi_bsc_1064", "retrieving_info": { "Fixed lidar ratio": {"value": "__polly_config_dict[LR1064]", "unit": "[Sr]"}, - "Lidar_calibration_constant_used": {"value": "__LCused[1064_total_FR]", "unit": ""} + "Lidar_calibration_constant_used": {"value": "__LCused[1064_total_FR][LC]", "unit": "", "method": "__LCused[1064_total_FR][method]"} }, "comment": "This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!" } @@ -214,7 +214,7 @@ "standard_name": "quasi_voldepol_532", "retrieving_info": { "Fixed lidar ratio": {"value": "__polly_config_dict[LR532]", "unit": "[Sr]"}, - "eta": {"value": "__etaused[532_FR]", "unit": ""} + "eta": {"value": "__etaused[532_FR][eta]", "unit": "", "method": "__etaused[532_FR][method]"} }, "comment": "This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin." } @@ -233,7 +233,7 @@ "standard_name": "quasi_pardepol_532", "retrieving_info": { "Fixed lidar ratio": {"value": "__polly_config_dict[LR532]", "unit": "[Sr]"}, - "eta": {"value": "__etaused[532_FR]", "unit": ""} + "eta": {"value": "__etaused[532_FR][eta]", "unit": "", "method": "__etaused[532_FR][method]"} }, "comment": "This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin." } diff --git a/ppcpy/config/json2nc-mapper_vol_depol.json b/ppcpy/config/json2nc-mapper_vol_depol.json index daed019..29457c2 100644 --- a/ppcpy/config/json2nc-mapper_vol_depol.json +++ b/ppcpy/config/json2nc-mapper_vol_depol.json @@ -105,7 +105,7 @@ "plot_range": "0.0, 0.3", "plot_scale": "linear", "retrieving_info": { - "eta": {"value": "__etaused[355_FR]", "unit": ""} + "eta": {"value": "__etaused[355_FR][eta]", "unit": "", "method": "__etaused[355_FR][method]"} }, "comment": "The depolarization ratio was calibrated with \\Delta 90\\circ method." } @@ -126,7 +126,7 @@ "plot_range": "0.0, 0.3", "plot_scale": "linear", "retrieving_info": { - "eta": {"value": "__etaused[532_FR]", "unit": ""} + "eta": {"value": "__etaused[532_FR][eta]", "unit": "", "method": "__etaused[532_FR][method]"} }, "comment": "The depolarization ratio was calibrated with \\Delta 90\\circ method." } @@ -147,7 +147,7 @@ "plot_range": "0.0, 0.3", "plot_scale": "linear", "retrieving_info": { - "eta": {"value": "__etaused[1064_FR]", "unit": ""} + "eta": {"value": "__etaused[1064_FR][eta]", "unit": "", "method": "__etaused[1064_FR][method]"} }, "comment": "The depolarization ratio was calibrated with \\Delta 90\\circ method." } diff --git a/ppcpy/interface/picassoProc.py b/ppcpy/interface/picassoProc.py index e8abdc3..febaa7a 100644 --- a/ppcpy/interface/picassoProc.py +++ b/ppcpy/interface/picassoProc.py @@ -3,6 +3,7 @@ import re import numpy as np import logging +from pathlib import Path import ppcpy.misc.pollyChannelTags as pollyChannelTags import ppcpy.preprocess.pollyPreprocess as pollyPreprocess import ppcpy.qc.pollySaturationDetect as pollySaturationDetect @@ -345,35 +346,128 @@ def SaturationDetect(self): data_cube = self, sigSaturateThresh = self.polly_config_dict['saturate_thresh']) - def polarizationCaliD90(self, db_path:str=None): - """Calibration with the Delta-90 method. + + def polarizationCaliD90(self, db_path:str=None, collect_debug:bool=False, **defaults): + """Calculate/Estimate and select optimal Depol calibration constants with + the Delta-90° method. + + The function calculates/estimates Depol calibration constants (DCs) through + the Delta-90° method. If configured, DCs from a defined database are also loaded. + The optimal DC per channel is chosen as the DC with the lowest relative error + from both the calculated/estimated and loaded DCs. If no DC is found for a + channel its default value is taken. Parameters ---------- - db_path : str - Path to database to read DC values from in case none - are successfully retrieved. Default is None. - + db_path : str, optional + Path to calibration constant database. If ``db_path`` is ``None`` + and config variable ``flagUsePreviousDepolCali`` is enabled, + a standard path is constructed from the config variable + ``calibrationDB``. Default is ``None``. + collect_debug : bool, optional + If True, collects debug information. Default is False. + polCaliEta355 : float, optional + Default depol calibration constant at 355 nm. + polCaliEtaStd355 : float, optional + Default depol calibration constant error at 355 nm. + polCaliEta532 : float, optional + Default depol calibration constant at 532 nm. + polCaliEtaStd532 : float, optional + Default depol calibration constant error at 532 nm. + polCaliEta1064 : float, optional + Default depol calibration constant at 1064 nm. + polCaliEtaStd1064 : float, optional + Default depol calibration constant error at 1064 nm. + + Attributes + ---------- + self.pol_cali['D90' & 'D90_db'] : dict + All retrieved and read delta-90° depol calibration constants + and retrieval information. + self.pol_cali['default'] : dict + Default depol calibration constants and retrieval information. + self.etaused : dict + The optimal depol calibration constants per channel. + + Notes + ----- + The function uses the following configuration flags: + + - ``flagLCCalibration``: Enables or disables calculation of Depol + calibration constants through Delta-90° method. + + - ``flagUsePreviousLC``: Enables or disables loading previously determined + Depol calibration constants from the calibration + database. + + All calibration constants are stored in ``self.pol_cali``. + The selected optimal constants are stored in ``self.etaused``. + + Default values are by standard taken from their config variable but can be + overwritten if passed as an input to this function. + The stuff that starts here in the matlab version https://github.com/PollyNET/Pollynet_Processing_Chain/blob/5efd7d35596c67ef8672f5948e47d1f9d46ab867/lib/interface/picassoProcV3.m#L442 + + **History** + + - xxxx-xx-xx: First edition by ... + - 2026-10-08: Consider retrieved and database constants equally + and added default calibration constants. + """ + # Load GHK-parameters polarization.loadGHK(self) - self.pol_cali['D90'] = polarization.calibrateGHK(self) - isUsable = [element['status'] for key, val in self.pol_cali['D90'].items() for element in val] - - if np.sum(isUsable) > 0: - logging.info("Using retieved polarization calibration constants.") - self.etaused = select.single_best(self.pol_cali['D90'], 'eta', 'eta_std') - elif db_path is not None: - logging.warning("Can not retieve viable polarization calibration constants, uses constants form the database.") + + if self.polly_config_dict['flagDepolCali']: + logging.info("Performing Delta-90° Depol calibration ...") + + # Estimate Depol calibration constants + self.pol_cali['D90'] = polarization.calibrateGHK( + self, + collect_debug=collect_debug + ) + all_D90_pol_calis = self.pol_cali['D90'].copy() + + # Check retrieval + isUsable = sum([True for val in self.pol_cali['D90'].values() for element in val if element]) + if not isUsable: + logging.warning("No viable depol calibration constants was retrieved from the measurement.") + + else: + logging.warning("Delta-90° Depol calibration calibration is turned off.") + all_D90_pol_calis = {} + + if self.polly_config_dict['flagUsePreviousDepolCali']: + if db_path is None: + base_dir = Path(self.picasso_config_dict['results_folder']) + db_path = base_dir.joinpath(self.device, self.polly_config_dict['calibrationDB']) + logging.info(f"Loading Depol calibration constants from database: {db_path}") + + # Load DCs from database table_name = 'depol_calibration_constant' ts_interval = self.retrievals_highres['time'][0], self.retrievals_highres['time'][-1] - self.pol_cali['D90_db'] = sql_db.get_from_sql_db(db_path, table_name, ts_interval)['D90_db'] - self.etaused = select.single_best(self.pol_cali['D90_db'], 'eta', 'eta_std') - else: - logging.critical("Can not retieve viable polarization calibration constants, and no database detected.") - raise ValueError("Can not retieve viable polarization calibration constants, and no database detected.") + db_table = sql_db.get_from_sql_db(db_path, table_name, ts_interval) + if 'D90_db' in db_table: + self.pol_cali['D90_db'] = db_table['D90_db'] + + # Combine retrieved and database D90 DCs + for ch in self.pol_cali['D90_db'].keys(): + all_D90_pol_calis[ch] = all_D90_pol_calis.get(ch, []) \ + + self.pol_cali['D90_db'][ch] + + # Load default DCs + self.pol_cali['default'] = polarization.loadDefaults(self, **defaults) + + # Select optimal DCs + self.etaused = select.single_best(self.pol_cali['default'], 'eta', 'eta_std', 'method', relative=True) |\ + select.single_best(all_D90_pol_calis, 'eta', 'eta_std', 'method', relative=True) + + # Check if default DCs were used + defaultDCsUsed = [ch for ch in self.etaused if self.etaused[ch]['method'] == 'default'] + if defaultDCsUsed: + logging.warning(f"Deafault Depol calibration constant used for channels: {defaultDCsUsed}") def cloudScreen(self, collect_debug:bool=False): @@ -656,40 +750,136 @@ def Angstroem(self): self.retrievals_profile[ret_prof_name] = angstroem.ae_cldFreeGrps( self, ret_prof_name) - def LidarCalibration(self, db_path:str=None, collect_debug:bool=False): - """calculate the lidar constant - .. TODO:: Find out how we prioritise raman, klett, and database retrieved LC... + def LidarCalibration(self, db_path:str=None, collect_debug:bool=False, **defaults): + """Calculate/Estimate and select optimal lidar calibration constants. + + The function calculates/estimates lidar calibration constants (LCs) from + both klett and Raman retrieved optical profiles. If configured, LCs from + a defined database are also loaded. The optimal LC per channel is chosen as + the LC with the lowest relative error per channel form both the + calculated/estimated and loaded LCs. If no LC is found for a channel its + default value is taken. + + Parameters + ---------- + db_path : str or Path, optional + Path to calibration constant database. If ``db_path`` is ``None`` + and config variable ``flagUsePreviousLC`` is enabled, a standard + path is constructed from the config variable ``calibrationDB``. + Default is ``None``. + collect_debug : bool, optional + If ``True``, collects debug information. Default is ``False``. + LC : list, optional + Default Lidar constant value per channel. + LCStd : list, optional + Default Lidar constant error per channel. + + Attributes + ---------- + self.LC['klett' & 'klett_db'] : dict + All retrieved and read Klett lidar calibration constants + and retrieval information. + self.LC['raman' & 'raman_db'] : dict + All retrieved and read Raman lidar calibration constants + and retrieval information. + self.LC['default'] : dict + Default lidar calibration constants and retrieval information. + self.LCused : dict + The optimal lidar calibration constant per channel. + + Notes + ----- + The function uses the following configuration flags: + + - ``flagLCCalibration``: Enables or disables calculation of Lidar + calibration constants from klett and Raman retrieval. + + - ``flagUsePreviousLC``: Enables or disables loading previously determined + lidar calibration constants from the calibration + database. + + All calibration constants are stored in ``self.LC``. + The selected optimal constants are stored in ``self.LCused``. + + Default values are by standard taken from their config variable but can be + overwritten if passed as an input to this function. + + **History** + + - xxxx-xx-xx: First edition by ... + - 2026-10-08: Consider retrieved and database constants equally + and added default calibration constants. + """ - self.LC['klett'] = lidarconstant.lc_for_cldFreeGrps( - self, - retrieval='klett', - collect_debug=collect_debug - ) - self.LC['raman'] = lidarconstant.lc_for_cldFreeGrps( - self, - retrieval='raman', - collect_debug=collect_debug - ) - logging.info("Choosing best LC per channel...") - if db_path is None: - logging.info("No database path found. Using retrieved LC values.") - self.LC['klett_db'] = {} - self.LC['raman_db'] = {} + if self.polly_config_dict['flagLCCalibration']: + logging.info("Performing Lidar calibration ...") + + # Estimate Lidar constant from klett retrieval + self.LC['klett'] = lidarconstant.lc_for_cldFreeGrps( + self, + retrieval='klett', + collect_debug=collect_debug + ) + all_klett_LCs = self.LC['klett'].copy() + + # Estimate Lidar constant from raman retrieval + self.LC['raman'] = lidarconstant.lc_for_cldFreeGrps( + self, + retrieval='raman', + collect_debug=collect_debug + ) + all_raman_LCs = self.LC['raman'].copy() + + # Check retrieval + isUsable = sum([True for val in list(self.LC['klett'].values()) \ + + list(self.LC['raman'].values()) for element in val if element]) + if not isUsable: + logging.warning("No viable Lidar calibration constants was retrieved from the measurement.") + else: - logging.info("Database LC values will be used when no retrieved ones are available.") - # db_path = self.polly_config_dict['calibrationDB'] + logging.warning("Lidar calibration is turned off.") + all_klett_LCs, all_raman_LCs = {}, {} + + if self.polly_config_dict['flagUsePreviousLC']: + if db_path is None: + base_dir = Path(self.picasso_config_dict['results_folder']) + db_path = base_dir.joinpath(self.device, self.polly_config_dict['calibrationDB']) + logging.info(f"Loading Lidar calibration constants from database: {db_path}") + + # Load LCs from database table_name = 'lidar_calibration_constant' ts_interval = self.retrievals_highres['time'][0], self.retrievals_highres['time'][-1] - self.LC['klett_db'] = sql_db.get_from_sql_db(db_path, table_name, ts_interval)['klett_db'] - self.LC['raman_db'] = sql_db.get_from_sql_db(db_path, table_name, ts_interval)['raman_db'] - - # Prioritise Raman retrieved LCs but use Klett retrieved ones when no Raman retrieval exists. - self.LCused = select.single_best(self.LC['klett_db'], 'LC', 'LCStd', relative=True) |\ - select.single_best(self.LC['raman_db'], 'LC', 'LCStd', relative=True) |\ - select.single_best(self.LC['klett'], 'LC', 'LCStd', relative=True) |\ - select.single_best(self.LC['raman'], 'LC', 'LCStd', relative=True) + db_table = sql_db.get_from_sql_db(db_path, table_name, ts_interval) + if 'klett_db' in db_table: + self.LC['klett_db'] = db_table['klett_db'] + + # Combine retrieved and database Klett LCs + for ch in self.LC['klett_db']: + all_klett_LCs[ch] = all_klett_LCs.get(ch, []) \ + + self.LC['klett_db'][ch] + + if 'raman_db' in db_table: + self.LC['raman_db'] = db_table['raman_db'] + + # Combine retrieved and database Raman LCs + for ch in self.LC['raman_db']: + all_raman_LCs[ch] = all_raman_LCs.get(ch, []) \ + + self.LC['raman_db'][ch] + + # Load default LCs + self.LC['default'] = lidarconstant.loadDefaults(self, **defaults) + + # Select optimal LCs + self.LCused = select.single_best(self.LC['default'], 'LC', 'LCStd', 'method', relative=True) |\ + select.single_best(all_klett_LCs, 'LC', 'LCStd', 'method', relative=True) |\ + select.single_best(all_raman_LCs, 'LC', 'LCStd', 'method', relative=True) + + # Check if default LCs were used + defaultLCsUsed = [ch for ch in self.LCused if self.LCused[ch]['method'] == 'default'] + if defaultLCsUsed: + logging.warning(f"Deafault Lidar calibration constant used for channels: {defaultLCsUsed}") def attBsc_volDepol(self): diff --git a/ppcpy/io/sql_interaction.py b/ppcpy/io/sql_interaction.py index 91cf7ee..404a912 100644 --- a/ppcpy/io/sql_interaction.py +++ b/ppcpy/io/sql_interaction.py @@ -3,6 +3,7 @@ from datetime import datetime, timezone, timedelta from collections import defaultdict import pandas as pd +import numpy as np import ppcpy.misc.helper as helper from ppcpy.misc.helper import default_to_regular @@ -31,7 +32,16 @@ def get_from_sql_db(db_path:str, table_name:str, ts_interval:list[str]) -> dict: ------- dict in calibration storage format + + **History** + + - 2026-03-28: First edition by Radenz. + - 2026-05-08: Added check for table existence. Added retrieval methods. + """ + + ret = {} + delta = timedelta(hours=24) start = ( datetime.fromtimestamp(ts_interval[0], timezone.utc) - delta @@ -40,13 +50,20 @@ def get_from_sql_db(db_path:str, table_name:str, ts_interval:list[str]) -> dict: datetime.fromtimestamp(ts_interval[0], timezone.utc) + delta ).strftime("%Y-%m-%d %H:%M") with sqlite3.connect(db_path) as conn: + exists = conn.execute( + "SELECT 1 FROM sqlite_master WHERE type='table' AND name=?;", + (table_name,) + ).fetchone() + + if exists is None: + logging.warning(f"Database {db_path} not found.") + return ret + df = pd.read_sql_query( f'SELECT * FROM {table_name} WHERE cali_start_time BETWEEN ? AND ?;', conn, params=(start, end)) conn.close() - ret = {} - if table_name == 'lidar_calibration_constant': d = defaultdict(list) for index, row in df[df.cali_method == 'Raman_Method'].iterrows(): @@ -54,9 +71,11 @@ def get_from_sql_db(db_path:str, table_name:str, ts_interval:list[str]) -> dict: d[k].append({ 'LC': row['liconst'], 'LCStd': row['uncertainty_liconst'], 'time_start': int(string_to_ts(row['cali_start_time'])), - 'time_end': int(string_to_ts(row['cali_stop_time'])), + 'time_end': int(string_to_ts(row['cali_stop_time'])), + 'method': 'raman_db', }) ret['raman_db'] = default_to_regular(d) + logging.info(f"Loaded {len(ret['raman_db'])} lines from table 'lidar_calibration_constant' with 'Raman_Method'.") d = defaultdict(list) for index, row in df[df.cali_method == 'Klett_Method'].iterrows(): @@ -64,10 +83,11 @@ def get_from_sql_db(db_path:str, table_name:str, ts_interval:list[str]) -> dict: d[k].append({ 'LC': row['liconst'], 'LCStd': row['uncertainty_liconst'], 'time_start': int(string_to_ts(row['cali_start_time'])), - 'time_end': int(string_to_ts(row['cali_stop_time'])), + 'time_end': int(string_to_ts(row['cali_stop_time'])), + 'method': 'klett_db', }) - ret['klett_db'] = default_to_regular(d) + logging.info(f"Loaded {len(ret['klett_db'])} lines from table 'lidar_calibration_constant' with 'Klett_Method'.") if table_name == 'depol_calibration_constant': d = defaultdict(list) @@ -76,10 +96,13 @@ def get_from_sql_db(db_path:str, table_name:str, ts_interval:list[str]) -> dict: d[k].append({ 'eta': row['depol_const'], 'eta_std': row['uncertainty_depol_const'], 'time_start': int(string_to_ts(row['cali_start_time'])), - 'time_end': int(string_to_ts(row['cali_stop_time'])), + 'time_end': int(string_to_ts(row['cali_stop_time'])), + 'method': 'D90_db', + # 'status': 1, }) ret['D90_db'] = default_to_regular(d) - + logging.info(f"Loaded {len(ret['D90_db'])} lines from table 'depol_calibration_constant'.") + return ret def prepare_for_sql_db_writing(data_cube, parameter:str, method:str) -> list[tuple]: @@ -105,15 +128,14 @@ def prepare_for_sql_db_writing(data_cube, parameter:str, method:str) -> list[tup elif method == 'klett': method_db = 'Klett_Method' - print(data_cube.LC.keys()) if parameter == 'LC': - for e in data_cube.LC[method].keys(): + for e in data_cube.LC.get(method, {}).keys(): wv, pol, tel = helper.get_wv_pol_telescope_from_dictkeyname(e) tel_db = mapping_inverse[tel] for line in data_cube.LC[method][e]: LC = line['LC'] LC_std = line['LCStd'] - LC_is_used = True if LC == data_cube.LCused[e] else False + LC_is_used = True if LC == data_cube.LCused[e]['LC'] else False start_unix = line['time_start'] stop_unix = line['time_end'] start = datetime.fromtimestamp(start_unix, timezone.utc).strftime("%Y-%m-%d %H:%M:%S") @@ -123,13 +145,13 @@ def prepare_for_sql_db_writing(data_cube, parameter:str, method:str) -> list[tup wv, str(data_cube.rawfile), data_cube.device, method_db, tel_db)) elif parameter == 'DC': - for e in data_cube.pol_cali['D90'].keys(): + for e in data_cube.pol_cali.get('D90', {}).keys(): wv, tel = e.split('_') tel_db = mapping_inverse[tel] for line in data_cube.pol_cali['D90'][e]: eta = line['eta'] eta_std = line['eta_std'] - eta_is_used = True if eta == data_cube.etaused[e] else False + eta_is_used = True if eta == data_cube.etaused[e]['eta'] else False start_unix = line['time_start'] stop_unix = line['time_end'] start = datetime.fromtimestamp(start_unix, timezone.utc).strftime("%Y-%m-%d %H:%M:%S") @@ -137,6 +159,7 @@ def prepare_for_sql_db_writing(data_cube, parameter:str, method:str) -> list[tup rows_to_insert.append(( str(start), str(stop), float(eta), float(eta_std), eta_is_used, wv, tel_db, str(data_cube.rawfile), data_cube.device)) + return rows_to_insert def setup_empty(db_path:str, table_name:str, column_names:list[str], data_types:list[str], unique:str=''): diff --git a/ppcpy/io/write2nc.py b/ppcpy/io/write2nc.py index 824d641..3d98cc9 100644 --- a/ppcpy/io/write2nc.py +++ b/ppcpy/io/write2nc.py @@ -88,10 +88,10 @@ def write_channelwise_2_nc_file(data_cube, root_dir=root_dir, prod_ls=[]): ## update variable attribute if "eta" in json_nc_mapping_dict['variables'][v]['attributes'].keys(): wv, t, tel = re.findall(r"(\d{3,4})_(\w+)_(\w+)", v)[0] - json_nc_mapping_dict['variables'][v]['attributes']['eta'] = data_cube.etaused[f'{wv}_{tel}'] + json_nc_mapping_dict['variables'][v]['attributes']['eta'] = data_cube.etaused[f'{wv}_{tel}']['eta'] if "Lidar_calibration_constant_used" in json_nc_mapping_dict['variables'][v]['attributes'].keys(): LC_used_key = v.split("attBsc_")[-1] - json_nc_mapping_dict['variables'][v]['attributes']['Lidar_calibration_constant_used'] = data_cube.LCused[LC_used_key] + json_nc_mapping_dict['variables'][v]['attributes']['Lidar_calibration_constant_used'] = data_cube.LCused[LC_used_key]['LC'] ### remove empty key-value-pairs if json_nc_mapping_dict['variables'][v]['data'] is None: json2nc_mapping.remove_variable_from_json_dict_mapper(data_dict=json_nc_mapping_dict, key_to_remove=v) diff --git a/ppcpy/qc/transCor.py b/ppcpy/qc/transCor.py index 46ccf1d..3622c7c 100644 --- a/ppcpy/qc/transCor.py +++ b/ppcpy/qc/transCor.py @@ -30,14 +30,15 @@ def transCorGHK_cube(data_cube, signal='BGCor'): print('G', config_dict['G'][flagt], config_dict['G'][flagc]) print('H', config_dict['H'][flagt], config_dict['H'][flagc]) - print('polCaliEta', data_cube.etaused[f'{wv}_{tel}']) + print('polCaliEta', data_cube.etaused[f'{wv}_{tel}']['eta']) # similar to voldepol_2d vdr, vdrStd = depolarization.calc_profile_vdr( sigBGCor_total, sigBGCor_cross, config_dict['G'][flagt], config_dict['G'][flagc], config_dict['H'][flagt], config_dict['H'][flagc], - data_cube.etaused[f'{wv}_{tel}'], config_dict[f'voldepol_error_{wv}'], + data_cube.etaused[f'{wv}_{tel}']['eta'], + config_dict[f'voldepol_error_{wv}'], ) sigTCor_total, bgTCor_total = transCor_E16_channel( diff --git a/ppcpy/retrievals/depolarization.py b/ppcpy/retrievals/depolarization.py index fbff880..806388f 100644 --- a/ppcpy/retrievals/depolarization.py +++ b/ppcpy/retrievals/depolarization.py @@ -63,12 +63,13 @@ def voldepol_cldFreeGrps(data_cube, ret_prof_name): # data_cube.retrievals_highres[f'BG{signal}'][slice(*cldFree),data_cube.gf(wv, 'total', tel)]), axis=0) sigc = np.squeeze(data_cube.retrievals_profile[f'sig{signal}'][i,:,flagc]) - print(channel, data_cube.etaused[f'{wv}_{tel}']) + print(channel, data_cube.etaused[f'{wv}_{tel}']['eta']) vdr, vdrStd = calc_profile_vdr( sigt, sigc, config_dict['G'][flagt], config_dict['G'][flagc], config_dict['H'][flagt], config_dict['H'][flagc], - data_cube.etaused[f'{wv}_{tel}'], config_dict[f'voldepol_error_{wv}'], + data_cube.etaused[f'{wv}_{tel}']['eta'], + config_dict[f'voldepol_error_{wv}'], window=config_dict[f'smoothWin_{retrieval}_{wv}'] ) opt_profiles[i][channel]['vdr'] = vdr @@ -90,7 +91,8 @@ def voldepol_cldFreeGrps(data_cube, ret_prof_name): vdr, vdrStd = calc_profile_vdr( sigt, sigc, config_dict['G'][flagt], config_dict['G'][flagc], config_dict['H'][flagt], config_dict['H'][flagc], - data_cube.etaused[f'{wv}_{tel}'], config_dict[f'voldepol_error_{wv}'], + data_cube.etaused[f'{wv}_{tel}']['eta'], + config_dict[f'voldepol_error_{wv}'], window=1 ) mdr, mdrStd, flgaDeftMdr = get_MDR( diff --git a/ppcpy/retrievals/highres.py b/ppcpy/retrievals/highres.py index 42d83be..a2800a9 100644 --- a/ppcpy/retrievals/highres.py +++ b/ppcpy/retrievals/highres.py @@ -44,7 +44,7 @@ def attbsc_2d(data_cube, nr:bool=True, collect_debug:bool=False): else: logging.info(f'{channel} skipped at attbsc_2d') continue - attBsc = sig * ranges2d / data_cube.LCused[channel] + attBsc = sig * ranges2d / data_cube.LCused[channel]['LC'] attBsc[data_cube.retrievals_highres['depCalMask'], :] = np.nan data_cube.retrievals_highres[f"attBsc_{channel}"] = attBsc @@ -68,7 +68,7 @@ def attbsc_2d(data_cube, nr:bool=True, collect_debug:bool=False): logging.info(f'{channel} skipped at attbsc_2d OL') continue - attBsc = sig * ranges2d / data_cube.LCused[channel] + attBsc = sig * ranges2d / data_cube.LCused[channel]['LC'] attBsc[data_cube.retrievals_highres['depCalMask'], :] = np.nan data_cube.retrievals_highres[f"attBsc_{wv}_{t}_OC"] = attBsc @@ -109,7 +109,8 @@ def voldepol_2d(data_cube): vdr, vdrStd = depolarization.calc_profile_vdr( sigt, sigc, config_dict['G'][flagt], config_dict['G'][flagc], config_dict['H'][flagt], config_dict['H'][flagc], - data_cube.etaused[f'{wv}_{tel}'], config_dict[f'voldepol_error_{wv}'], + data_cube.etaused[f'{wv}_{tel}']['eta'], + config_dict[f'voldepol_error_{wv}'], window=1) vdr[data_cube.retrievals_highres['depCalMask'], :] = np.nan data_cube.retrievals_highres[f"voldepol_{wv}_total_{tel}"] = vdr diff --git a/ppcpy/retrievals/quasi.py b/ppcpy/retrievals/quasi.py index 0c0de04..c3f4e21 100644 --- a/ppcpy/retrievals/quasi.py +++ b/ppcpy/retrievals/quasi.py @@ -64,7 +64,7 @@ def quasi_pdr(data_cube, wvs:list=[532], version:str='V1'): sigt=sigt, sigc=sigc, Gt=config_dict['G'][flagt], Gr=config_dict['G'][flagc], Ht=config_dict['H'][flagt], Hr=config_dict['H'][flagc], - eta=data_cube.etaused[f"{wv}_{tel}"], + eta=data_cube.etaused[f"{wv}_{tel}"]['eta'], voldepol_error=config_dict[f'voldepol_error_{wv}'], window=1 )