implement the cal constants reading from database

docs/source/api_doc/calibration.rst

@@ -24,3 +24,10 @@ ppcpy.calibration.rayleighfit
 
 .. automodule:: ppcpy.calibration.rayleighfit
    :members:
+
+*************************
+ppcpy.calibration.select
+*************************
+
+.. automodule:: ppcpy.calibration.select
+   :members:

docs/source/userguide/cal_constants.rst

@@ -0,0 +1,117 @@
+
+*********************
+Calibration Constants
+*********************
+
+Polarization Calibration
+=========================
+
+The delta-90 polarization calibration for the detected time periods is done with :py:meth:`~ppcpy.interface.picassoProc.PicassoProc.polarizationCaliD90`.
+Resulting etas are stored in ``data_cube.pol_cali`` for all the time intervals and ``data_cube.etaused`` for the one with the lowest standard deviation.
+
+
+.. code-block:: python
+
+    >>> data_cube.polarizationCaliD90()
+    >>> data_cube.pol_cali
+    {'D90': 
+      {'355_FR': [
+        {'eta': 10.602253607032777,
+          'eta_std': 0.36947564349830514,
+          'time_start': 1755916590,
+          'time_end': 1755916800,
+          'status': 1},
+        {'eta': 10.45194101666947,
+          'eta_std': 0.27254858255926373,
+          'time_start': 1755968190,
+          'time_end': 1755968400,
+          'status': 1}, 
+    ... ]
+    >>> data_cube.etaused
+    {'355_FR': np.float64(10.754878969562773),
+     '532_FR': np.float64(6.852638549685832)}
+
+A method for polarization calibration at the reference height is also implemented, but not extensively tested or used operationally: :py:meth:`~ppcpy.interface.picassoProc.PicassoProc.polarizationCaliMol`
+
+
+Lidar (absolute) Calibration
+==============================
+
+After retrieval of the optical profiles, the Lidar calibration constants can be calculated with :py:meth:`~ppcpy.interface.picassoProc.PicassoProc.LidarCalibration`.
+Similarly to the polarization calibration the values are stored in ``data_cube.LC`` with values with the lowest standard deviation are selected into ``data_cube.LCused``.
+
+
+.. code-block:: python
+
+    >>> data_cube.LidarCalibration()
+    >>> data_cube.LC
+    {'klett': {
+      '532_total_FR': [
+        {'LC': np.float64(76651241245552.4),
+         'LCStd': np.float64(38384910061.7245),
+         'time_start': 1755907200, 'time_end': 1755910770},
+        {'LC': np.float64(79019754266917.53),
+         'LCStd': np.float64(20060545795.042507),
+         'time_start': 1755910800, 'time_end': 1755916200},
+        {'LC': np.float64(77383458139465.62),
+         'LCStd': np.float64(30826575129.720245),
+         'time_start': 1755916830, 'time_end': 1755920400},
+        ... ],
+      '355_total_FR': [
+        {'LC': np.float64(61221519590893.43),
+         'LCStd': np.float64(15266830837.91351),
+         'time_start': 1755907200, 'time_end': 1755910770},
+        ... ]},
+    'raman': {
+      '355_total_FR': [
+        {'LC': np.float64(57637355636083.1),
+         'LCStd': np.float64(42772775199.508766),
+         'time_start': 1755907200, 'time_end': 1755910770},
+        {'LC': np.float64(56121416303705.17),
+         'LCStd': np.float64(65519957166.51127),
+         'time_start': 1755910800, 'time_end': 1755916200},
+        ... ]},
+    }
+    >>> data_cube.LCused
+    {'532_total_FR': np.float64(47109585570330.64),
+     '355_total_FR': np.float64(51991322600216.1),
+     '1064_total_FR': np.float64(418368070699337.2),
+     '532_total_NR': np.float64(6649087258385.09),
+     '355_total_NR': np.float64(9550053707341.717),
+     '387_total_FR': np.float64(51991322600216.1),
+     '607_total_FR': np.float64(47109585570330.64),
+     '607_total_NR': np.float64(6649087258385.09),
+     '387_total_NR': np.float64(9550053707341.717)}
+
+
+Writing/Reading Database
+==========================
+
+A sqlite database for storage of the retrieved calibration data is defined in ``data_cube.polly_config_dict['calibrationDB']`` (originally in the respective config file): :py:meth:`~ppcpy.interface.picassoProc.PicassoProc.write_2_sql_db`
+
+.. code-block:: python
+
+    data_cube.write_2_sql_db(parameter='LC', method='raman')
+    data_cube.write_2_sql_db(parameter='LC', method='klett')
+    data_cube.write_2_sql_db(parameter='DC')
+
+Storing multiple times will update the calibration constants, if exactly the same time interval is used for the optical profile.
+The calibration constants can also be read back into picassopy: :py:meth:`~ppcpy.interface.picassoProc.PicassoProc.read_calibration_db`
+
+.. code-block:: python
+
+    data_cube.read_calibration_db()
+
+The values are then stored into ``data_cube.LC['raman_db']``, ``data_cube.LC['klett_db']`` and ``data_cube.pol_cali['D90_db']``, respectively.
+
+A basic visualization of the calibration constants is also available by :py:func:`~ppcpy.calibration.select.plot_cals`
+
+.. code-block:: python
+
+    from ppcpy.calibration.select import plot_cals
+    plot_cals(data_cube.pol_cali, 'eta', used=data_cube.etaused)
+    plot_cals(data_cube.LC, 'LC', used=data_cube.LCused)
+
+
+.. image:: img/plot_cals_example.png
+    :width: 55%

docs/source/userguide/index.rst

@@ -12,5 +12,7 @@ User Guide
    configfiles.rst
    preprocess.rst
 
+   cal_constants.rst
+
    meteodata.rst
    todo.rst

docs/source/userguide/preprocess.rst

@@ -15,10 +15,10 @@ Each of these quantities has to be set for each channel (mostly depending on the
 
 .. code-block:: python
 
-    data_cube.polly_config_dict['first_range_gate_indx']
-    >>> [252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252]
-    data_cube.polly_config_dict['first_range_gate_height']
-    >>> [3.75, 3.75, 3.75, 3.75, 3.75, 3.75, 3.75, 3.75, 6.75, 6.75, 6.75, 6.75, 6.75]
+    >>> data_cube.polly_config_dict['first_range_gate_indx']
+    [252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252]
+    >>> data_cube.polly_config_dict['first_range_gate_height']
+    [3.75, 3.75, 3.75, 3.75, 3.75, 3.75, 3.75, 3.75, 6.75, 6.75, 6.75, 6.75, 6.75]
 
 To check the correct setting of the bin height, the ``preproSignal`` and the ``sigBGCor`` can be compared close to the trigger
 
@@ -65,7 +65,7 @@ Where the ``preproSignal`` (left) still contains the pretrigger and the laser pu
 that contains a 0 count rate.
 
 .. image:: img/first_range_gate_indx.png
-  :width: 100%
+    :width: 100%
 
 
 The ``first_range_gate_height`` affects the range corrected signal as well as the mapping from time to range. In the original labview
@@ -110,6 +110,6 @@ Later, when the range corrected signal is calculated, the the count rates are co
 
 
 .. image:: img/first_range_gate_height_correction.png
-  :width: 100%
+    :width: 100%
 
 

ppcpy/calibration/lidarconstant.py

@@ -4,8 +4,10 @@
 
 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}
 
@@ -40,9 +42,10 @@ def lc_for_cldFreeGrps(data_cube, retrieval:str) -> list:
     print('LCMeanWindow', config_dict['LCMeanWindow'], 
           'LCMeanMinIndx', config_dict['LCMeanMinIndx'],
           'LCMeanMaxIndx', config_dict['LCMeanMaxIndx'])
-    LCs = [{} for i in range(len(data_cube.clFreeGrps))]
+    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]
 
@@ -94,7 +97,10 @@ def lc_for_cldFreeGrps(data_cube, retrieval:str) -> list:
                 print(f'Can not find a stable LC value, skipping {channel} {cldFree}')
                 continue
 
-            LCs[i][channel] = {'LC': LC_stable, 'LCStd': LC_stable * LCStd}
+            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)] 
@@ -124,43 +130,10 @@ def lc_for_cldFreeGrps(data_cube, retrieval:str) -> list:
                     print(f'Can not find a stable LC value, skipping {channel} {cldFree}')
                     continue
 
