Skip to content

Commit e85bef0

Browse files
HolgerPollyNetHolgerPollyNet
authored
Merge pull request #297 from PollyNET/QC_profiles
Cleaned up scattering calculations in low mode
2 parents db03123 + 13eec5b commit e85bef0

9 files changed

Lines changed: 1153 additions & 954 deletions

‎lib/config/polly_global_config.json‎

Lines changed: 299 additions & 299 deletions
Large diffs are not rendered by default.

‎lib/interface/picassoProcV3.m‎

Lines changed: 431 additions & 651 deletions
Large diffs are not rendered by default.

‎lib/io/pollySaveProfiles_QC.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ function pollySaveProfiles_QC(data)
2020
flagCh1064FR = PollyConfig.isFR & PollyConfig.is1064nm & PollyConfig.isTot;
2121

2222
missing_value = -999;
23-
23+
%%%%%%%%%%%%%%%%%%%%%here QC starts%%%%%%%%%%%%%%%%%%%%%
2424
for iGrp = 1:size(data.clFreGrps, 1)
2525
no_fill_low_profile =false;
2626
% QC flags will yet be implemetned only for standard products

‎lib/qc/pollyOLCor.m‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -150,7 +150,8 @@
150150
height, height(p.Results.normRange));
151151
bgOLCor = sigGlue(bgFR, p.Results.bgNR, p.Results.signalRatio, height, height(p.Results.normRange));
152152
end
153-
153+
case 4
154+
% signal gleuing according to maria T script
154155

155156
otherwise
156157
error('Unknown overlap correction mode %d', p.Results.overlapCorMode);

‎lib/retrievals/pollyRamanBsc.m‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -70,7 +70,7 @@
7070

7171
if wavelength == 355
7272
ext_aer_factor = (355 / 387) .^ angstroem;
73-
ext_mol_factor = (355 / 387) ^ 4;
73+
ext_mol_factor = (355 / 387) ^ 4; %ist nicht konsistent sollte explizit gerechnet werden
7474
elseif wavelength == 532
7575
ext_aer_factor = (532 / 607) .^ angstroem;
7676
ext_mol_factor = (532 / 607) ^ 4;
@@ -166,7 +166,7 @@
166166
aer_vr_OD = nansum(ext_aer_raman(1:refIndx)) * dH - ...
167167
nancumsum(ext_aer_raman) * dH;
168168

169-
hIndx = false(1, length(height));
169+
hIndx = false(1, length(height));
170170
hIndx(HRefIndx(1):HRefIndx(end)) = true;
171171

172172
% calculate the signal ratio at the reference height

‎lib/retrievals/pollyRamanBsc_smart.m‎

