PollyNET / Pollynet_Processing_Chain

NRT lidar data processing program for multiwavelength polarization Raman lidar network (PollyNET)
https://polly.tropos.de/
GNU General Public License v3.0
19 stars 8 forks source link

Overlap calculation: Raman method missing #204

Open martin-rdz opened 1 year ago

martin-rdz commented 1 year ago

When lookin into the arielle processing, @HolgerPollyNet and me realized, that one (documented) overlap calculation method is missing. See Code of overlapCalc

The documentation reads

%    overlapCalMode: integer
%        overlap calculation mode.
%        1: signal ratio between near-range and far-range channel
%        2: Raman method

Additionally, that contradictory to the described behaviour of pollyOVLCalc

%    overlapCalMode: numeric
%        overlap calculation mode. (default: 1)
%        0: no overlap correction
%        1:overlap correction with using the default overlap function
%        2: overlap correction with using the calculated overlap function
%        3: overlap correction with gluing near-range and far-range signal

Is there an easy fix, allowing rapid implementation or would that need extended development? @HolgerPollyNet again, what overlap estimates are used in Verlauf?

Best, Martin

HolgerPollyNet commented 1 year ago

@cristoferjimenez Do you have already a Matlab code for this issue?

HolgerPollyNet commented 1 year ago

We got an old code from Birgit, is somebody volunteering to implement this? @AthenaAugousta @gast-ben @cristoferjimenez @Moritz-TROPOS @ulysses78 @Ronny-TROPOS @griesche

% overlap.m % % nach Wandinger 2002, Experimental determination of the lidar % overlap profile with Raman lidar, AO, Vol. 41, no 3, Jan 2002 % %
% 12/06 BHeese
% % overlapp = p_R(z)z2 / (Cbeta_R(z)T_0(z) * T_R(z))

clear p_ave_elast_ana p_ave_raman_ana clear m_ave_elast m_ave_raman clear T_0 T_r

ray_fak_387 = (355/387)^4.085; beta_mol_387= beta_molray_fak_387; ray_fak_532 = (355/532)^4.085; beta_mol_532 = beta_molray_fak_532;

%plot(beta_mol) %hold on %plot(beta_mol_387_n,'r') %plot(beta_mol_532,'g') %plot(beta_mol_532_n,'g') % RefBin = RefBin_ana1; deltar = 7.5e-3; % T_0 = ones(RefBin,1); T_r = ones(RefBin,1); % Const = 25.3; % am 27.05.06 % Const = 15.4; %am 15.5.06 % Const = 38; % for i = 21:RefBin
% Raman Particle extinction p_ave_raman_ana(i) = 0.5(real(aero_ext_raman_ana(i)) + ... real(aero_ext_raman_ana(i+1))); % Raman molecular extinction m_ave_raman(i) = 0.5(ray_ext_387(i)+ray_ext_387(i+1));
% Elastic particle extinction p_ave_elast_ana(i) = 0.5(real(aero_ext_raman_ana(i)) + real(aero_ext_raman_ana(i+1))); % Elastic molecular extinction m_ave_elast(i) = 0.5(ray_ext_355(i)+ray_ext_355(i+1)); % % Transmissionen T_0(i+1) = T_0(i) exp(-(p_ave_elast_ana(i) + m_ave_elast(i))deltar); T_r(i+1) = T_r(i) exp(-(p_ave_raman_ana(i) + m_ave_raman(i))deltar);

% T_0(i) = exp(-(ray_ext_355 + aero_ext_raman_PC(i))); % T_r(i) = exp(-(ray_ext_387(i) + aero_ext_raman_PC(i)(355/387))); %
overl(i) = pr2_ana2_int(i)/(Const
beta_raman_ana(i)T_0(i)T_r(i)); end % overl_sm = smooth(overl(:),21,'sgolay',2);