-                LCs[i][f"{wv_r}_{t}_{tel}"] = {'LC': LC_r_stable, 'LCStd': LC_r_stable * LCStd_r}
-
-    return LCs
-
-
-def get_best_LC(LCs:list) -> dict:
-    """Get lidar constant with the lowest standard deviation.
-
-    Parameters
-    ----------
-    LCs : list
-        Lidar constant for each channel per cloud free period.
-
-    Returns
-    -------
-    LCused : dict
-        Lidar constants with lowest standard deviation per channel.
-
-    Notes
-    -----
-    - Since ``LC = LC_stable`` and ``LCStd = LC_stable * LC_Std`` so will any negative LC also have
-      a negative LCStd, and thus be chosen as the best LC.
-
-    **History**
-
-    - 2026-02-16: Added additional checks to hinder negative LCs to be chosen.
-
-    """
-    # list comprehension for nested list
-    all_channels = set([k for e in LCs for k in e.keys()])
-
-    LCused = {}
-    for channel in all_channels:
-        lcs = np.array([e[channel]['LC'] for e in LCs if channel in e and e[channel]['LC'] >= 0])
-        lcsstd = np.array([e[channel]['LCStd'] for e in LCs if channel in e and e[channel]['LCStd'] >= 0])
-
-        LCused[channel] = lcs[np.argmin(lcsstd)]
-    return LCused
+                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)
 