Lines changed: 120 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,120 @@
1+
function [ beta_aer, LR ] = pollyRamanBsc_smart(height, sigElastic, sigVRN2, ext_aer, angstroem, ext_mol, beta_mol,ext_mol_raman, beta_mol_inela, HRef, wavelength, betaRef, window_size, flagSmoothBefore, el_lambda, inel_lambda)
2+
% POLLYRAMANBSC calculate the aerosol backscatter coefficient with Raman method.
3+
%Modified by HB for cinsitency in 2024
4+
% USAGE:
5+
% [ beta_aer, LR ] = pollyRamanBsc(height, sigElastic, sigVRN2, ext_aer, ext_mol, beta_mol, HRef, wavelength, betaRef, window_size, flagSmoothBefore)
6+
%
7+
% INPUTS:
8+
% height: array
9+
% height. [m]
10+
% sigElastic: array
11+
% elastic photon count signal.
12+
% sigVRN2: array
13+
% N2 vibration rotational raman photon count signal.
14+
% ext_aer: array
15+
% aerosol extinction coefficient. [m^{-1}]
16+
% angstroem: array
17+
% aerosol angstroem exponent.
18+
% ext_mol: array
19+
% molecular extinction coefficient. [m^{-1}]
20+
% beta_mol: array
21+
% molecular backscatter coefficient. [m^{-1}Sr^{-1}]
22+
% HRef: 2 element array
23+
% reference region. [m]
24+
% wavelength: integer
25+
% wavelength of the corresponding elastic signal. [nm]
26+
% betaRef: float
27+
% aerosol backscatter coefficient at the reference region.
28+
% [m^{-1}Sr^{-1}]
29+
% window_size: integer
30+
% number of the bins of the sliding window for the signal smooth.
31+
% [default: 40]
32+
% flagSmoothBefore: logical
33+
% flag bit to control whether to smooth the signal before or after
34+
% calculating the signal ratio.
35+
% el_lambda: integer
36+
% elastic wavelength in nm
37+
% inel_lambda: interger
38+
% inelastic wavelength in nm
39+
% OUTPUTS:
40+
% beta_aer: array
41+
% aerosol backscatter coefficient. [m^{-1}Sr^{-1}]
42+
% LR: array
43+
% aerosol lidar ratio.
44+
%
45+
% REFERENCES:
46+
% netcdf-florian
47+
% Ansmann, A., et al. (1992). "Independent measurement of extinction and backscatter profiles in cirrus clouds by using a combined Raman elastic-backscatter lidar." Applied optics 31(33): 7113-7131.
48+
%
49+
% HISTORY:
50+
% - 2018-01-02: First edition by Zhenping.
51+
% - 2018-07-24: Add the ext_mol_factor and ext_aer_factor for wavelength of 1064nm
52+
% - 2018-09-04: Change the smoothing order. Previous way is smoothing the signal. This will create large drift at the signal ridges.
53+
% - 2018-09-05: Keep the original smoothing order for 355, which makes the retrieving results at the far range bins quite stable.
54+
% - 2024-11-12: Modified by HB for cinsitency in 2024
55+
% .. Authors: - HolgerPollyNet
56+
57+
% check the inputs
58+
if ~ (nargin >= 9)
59+
error('Not enough input arguments.');
60+
end
61+
62+
if ~ exist('window_size', 'var')
63+
window_size = 40;
64+
end
65+
66+
if (HRef(1) >= height(end)) || (HRef(end) <= height(1))
67+
error('HRef is out of range.');
68+
end
69+
70+
if ~ exist('flagSmoothBefore', 'var')
71+
flagSmoothBefore = true;
72+
end
73+
74+
75+
ext_aer_factor = (el_lambda / inel_lambda) .^ angstroem;
76+
77+
dH = height(2) - height(1); % height resolution. [m]
78+
79+
% find the index of the HRef region and the midpoint of HRef
80+
HRefIndx = [fix((HRef(1) - height(1)) / dH) + 1, ...
81+
fix((HRef(end) - height(1)) / dH) + 1];
82+
refIndx = fix((mean(HRef) - height(1)) / dH) + 1;
83+
84+
% calculate the extinction coefficient at inelastic wavelength.
85+
ext_aer_raman = ext_aer .* ext_aer_factor;
86+
87+
% calculate the optical depth from any point to refIndx
88+
mol_el_OD = nansum(ext_mol(1:refIndx)) * dH - nancumsum(ext_mol) * dH;
89+
mol_vr_OD = nansum(ext_mol_raman(1:refIndx)) * dH - ...
90+
nancumsum(ext_mol_raman) * dH;
91+
aer_el_OD = nansum(ext_aer(1:refIndx)) * dH - nancumsum(ext_aer) * dH;
92+
aer_vr_OD = nansum(ext_aer_raman(1:refIndx)) * dH - ...
93+
nancumsum(ext_aer_raman) * dH;
94+
95+
hIndx = false(1, length(height));
96+
hIndx(HRefIndx(1):HRefIndx(end)) = true;
97+
98+
% calculate the signal ratio at the reference height
99+
elMean = sigElastic(hIndx) ./ ...
100+
(beta_mol(hIndx) + betaRef);
101+
vrMean = sigVRN2(hIndx) ./ (beta_mol(hIndx));
102+
103+
% calculate the aerosol backscatter coefficient.
104+
if ~ flagSmoothBefore
105+
beta_aer = transpose(smoothWin(((sigElastic ./ sigVRN2) ...
106+
.* (nanmean(vrMean) ./ nanmean(elMean)) .* ...
107+
exp(mol_vr_OD - mol_el_OD + aer_vr_OD - aer_el_OD) - 1) .* ...
108+
beta_mol, window_size, 'moving'));
109+
else
110+
beta_aer = ...
111+
((transpose(smoothWin(sigElastic, window_size, 'moving')) ./ ...
112+
transpose(smoothWin(sigVRN2, window_size, 'moving'))) ...
113+
.* (nanmean(vrMean) ./ nanmean(elMean)) .* ...
114+
exp(mol_vr_OD - mol_el_OD + aer_vr_OD - aer_el_OD) - 1) .* beta_mol;
115+
end
116+
LR = ext_aer ./ beta_aer;
117+
118+
119+
120+
end