% axis([0 alt(500) 0 1.1 ]) figure(1) set(gcf,'position',[scrsz(3)-0.95scrsz(3) scrsz(4)-0.95scrsz(4) scrsz(3)-0.6scrsz(3) scrsz(4)-0.15scrsz(4)]);
% subplot(1,2,1) subplot('Position',[0.1 0.08 0.37 0.85]); title(['Overlap function on ' datum],'fontsize',[9]) ylabel('range / km','fontsize',[12]) axis([-.05 1.5 0 range(500)]); box on hold on plot(overl(1:500),range(1:500)) hold on x=ones(500); plot(x, range(1:500),'r') legend('overlap', 'unity') grid on

subplot(1,2,2) subplot('Position',[0.55 0.08 0.37 0.85]); title([' ' filename(1,8:9) ':' filename(1,10:11) ' - ' ... filename(nfiles,8:9) ':' filename(nfiles,10:11) ' UTC'],'fontsize',[9]) axis([0 10 0 range(500)]); box on hold on plot(pr2_ana1_int(1:500)./overl(1:500)',range(1:500)) hold on plot(pr2_ana1_int(1:500),range(1:500),'r') legend('RC-Signal/overlap', 'RC-Signal') grid on

cristoferjimenez commented 1 year ago

I will try to implement this module into the Pollynet chain during the next days. From Birgit's code it seems it should go straigth forward

HolgerPollyNet commented 1 year ago

@cristoferjimenez Your code is not yet working on the server (example cabo verde data):

[2023-05-31 13:23:48] Start retrieving aerosol optical properties. Meteorological file : . [2023-05-31 13:23:48] Start overlap estimation. Attempted to access aerBsc355_raman_mean(5); index out of bounds because numel(aerBsc355_raman_mean)=1.

Error in picassoProcV3 (line 1953) aerBsc355_raman_mean(1:5)=aerBsc355_raman_mean(5); %replace first 5 bins by a constant value

Error in picassoProcTodolist (line 62) reportTmp = picassoProcV3(pollyDataFile, pollyDataTasks.pollyType{iTask}, ...

Error in picassoProcHistoryData (line 69) picassoProcTodolist(PicassoConfigFile);

cristoferjimenez commented 1 year ago

I thought it won't be a problem since aerBsc355_raman should be at least a NaN-matrix. But of course., if clFreGrps is empty the aerBsc355_raman will be of size zeroxheight_length... I will consider this for the second version.

I think for now you can just comment the line or put a condition...

Like: if size(aerBsc355_raman,1)>0 aerBsc355_raman_mean(1:5)=aerBsc355_raman_mean(5); %replace first 5 bins by a constant value end

HolgerPollyNet commented 1 year ago

Hi Cristofer, I think it would be a good exercise if you can do the testing and final implementation by your own ( I will start for now, but I guess the issue will not be fully solved yet). @ulysses78 can he work directly on the server in his home directory, maybe giving his user full read rights?

ulysses78 commented 1 year ago

But Cristofer has to have a matlab-license on the server if he wants to use picasso on the rsd2-server. Maybe a good solution would be to work as user Bildermacher2 and create your own experimental polly-directory.

HolgerPollyNet commented 1 year ago

I thought it won't be a problem since aerBsc355_raman should be at least a NaN-matrix. But of course., if clFreGrps is empty the aerBsc355_raman will be of size zeroxheight_length... I will consider this for the second version.

I think for now you can just comment the line or put a condition...

Like: if size(aerBsc355_raman,1)>0 aerBsc355_raman_mean(1:5)=aerBsc355_raman_mean(5); %replace first 5 bins by a constant value end

Still issues with this solution:

[2023-06-02 07:16:03] Start overlap estimation. Error using .* Matrix dimensions must agree.

Error in picassoProcV3 (line 1962) trans387 = exp(-cumsum((mExt387+PollyConfig.LR355aerBsc355_raman_mean(355/387)^PollyConfig.angstrexp) .* [data.distance0(1), diff(data.distance0)]));