ppcpy/calibration/polarization.py

@@ -6,6 +6,7 @@
 
 from ppcpy.misc.helper import uniform_filter
 from ppcpy.retrievals.collection import calc_snr
+from ppcpy.misc.helper import default_to_regular
 
 
 # Helper functions
@@ -176,7 +177,8 @@ def calibrateGHK(data_cube):
                 data_cube.polly_config_dict[f'rel_std_dminus_{wv}'],
                 data_cube.polly_config_dict[f'depol_cal_segmentLen_{wv}'],
                 data_cube.polly_config_dict[f'depol_cal_smoothWin_{wv}'], collect_debug=False)
-            logging.info(f'pol_cali_{wv}',  pol_cali[f'{wv}_{tel}'])
+            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')
 
@@ -371,8 +373,10 @@ def depol_cali_ghk(signal_t, bg_t, signal_x, bg_x, time, pol_cali_pang_start_tim
     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': pol_cali_eta, 'eta_std': pol_cali_eta_std, 'time_start': pol_cali_start_time, 'time_end': pol_cali_stop_time, 'status': 1}
-    results['eta_best'] = pol_cali_eta[np.argmin(pol_cali_eta_std)]
+    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)]
+
     if collect_debug:
         results['global_attri'] = dict(global_attri)
     return results
@@ -422,17 +426,13 @@ def analyze_segments(dplus, dminus, segment_len, rel_std_dplus, rel_std_dminus):
     data.polCaliEta532=data.polCali532Attri.polCaliEta(index_min);
 """
 
-def default_to_regular(d):
-    """ """
-    if isinstance(d, defaultdict):
-        d = {k: default_to_regular(v) for k, v in d.items()}
-    return d
-
-
 def calibrateMol(data_cube):
     """calibrate the polarization with the molecular signal
 
     Converted from the matlab code to the best knowledge, but not cross-validated yet
+
+    .. todo::
+        had to calculate TR_t, TR_c again when calling depol_cali_mol()
 
     """
 
@@ -469,22 +469,18 @@ def calibrateMol(data_cube):
                     background_t=bg_total, 
                     signal_c=sigBGCor_cross[:, slice(*refHInd)],
                     background_c=bg_cross,
-                    TR_t=np.squeeze(data_cube.polly_config_dict['TR'][data_cube.gf(wv, t, tel)]),
+                    TR_t=onemx_onepx(np.squeeze(data_cube.polly_config_dict['H'][data_cube.gf(wv, t, tel)])),
                     TR_t_std=0,
-                    TR_c=np.squeeze(data_cube.polly_config_dict['TR'][data_cube.gf(wv, 'cross', tel)]),
+                    TR_c=onemx_onepx(np.squeeze(data_cube.polly_config_dict['H'][data_cube.gf(wv, 'cross', tel)])),
                     TR_c_std=0,
                     minSNR=10,
                     mdr=config_dict[f'molDepol{wv}'],
                     mdrStd=config_dict[f'molDepolStd{wv}'],
                 )
+                ret['time_start'] = int(cldFreeTime[0])
+                ret['time_end'] = int(cldFreeTime[1])
                 if not ret['status'] == 0:
-                    pol_cali[f'{wv}_{tel}']['eta'].append(ret['eta'])
-                    pol_cali[f'{wv}_{tel}']['eta_std'].append(ret['eta_std'])
-                    pol_cali[f'{wv}_{tel}']['fac'].append(ret['fac'])
-                    pol_cali[f'{wv}_{tel}']['fac_std'].append(ret['fac_std'])
-                    pol_cali[f'{wv}_{tel}']['time_start'].append(int(cldFreeTime[0]))
-                    pol_cali[f'{wv}_{tel}']['time_end'].append(int(cldFreeTime[1]))
-                    pol_cali[f'{wv}_{tel}']['status'] = 1
+                    pol_cali[f'{wv}_{tel}'].append(ret)
 
     return default_to_regular(pol_cali)