‎lib/retrievals/pollyRamanBsc_smart_MC.m‎

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
function [beta_aer, aerBscStd, LR] = pollyRamanBsc_smart_MC(height, sigElastic, sigVRN2, ext_aer, angstroem, ext_mol, beta_mol, ext_mol_raman, beta_mol_inela, HRef, betaRef, window_size, flagSmoothBefore, el_lambda, inel_lambda, bgElastic, bgVRN2,sigma_ext_aer, sigma_angstroem, MC_count, method)
2+
3+
% POLLYRAMANBSCSTD calculate uncertainty of aerosol backscatter coefficient with Monte-Carlo simulation.
4+
%
5+
% USAGE:
6+
% [ aerBscStd ] = pollyRamanBscStd(height, sigElastic, bgElastic, ...
7+
% sigVRN2, bgVRN2, ext_aer, sigma_ext_aer, angstroem, ...
8+
% sigma_angstroem, ext_mol, beta_mol, HRef, wavelength, ...
9+
% betaRef, smoothWindow, MC_count, method, flagSmoothBefore)
10+
%
11+
% INPUTS:
12+
% height: array
13+
% height. [m]
14+
% sigElastic: array
15+
% elastic photon count signal.
16+
% bgElastic: array
17+
% background of elastic signal.
18+
% sigVRN2: array
19+
% N2 vibration rotational raman photon count signal.
20+
% bgVRN2: array
21+
% background of N2 vibration rotational signal.
22+
% ext_aer: array
23+
% aerosol extinction coefficient. [m^{-1}]
24+
% sigma_ext_aer: array
25+
% uncertainty of aerosol extinction coefficient. [m^{-1}]
26+
% angstroem: array
27+
% aerosol angstroem exponent.
28+
% ext_mol: array
29+
% molecular extinction coefficient. [m^{-1}]
30+
% beta_mol: array
31+
% molecular backscatter coefficient. [m^{-1}Sr^{-1}]
32+
% HRef: 2 element array
33+
% reference region. [m]
34+
% wavelength: integer
35+
% wavelength of the corresponding elastic signal. [nm]
36+
% betaRef: float
37+
% aerosol backscatter coefficient at the reference region. [m^{-1}Sr^{-1}]
38+
% smoothWindow: integer or n*3 matrix
39+
% number of the bins of the sliding window for the signal smooth. [default: 40]
40+
% MC_count: scalar or matrix
41+
% samples for each error source. [samples_angstroem, samples_aerExt, samples_signal, samples_aerBscRef]
42+
% method: char
43+
% computational method. ['monte-carlo' or 'analytical']
44+
% flagSmoothBefore: logical
45+
% flag to control the smooth order.
46+
%
47+
% OUTPUTS:
48+
% aerBscStd: array
49+
% uncertainty of aerosol backscatter coefficient. [m^{-1}Sr^{-1}]
50+
%
51+
% HISTORY:
52+
% - 2021-07-16: first edition by Zhenping
53+
%
54+
% .. Authors: - zhenping@tropos.de
55+
56+
if ~ exist('method', 'var')
57+
method = 'monte-carlo';
58+
end
59+
60+
if ~ exist('flagSmoothBefore', 'var')
61+
flagSmoothBefore = true;
62+
end
63+
64+
if ~ exist('MC_count', 'var')
65+
MC_count = 3;
66+
end
67+
68+
if isscalar(MC_count)
69+
MC_count = ones(1, 4) * MC_count;
70+
end
71+
72+
if prod(MC_count) > 1e5
73+
warning('Too large sampling for monte-carlo simulation.');
74+
aerBscStd = NaN(size(sigElastic));
75+
return;
76+
end
77+
78+
[ beta_aer, LR ] = pollyRamanBsc_smart(height, sigElastic, sigVRN2, ext_aer, angstroem, ext_mol, beta_mol,ext_mol_raman, beta_mol_inela, HRef, el_lambda, betaRef, window_size, flagSmoothBefore, el_lambda, inel_lambda);
79+
80+
if strcmpi(method, 'monte-carlo')
81+
hRefIndx = (height >= HRef(1)) & (height < HRef(2));
82+
% rel_std_betaRef = std(sigElastic(hRefIndx)./sigVRN2(hRefIndx)) / mean(sigElastic(hRefIndx)./sigVRN2(hRefIndx)) * 0;
83+
rel_std_betaRef = 0.1;
84+
betaRefSample = transpose(sigGenWithNoise(betaRef, rel_std_betaRef*mean(beta_mol(hRefIndx)), MC_count(4), 'norm'));
85+
angstroemSample = transpose(sigGenWithNoise(angstroem, sigma_angstroem, MC_count(1), 'norm'));
86+
ext_aer_sample = transpose(sigGenWithNoise(ext_aer, sigma_ext_aer, MC_count(2), 'norm'));
87+
sigElasticSample = transpose(sigGenWithNoise(sigElastic, sqrt(sigElastic + bgElastic), MC_count(3), 'norm'));
88+
sigVRN2Sample = transpose(sigGenWithNoise(sigVRN2, sqrt(sigVRN2 + bgVRN2), MC_count(3), 'norm'));
89+
90+
aerBscSample = NaN(prod(MC_count), length(ext_aer));
91+
for iX = 1:MC_count(1)
92+
for iY = 1:MC_count(2)
93+
for iZ = 1:MC_count(3)
94+
for iM = 1:MC_count(4)
95+
aerBscSample(iM + MC_count(4)*(iZ - 1) + MC_count(4)*MC_count(3)*(iY - 1) + MC_count(4)*MC_count(3)*MC_count(2)*(iX - 1), :) = ...
96+
pollyRamanBsc_smart(height, sigElasticSample(iZ, :), sigVRN2Sample(iZ, :), ext_aer_sample(iY, :), angstroemSample(iX, :), ext_mol, beta_mol,ext_mol_raman, beta_mol_inela, HRef, el_lambda, betaRefSample(iM), window_size, flagSmoothBefore, el_lambda, inel_lambda);
97+
end
98+
end
99+
end
100+
end
101+
102+
aerBscStd = nanstd(aerBscSample, 1, 0);
103+
104+
elseif strcmpi(method, 'analytical')
105+
aerBscStd=NaN(length(beta_aer));
106+
%TODO: analytical error analysis for Raman Backscatter retrieval
107+
else
108+
aerBscStd=NaN(length(beta_aer));
109+
error('Unkown method to estimate the uncertainty.');
110+
end
111+
112+
end

