Skip to content

90 Refine calibration constant selection. #91

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
HavardStridBuholdt wants to merge 1 commit into
base: main
Choose a base branch
from
Open

90 Refine calibration constant selection. #91

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
179 changes: 138 additions & 41 deletions ppcpy/calibration/lidarconstant.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
----------
Expand All @@ -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')
Expand All @@ -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_'
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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(
Expand All @@ -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)

Loading