Skip to content
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

Record and Use Theta levels for the plot function and a few enhancements #16

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Record and Use Theta levels for the plot function and a few enhancements #16

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
69d0e0d
Create create_la_wmo_boxes_file.m
delphinedobler Aug 1, 2024
5bb39d0
DD (2024-08-01) 1
delphinedobler Aug 1, 2024
74cbc0e
DD (2024-08-01) - 2
delphinedobler Aug 1, 2024
e47e80d
DD (2024-07-02): 3.1, 3.2 and 3.3 Performance improvements
delphinedobler Aug 2, 2024
2fc4069
DD (2024-08-02) Record end of exec date and time
delphinedobler Aug 2, 2024
0c9d48d
Record and use of recorded theta levels (case of multiple calibration…
delphinedobler Sep 5, 2024
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
54 changes: 46 additions & 8 deletions matlab_codes/calculate_piecewisefit.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur

%
% Cecile Cabanes, June 2013 : calculate off-diagonal terms for error estimate: add horizontal covariance to track changes: see "change config 129"
% Delphine Dobler (DD), September 2024 :
% 4 - save computed theta levels in a file

%pn_float_dir='testfloats/';
%pn_float_name='robbins4900178';
Expand Down Expand Up @@ -156,10 +158,15 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur
return
end

% DD (2024/09-4): initialisation of arrays to be saved in a file
Theta_nseq=NaN*zeros(10,n_seq);
Index=NaN*zeros(10,n);
Plevels=NaN*zeros(10,n_seq);

% loop through sequences of calseries, ie if time series is broken into segments --
for i=1:n_seq
for iseq=1:n_seq

calindex = find(calseries==unique_cal(i));
calindex = find(calseries==unique_cal(iseq));
k = length(calindex);

% choose 10 float theta levels to use in the piecewise linear fit --------
Expand All @@ -180,6 +187,31 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur

[Theta, P, index, var_s_th, th] =...
find_10thetas( unique_SAL, unique_PTMP, unique_PRES, unique_la_ptmp, use_theta_gt, use_theta_lt, use_pres_gt, use_pres_lt, use_percent_gt);

% DD (2024/09-4): Save Theta-related information for plotting functions
Index(:,calindex)=index;
Plevels(:,iseq)=P;
[mvth,~]=size(var_s_th);
Theta_nseq(:,iseq)=Theta;
if iseq==1
Var_s_Thetas=NaN.*ones(mvth,n_seq);
Thetas=NaN.*ones(mvth,n_seq);
Var_s_Thetas(:,iseq)=var_s_th;
Thetas(:,iseq)=th;
else
[mVth,~]=size(Var_s_Thetas);
tmp1=Var_s_Thetas;
tmp2=Thetas;
mm=max(mVth,mvth);
Var_s_Thetas=NaN.*ones(mm,n_seq);
Thetas=NaN.*ones(mm,n_seq);
Var_s_Thetas(1:mVth,:)=tmp1;
Var_s_Thetas(1:mvth,iseq)=var_s_th;
Thetas(1:mVth,:)=tmp2;
Thetas(1:mvth,iseq)=th;
end
%


pp=find(isnan(index)==0);
if(isempty(pp)==0) % only proceed when there are valid levels ----
Expand Down Expand Up @@ -229,18 +261,18 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur
if isempty(breaks)
[xfit(calindex), pcond_factor(calindex), pcond_factor_err(calindex), time_deriv(calindex), ...
time_deriv_err(calindex), sta_mean(calindex), sta_rms(calindex), NDF, fitcoef, fitbreaks] = ...
fit_cond(x, y, err, covariance, 'max_no_breaks', max_breaks(i));
fit_cond(x, y, err, covariance, 'max_no_breaks', max_breaks(iseq));
else
breaks_in = breaks(i,:);
breaks_in = breaks(iseq,:);
breaks_in = breaks_in(find(isfinite(breaks_in)));
if isempty(max_breaks(i))
if isempty(max_breaks(iseq))
[xfit(calindex), pcond_factor(calindex), pcond_factor_err(calindex), time_deriv(calindex), ...
time_deriv_err(calindex), sta_mean(calindex), sta_rms(calindex), NDF, fitcoef, fitbreaks] = ...
fit_cond(x, y, err, covariance, 'breaks', breaks_in);
else
[xfit(calindex), pcond_factor(calindex), pcond_factor_err(calindex), time_deriv(calindex), ...
time_deriv_err(calindex), sta_mean(calindex), sta_rms(calindex), NDF, fitcoef, fitbreaks] = ...
fit_cond(x, y, err, covariance, 'breaks', breaks_in, 'max_no_breaks', max_breaks(i));
fit_cond(x, y, err, covariance, 'breaks', breaks_in, 'max_no_breaks', max_breaks(iseq));
end
end

Expand All @@ -266,16 +298,22 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur
sta_SAL1(:,calindex) = sw_salt( (sta_COND(:,calindex)+sta_COND_err(:,calindex))/sw_c3515, unique_PTMP, 0 );
sta_SAL_err(:,calindex) = abs(sta_SAL(:,calindex)-sta_SAL1(:,calindex));

fcoef(i,1:length(fitcoef)) = fitcoef;
fcoef(iseq,1:length(fitcoef)) = fitcoef;
if ~isempty(fitbreaks)
fbreaks(i,1:length(fitbreaks)) = fitbreaks;
fbreaks(iseq,1:length(fitbreaks)) = fitbreaks;
end

end

end %if there are valid levels
end %for each unique_cal

% DD (2024/09-4) : save theta levels information inside a file.
Theta=Theta_nseq;
ls_theta = strcat( po_system_configuration.FLOAT_CALIB_DIRECTORY, ...
pn_float_dir, 'selected_theta_', pn_float_name,'.mat');
save(ls_theta, 'Theta','Index','Plevels','Var_s_Thetas','Thetas')
% end of DD (2024/09-4)

% save calibration data --------------------------------

Expand Down
111 changes: 111 additions & 0 deletions matlab_codes/create_la_wmo_boxes_file.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
function create_la_wmo_boxes_file( po_system_configuration )

% function create_la_wmo_boxes_file( po_system_configuration )
%
% Delphine Dobler, August 2024

clim_ctd_=strcat(po_system_configuration.HISTORICAL_DIRECTORY,po_system_configuration.HISTORICAL_CTD_PREFIX);
clim_argo_=strcat(po_system_configuration.HISTORICAL_DIRECTORY,po_system_configuration.HISTORICAL_ARGO_PREFIX);
clim_bottle_=strcat(po_system_configuration.HISTORICAL_DIRECTORY,po_system_configuration.HISTORICAL_BOTTLE_PREFIX);
la_wmo_boxes_file=strcat( po_system_configuration.CONFIG_DIRECTORY, po_system_configuration.CONFIG_WMO_BOXES);

la_wmo_boxes=zeros(648,4);
la_wmo_boxes(:,1)=[1800;1700;1600;1500;1400;1300;1200;1100;1000;3000;3100;...
3200;3300;3400;3500;3600;3700;3800;1801;1701;1601;1501;...
1401;1301;1201;1101;1001;3001;3101;3201;3301;3401;3501;...
3601;3701;3801;1802;1702;1602;1502;1402;1302;1202;1102;...
1002;3002;3102;3202;3302;3402;3502;3602;3702;3802;1803;...
1703;1603;1503;1403;1303;1203;1103;1003;3003;3103;3203;...
3303;3403;3503;3603;3703;3803;1804;1704;1604;1504;1404;...
1304;1204;1104;1004;3004;3104;3204;3304;3404;3504;3604;...
3704;3804;1805;1705;1605;1505;1405;1305;1205;1105;1005;...
3005;3105;3205;3305;3405;3505;3605;3705;3805;1806;1706;...
1606;1506;1406;1306;1206;1106;1006;3006;3106;3206;3306;...
3406;3506;3606;3706;3806;1807;1707;1607;1507;1407;1307;...
1207;1107;1007;3007;3107;3207;3307;3407;3507;3607;3707;...
3807;1808;1708;1608;1508;1408;1308;1208;1108;1008;3008;...
3108;3208;3308;3408;3508;3608;3708;3808;1809;1709;1609;...
1509;1409;1309;1209;1109;1009;3009;3109;3209;3309;3409;...
3509;3609;3709;3809;1810;1710;1610;1510;1410;1310;1210;...
1110;1010;3010;3110;3210;3310;3410;3510;3610;3710;3810;...
1811;1711;1611;1511;1411;1311;1211;1111;1011;3011;3111;...
3211;3311;3411;3511;3611;3711;3811;1812;1712;1612;1512;...
1412;1312;1212;1112;1012;3012;3112;3212;3312;3412;3512;...
3612;3712;3812;1813;1713;1613;1513;1413;1313;1213;1113;...
1013;3013;3113;3213;3313;3413;3513;3613;3713;3813;1814;...
1714;1614;1514;1414;1314;1214;1114;1014;3014;3114;3214;...
3314;3414;3514;3614;3714;3814;1815;1715;1615;1515;1415;...
1315;1215;1115;1015;3015;3115;3215;3315;3415;3515;3615;...
3715;3815;1816;1716;1616;1516;1416;1316;1216;1116;1016;...
3016;3116;3216;3316;3416;3516;3616;3716;3816;1817;1717;...
1617;1517;1417;1317;1217;1117;1017;3017;3117;3217;3317;...
3417;3517;3617;3717;3817;7817;7717;7617;7517;7417;7317;...
7217;7117;7017;5017;5117;5217;5317;5417;5517;5617;5717;...
5817;7816;7716;7616;7516;7416;7316;7216;7116;7016;5016;...
5116;5216;5316;5416;5516;5616;5716;5816;7815;7715;7615;...
7515;7415;7315;7215;7115;7015;5015;5115;5215;5315;5415;...
5515;5615;5715;5815;7814;7714;7614;7514;7414;7314;7214;...
7114;7014;5014;5114;5214;5314;5414;5514;5614;5714;5814;...
7813;7713;7613;7513;7413;7313;7213;7113;7013;5013;5113;...
5213;5313;5413;5513;5613;5713;5813;7812;7712;7612;7512;...
7412;7312;7212;7112;7012;5012;5112;5212;5312;5412;5512;...
5612;5712;5812;7811;7711;7611;7511;7411;7311;7211;7111;...
7011;5011;5111;5211;5311;5411;5511;5611;5711;5811;7810;...
7710;7610;7510;7410;7310;7210;7110;7010;5010;5110;5210;...
5310;5410;5510;5610;5710;5810;7809;7709;7609;7509;7409;...
7309;7209;7109;7009;5009;5109;5209;5309;5409;5509;5609;...
5709;5809;7808;7708;7608;7508;7408;7308;7208;7108;7008;...
5008;5108;5208;5308;5408;5508;5608;5708;5808;7807;7707;...
7607;7507;7407;7307;7207;7107;7007;5007;5107;5207;5307;...
5407;5507;5607;5707;5807;7806;7706;7606;7506;7406;7306;...
7206;7106;7006;5006;5106;5206;5306;5406;5506;5606;5706;...
5806;7805;7705;7605;7505;7405;7305;7205;7105;7005;5005;...
5105;5205;5305;5405;5505;5605;5705;5805;7804;7704;7604;...
7504;7404;7304;7204;7104;7004;5004;5104;5204;5304;5404;...
5504;5604;5704;5804;7803;7703;7603;7503;7403;7303;7203;...
7103;7003;5003;5103;5203;5303;5403;5503;5603;5703;5803;...
7802;7702;7602;7502;7402;7302;7202;7102;7002;5002;5102;...
5202;5302;5402;5502;5602;5702;5802;7801;7701;7601;7501;...
7401;7301;7201;7101;7001;5001;5101;5201;5301;5401;5501;...
5601;5701;5801;7800;7700;7600;7500;7400;7300;7200;7100;...
7000;5000;5100;5200;5300;5400;5500;5600;5700;5800];

for i_hist_type=2:4
if i_hist_type==2
clim_dir=clim_ctd_;
clim_type="ctd_";
prefix=po_system_configuration.HISTORICAL_CTD_PREFIX;
end
if i_hist_type==4
clim_dir=clim_argo_;
clim_type="argo_";
prefix=po_system_configuration.HISTORICAL_ARGO_PREFIX;
end
if i_hist_type==3
clim_dir=clim_bottle_;
clim_type="bot_";
prefix=po_system_configuration.HISTORICAL_BOTTLE_PREFIX;
end

prefix_str=split(prefix,"/");
clim_dir=strcat(po_system_configuration.HISTORICAL_DIRECTORY,strjoin(prefix_str(1:end-1),"/"));

listing=dir(clim_dir);
n_files=size(listing,1);
for ifile=1:n_files
l_name=listing(ifile).name;
if strlength(l_name) < strlength(clim_type)
continue
end
if l_name(1:strlength(clim_type))== clim_type
box_no=strrep(l_name,clim_type,'');
box_no=strrep(box_no,'.mat','');
box_no=str2double(box_no);
i_box=(la_wmo_boxes(:,1)==box_no);
la_wmo_boxes(i_box,i_hist_type)=1;
end

end
end

save(la_wmo_boxes_file,'la_wmo_boxes')
30 changes: 18 additions & 12 deletions matlab_codes/get_region_ow.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@

function [ pa_grid_lat, pa_grid_long, pa_grid_dates ] = get_region_ow(pa_wmo_numbers, pa_float_name, po_config_data) ;
function [ pa_grid_lat, pa_grid_long, pa_grid_dates, pa_grid_box_id ] = get_region_ow(pa_wmo_numbers, pa_float_name, po_config_data)

% function [ pa_grid_lat, pa_grid_long, pa_grid_dates ] = get_region_ow( pa_wmo_numbers, pa_float_name, po_config_data ) ;
%
Expand All @@ -18,12 +18,16 @@
% Annie Wong, December 2007
% Breck Owens, December 2006
% C.Cabanes (2015) to save time: load only long, lat, dates and source data
% D.Dobler (DD), August 2024: to save more time later in retr_region_ow: also output
% corresponding boxes ids. Additional minor modifications: comment date_hist as it is unused, and
% correct indentation twice.

pa_grid_lat = [ ] ;
pa_grid_long = [ ] ;
pa_grid_dates = [ ] ;
pa_grid_lat = [ ] ;
pa_grid_long = [ ] ;
pa_grid_dates = [ ] ;
pa_grid_box_id = [ ] ;

[m,n]=size(pa_wmo_numbers);
[m,~]=size(pa_wmo_numbers);



Expand All @@ -40,7 +44,7 @@
%toc

not_use=[];
date_hist = changedates(lo_box_data.dates);
%date_hist = changedates(lo_box_data.dates);
lo_box_data.lat(not_use)=[];
lo_box_data.long(not_use)=[];
lo_box_data.dates(not_use)=[];
Expand All @@ -55,22 +59,24 @@
not_use=[];
for i=1:length(lo_box_data.lat)
profile=lo_box_data.source{i};
jj=findstr(profile,'_');
jj=findstr(profile,'_');
ref_float=profile(1:jj-1);
kk=findstr(pa_float_name, ref_float);
if(isempty(kk)==0)
not_use=[not_use,i];
not_use=[not_use,i];
end
end
lo_box_data.lat(not_use)=[];
lo_box_data.long(not_use)=[];
lo_box_data.dates(not_use)=[];
%-----------------------------------------------------------------------
end

pa_grid_lat = [ pa_grid_lat, lo_box_data.lat ] ;
pa_grid_long = [ pa_grid_long, lo_box_data.long ] ;
pa_grid_dates = [ pa_grid_dates, lo_box_data.dates ] ;

n_observations = size(lo_box_data.lat,2);
pa_grid_box_id = [ pa_grid_box_id, pa_wmo_numbers(ln_index,1)*ones(1,n_observations)] ;
pa_grid_lat = [ pa_grid_lat, lo_box_data.lat ] ;
pa_grid_long = [ pa_grid_long, lo_box_data.long ] ;
pa_grid_dates = [ pa_grid_dates, lo_box_data.dates ] ;
end
end
fclose('all');
Expand Down
11 changes: 9 additions & 2 deletions matlab_codes/map_data_grid.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,10 @@
% Bretherton, etal, 1975) is removed. The error estimate includes
% the contributions from both the mapping and the mean
% estimatation.
%
% Delphine Dobler (DD), August 2024:
% 3.3 - Performance: Inside map_data_grid, comment calculation
% of vdataerror as this is never used afterwards.
% **********************************************************************

[m,n]=size(v);
Expand Down Expand Up @@ -72,8 +76,11 @@
Cmd = s*covarxyt_pv(posdata, posdata, gs, ts, q, phi, pv);
vdata = Cmd*wght + vm;
% include error in mean (eqn 24 of Bretherton, etal)
vdataerror = repmat(sqrt(s-diag(Cmd*Cdd*Cmd') ...
+(1- sum(Cmd*Cdd,2)).^2/sum_Cdd),1,1);
%vdataerror = repmat(sqrt(s-diag(Cmd*Cdd*Cmd') ...
% +(1- sum(Cmd*Cdd,2)).^2/sum_Cdd),1,1);
% DD (2024/08-3.3): never used afterwards, save some time (Uncomment to
% compute)
vdataerror = [];

% map to posgrid -----

Expand Down
Loading