‎lib/retrievals/pollyRamanExt_smart.m‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
function [ ext_aer ] = pollyRamanExt_smart(height, sig, lambda_emit, ...
2+
lambda_Raman,alpha_molecular_elastic, alpha_molecular_Raman, number_density, angstrom, window_size, method, measure_error)
3+
% POLLYRAMANEXT etrieve the aerosol extinction coefficient with Raman method
4+
%
5+
% USAGE:
6+
% [ ext_aer ] = pollyRamanExt(height, sig, lambda_emit, lambda_Raman, angstrom, pressure, temperature, window_size, C, rh, method, measure_error)
7+
%
8+
% INPUTS:
9+
% height: array
10+
% height[m]
11+
% sig: array
12+
% measured raman signal. Unit: Photon Count
13+
% lambda_emit: float
14+
% the wavelength of the emitted laser beam.[nm]
15+
% lambda_Raman: float
16+
% the wavelength of raman sigal.[nm]
17+
% angstrom: float
18+
% the angstrom exponent for aerosol extinction coefficient
19+
% window_size: integer
20+
% window_size for smoothing the signal with sgolay filter.
21+
% order: integer
22+
% order of the implemented sgolay filter.
23+
% method: char
24+
% specify the method to calculate the slope of the signal. You can
25+
% choose from 'moving', 'smoothing' and 'chi2'.
26+
% measure_error: array
27+
% measurement error for each bin.
28+
% alpha_molecular_elastic: array
29+
% molecular scattering coefficient at emitted wavelength in m^-1 sr^-1
30+
% alpha_molecular_Raman: array
31+
% molecular scattering coefficient at Raman wavelength in m^-1 sr^-1
32+
% number_density: array
33+
% molecular number density
34+
% Mc_countNumber of MC iterations
35+
% OUTPUTS:
36+
% ext_aer: array
37+
% aerosol extinction coefficient [m^{-1}]
38+
%
39+
% REFERENCES:
40+
% https://bitbucket.org/iannis_b/lidar_processing
41+
%
42+
% Ansmann, A. et al. Independent measurement of extinction and backscatter profiles in cirrus clouds by using a combined Raman elastic-backscatter lidar. Applied Optics Vol. 31, Issue 33, pp. 7113-7131 (1992)
43+
%
44+
% HISTORY:
45+
% - 2021-05-31: first edition by Zhenping
46+
%
47+
% .. Authors: - zhenping@tropos.de
48+
49+
% default method is movingslope
50+
if ~ exist('method', 'var')
51+
method = 'movingslope';
52+
end
53+
54+
if ~ exist('measure_error', 'var')
55+
measure_error = zeros(size(sig));
56+
end
57+
58+
temp = number_density ./ (sig .* height.^2);
59+
temp(temp <= 0) = NaN;
60+
ratio = log(temp);
61+
62+
if strcmpi(method, 'moving') || strcmpi(method, 'movingslope')
63+
deriv_ratio = movingslope_variedWin(ratio, window_size) ./ ...
64+
[height(2) - height(1), diff(height)];
65+
elseif strcmpi(method, 'smoothing') || strcmpi(method, 'smooth')
66+
deriv_ratio = movingsmooth_variedWin(ratio, window_size) ./ [height(2) - ...
67+
height(1), diff(height)];
68+
elseif strcmpi(method, 'chi2')
69+
deriv_ratio = movingLinfit_variedWin(height, ratio, measure_error, ...
70+
window_size);
71+
else
72+
error('Please set a valid method for calculate the extinction coefficient.');
73+
end
74+
75+
ext_aer = (deriv_ratio - alpha_molecular_elastic - alpha_molecular_Raman) ./ ...
76+
(1 + (lambda_emit ./ lambda_Raman) .^ angstrom);
77+
end

0 commit comments

Comments
 (0)