diff --git a/matlab_codes/calculate_piecewisefit.m b/matlab_codes/calculate_piecewisefit.m index b62e258..0ce3f65 100755 --- a/matlab_codes/calculate_piecewisefit.m +++ b/matlab_codes/calculate_piecewisefit.m @@ -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'; @@ -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 -------- @@ -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 ---- @@ -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 @@ -266,9 +298,9 @@ 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 @@ -276,6 +308,12 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur 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 -------------------------------- diff --git a/matlab_codes/create_la_wmo_boxes_file.m b/matlab_codes/create_la_wmo_boxes_file.m new file mode 100644 index 0000000..c56922d --- /dev/null +++ b/matlab_codes/create_la_wmo_boxes_file.m @@ -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') diff --git a/matlab_codes/get_region_ow.m b/matlab_codes/get_region_ow.m index cb021cc..ab9fc4c 100755 --- a/matlab_codes/get_region_ow.m +++ b/matlab_codes/get_region_ow.m @@ -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 ) ; % @@ -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); @@ -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)=[]; @@ -55,11 +59,11 @@ 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)=[]; @@ -67,10 +71,12 @@ 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'); diff --git a/matlab_codes/map_data_grid.m b/matlab_codes/map_data_grid.m index 12e048f..30bf750 100755 --- a/matlab_codes/map_data_grid.m +++ b/matlab_codes/map_data_grid.m @@ -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); @@ -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 ----- diff --git a/matlab_codes/plot_diagnostics_ow.m b/matlab_codes/plot_diagnostics_ow.m index d945a36..638907d 100644 --- a/matlab_codes/plot_diagnostics_ow.m +++ b/matlab_codes/plot_diagnostics_ow.m @@ -4,12 +4,43 @@ function plot_diagnostics_ow( pn_float_dir, pn_float_name, po_system_configurati % % Annie Wong, 14 June 2011 % Breck Owens, October 2006 +% +% Delphine Dobler (DD), August 2024: +% 2 - create output directory if they do not exist already +% Delphine Dobler (DD), September 2024: +% 4.1 - Use of theta levels from those now saved during +% calculate_piecewisefit function call (also corrects graph 2 +% and 4 in case of several calibration series). +% 4.2 - For graph 6 and 8: Change indexes names by a more +% explicit terminology: usage of iseq, istep, ilevel, iplot instead of i,j,k ... +% Remove unused parameter lines. +% These changes were not marked done for the sake of lisibility +% 4.3 - Add a configuration parameter to +% select from eps or png type of graph format. +% 4.4 - Enhance Graph 3: curate legend content, split title in two +% lines, add a grid, position the legend outside the graph. +% 4.5 - Automalically set graph 5 y-axis limits from selected thetas +% and set the upper bound to 14°C. + %-------------------------------------------------------------------------- %pn_float_dir='uw/'; %pn_float_name='R5902134'; %po_system_configuration = load_configuration( 'ow_config.txt' ); +outdir=[ po_system_configuration.FLOAT_PLOTS_DIRECTORY pn_float_dir ]; +if not(isfolder(outdir)) + mkdir(outdir) +end + +% DD (2024/09 - 4.3) : Configuration parameter to select graph format +graph_format = 'png'; % 'png' or 'eps' +if strcmp(graph_format,'png') + graph_print_option = '-dpng'; +end +if strcmp(graph_format,'eps') + graph_print_option = '-depsc'; +end close all @@ -65,6 +96,17 @@ function plot_diagnostics_ow( pn_float_dir, pn_float_name, po_system_configurati use_pres_lt = lo_float_calseries.use_pres_lt; use_percent_gt = lo_float_calseries.use_percent_gt; +% DD (2024/09 - 4.1) : load the saved theta levels: +ls_theta = strcat( po_system_configuration.FLOAT_CALIB_DIRECTORY, ... + pn_float_dir, 'selected_theta_', pn_float_name,'.mat'); +theta_levels=load(ls_theta); +tlevels=theta_levels.Theta; +plevels=theta_levels.Plevels; +index=theta_levels.Index; +index(index==0)=NaN; +var_s_Thetalevels=theta_levels.Var_s_Thetas; +Thetalevels=theta_levels.Thetas; + % load the station by station fits load(fullfile( po_system_configuration.FLOAT_CALIB_DIRECTORY, pn_float_dir,... strcat( po_system_configuration.FLOAT_CALIB_PREFIX, pn_float_name, po_system_configuration.FLOAT_CALIB_POSTFIX ) ),'-regexp','^sta') @@ -117,7 +159,7 @@ function plot_diagnostics_ow( pn_float_dir, pn_float_name, po_system_configurati plot(LONG,LAT,'r-'); hold on plot(x,y,'b.') -legend('float','historical points','Location','Best') +%legend('float','historical points','Location','Best') plot(coastdata_x,coastdata_y,'k.-'); for i=1:n @@ -152,7 +194,7 @@ function plot_diagnostics_ow( pn_float_dir, pn_float_name, po_system_configurati drawnow set(gcf,'papertype','usletter','paperunits','inches','paperorientation','portrait','paperposition',[.25,.75,8,9.5]); -print('-depsc', strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_1.eps')); +print(graph_print_option, strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_1.',graph_format)); % plot the uncalibrated theta-S curves from the float (figure 2) -------- @@ -185,7 +227,9 @@ function plot_diagnostics_ow( pn_float_dir, pn_float_name, po_system_configurati set(qq,'color',c(i,:)) ; end -[tlevels, plevels, index, var_s_Thetalevels, Thetalevels] = find_10thetas( SAL, PTMP, PRES, la_ptmp, use_theta_gt, use_theta_lt, use_pres_gt, use_pres_lt, use_percent_gt); +% DD (2024/09 - 4.1) : correction when several calibration series: use of saved +% theta levels from calculate_piecewisefit function (loaded herebefore in the loading steps). +%[tlevels, plevels, index, var_s_Thetalevels, Thetalevels] = find_10thetas( SAL, PTMP, PRES, la_ptmp, use_theta_gt, use_theta_lt, use_pres_gt, use_pres_lt, use_percent_gt); for i=1:n b = find( isnan(index(:,i))==0 ); @@ -209,7 +253,7 @@ function plot_diagnostics_ow( pn_float_dir, pn_float_name, po_system_configurati drawnow title( strcat( title_floatname, ' uncalibrated float data (-) and mapped salinity (o) with objective errors' ) ); set(gcf,'papertype','usletter','paperunits','inches','paperorientation','portrait','paperposition',[.25,.75,8,9.5]); -print('-depsc', strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_2.eps')); +print(graph_print_option, strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_2.',graph_format)); % calibration curve (figure 3) -------------------------- @@ -248,45 +292,48 @@ function plot_diagnostics_ow( pn_float_dir, pn_float_name, po_system_configurati set(gcf,'defaultlinelinewidth',2) set(gcf,'defaultaxesfontsize',16) subplot(2,1,1) -plot(PROFILE_NO, pcond_factor, 'b-'); hold on -plot(PROFILE_NO, pcond_factor, 'g-'); -% plot station by station fit +% DD (2024/09 - 4.4): correction on legend display +p1=errorbar(PROFILE_NO, pcond_factor, 2*pcond_factor_err,'b'); +p2=errorbar(PROFILE_NO, pcond_factor, pcond_factor_err,'g*-'); ok = find(isfinite(sta_mean)); -plot(PROFILE_NO(ok), sta_mean(ok), 'r-'); -legend('2 x cal error','1 x cal error','1-1 profile fit', 'Location', 'Best'); -errorbar(PROFILE_NO, pcond_factor, 2*pcond_factor_err,'b') -errorbar(PROFILE_NO, pcond_factor, pcond_factor_err,'g*-') -plot(PROFILE_NO(ok), sta_mean(ok), 'r-'); - +p3=plot(PROFILE_NO(ok), sta_mean(ok), 'r-'); plot( [0, max(PROFILE_NO)+1], [1,1], 'k-') +legend([p1 p2 p3],{'2 x cal error','1 x cal error','1-1 profile fit'}, 'Location', 'bestoutside'); +grid on + axis([ 0, max(PROFILE_NO)+1, min([pcond_factor-pcond_factor_err,1])-.0004, max([pcond_factor+pcond_factor_err,1])+.0004 ]) + set(gca,'FontSize',12) ylabel('r') % multiplicative term has no units -title( strcat(title_floatname, ' potential conductivity (mmho/cm) multiplicative correction r with errors') ); +% DD (2024/09 - 4.4): split title on two lines +%title( strcat(title_floatname, ' potential conductivity (mmho/cm) multiplicative correction r with errors') ); +title({strcat(title_floatname, ' potential conductivity (mmho/cm)');'multiplicative correction r with errors'}); subplot(2,1,2) -plot(PROFILE_NO, avg_Soffset, 'b-'); hold on -plot(PROFILE_NO, avg_Soffset, 'g-'); -% Plot station by station fit + +% DD (2024/09 - 4.4): correction on legend display +p1=errorbar(PROFILE_NO, avg_Soffset, 2*avg_Soffset_err,'b'); +p2=errorbar(PROFILE_NO, avg_Soffset, avg_Soffset_err,'go-'); ok = find(isfinite(avg_Staoffset)); -plot(PROFILE_NO(ok), avg_Staoffset(ok), 'r-'); -legend('2 x cal error','1 x cal error','1-1 profile fit', 'Location', 'Best'); -errorbar(PROFILE_NO, avg_Soffset, 2*avg_Soffset_err,'b') -errorbar(PROFILE_NO, avg_Soffset, avg_Soffset_err,'go-') -plot(PROFILE_NO(ok), avg_Staoffset(ok), 'r-'); +p3=plot(PROFILE_NO(ok), avg_Staoffset(ok), 'r-'); +plot( [0, max(PROFILE_NO)+1], [0,0], 'k-') +legend([p1 p2 p3],'2 x cal error','1 x cal error','1-1 profile fit', 'Location', 'bestoutside'); +grid on axis([ 0, max(PROFILE_NO)+1, min([avg_Soffset-avg_Soffset_err,0])-.02, max([avg_Soffset+avg_Soffset_err,0])+.02 ]) -plot( [0, max(PROFILE_NO)+1], [0,0], 'k-') + set(gca,'FontSize',12) xlabel('float profile number'); -ylabel('\Delta S') -title( strcat(title_floatname, ' vertically-averaged salinity (PSS-78) additive correction \Delta S with errors') ); +ylabel('\Delta S (PSS-78)') +% DD (2024/09 - 4.4): split title on two lines +%title( strcat(title_floatname, ' vertically-averaged salinity (PSS-78) additive correction \Delta S with errors') ); +title({strcat(title_floatname, ' vertically-averaged salinity (PSS-78)');'additive correction \Delta S with errors'}); drawnow set(gcf,'papertype','usletter','paperunits','inches','paperorientation','portrait','paperposition',[.25,.75,8,9.5]); -print('-depsc', strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_3.eps')); +print(graph_print_option, strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_3.',graph_format)); % plot the calibrated theta-S curves from the float (figure 4) -------------------------- @@ -343,7 +390,7 @@ function plot_diagnostics_ow( pn_float_dir, pn_float_name, po_system_configurati drawnow title( strcat(title_floatname, ' calibrated float data (-) and mapped salinity (o) with objective errors' ) ); set(gcf,'papertype','usletter','paperunits','inches','paperorientation','portrait','paperposition',[.25,.75,8,9.5]); -print('-depsc', strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_4.eps')); +print(graph_print_option, strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_4.',graph_format)); % Brian King's plot: salinity anomaly time series on theta levels (figure 5) ------------ @@ -355,7 +402,9 @@ function plot_diagnostics_ow( pn_float_dir, pn_float_name, po_system_configurati fl.useqc = '0'; fl.plot = 1; -fl.yaxes = [2 5 20]; +% DD (2024/09 - 4.5) : automatically set y-axis limits using tlevels. +%fl.yaxes = [2 5 20]; +fl.yaxes = [floor(min(tlevels,[],'all')) ceil(max(tlevels,[],'all')) ceil(max(TEMP))]; d.PSAL = SAL; d.TEMP = TEMP; d.PRES = PRES; @@ -371,59 +420,62 @@ function plot_diagnostics_ow( pn_float_dir, pn_float_name, po_system_configurati drawnow set(gcf,'papertype','usletter','paperunits','inches','paperorientation','portrait','paperposition',[.25,.5,8,10]); -print('-depsc', strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_5.eps')); +print(graph_print_option, strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_5.',graph_format)); % plot salinity time series on theta levels with the smallest S variance (figure 6) ------------ % CC changes 06/23 figure 6 to take into account the splitting of the time series -% add l381 to l404 +% DD (2024/09 - 4.1): updated to account for saved thetas instead of re-computing them (see loading step +% here above and calculate_piecewisefit function) + change i,j indexes by +% more explicit names (iplot, iseq, icycle). + unique_cal = unique(lo_float_calseries.calseries); +unique_cal(unique_cal==0)=[]; % as is done in calculate_piecewisefit n_seq = length(unique_cal); -for iy=1:n_seq - if unique_cal(iy)>0 - calindex = find(lo_float_calseries.calseries==unique_cal(iy)); - k = length(calindex); - - % choose the two most stable theta level for each part of the time series -------- - unique_SAL = SAL(:, calindex); - unique_PTMP = PTMP(:, calindex); - unique_PRES = PRES(:, calindex); - unique_la_ptmp = la_ptmp(:, calindex); - unique_mapped_sal = mapped_sal(:, calindex); - unique_mapsalerrors = mapsalerrors(:, calindex); - - ten_SAL = NaN.*ones(10,k); - ten_PTMP = NaN.*ones(10,k); - ten_PRES = NaN.*ones(10,k); - ten_mapped_sal = NaN.*ones(10,k); - ten_mapsalerrors = NaN.*ones(10,k); - - - [tlevels, plevels, index, var_s_Thetalevels, Thetalevels] = 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); - %[tlevels, plevels, index, var_s_Thetalevels, Thetalevels] = find_10thetas( SAL, PTMP, PRES, la_ptmp, use_theta_gt, use_theta_lt, use_pres_gt, use_pres_lt, use_percent_gt); CC changes 06/23 +for iseq=1:n_seq + + tplot=[1:2]; + nplot=length(tplot); + Sint=NaN.*ones(nplot,n); + Smap=NaN.*ones(nplot,n); + Smaperr=NaN.*ones(nplot,n); + Scal=NaN.*ones(nplot,n); + Scalerr=NaN.*ones(nplot,n); + Thetalevel_indexes=NaN.*ones(nplot,n); + + % CC full time series is plotted on specified theta + trimPRES=PRES; % use only manually specified THETA & PRES range --- + trimSAL=SAL; + trimPTMP=PTMP; + trim_mapped_sal=mapped_sal; + trim_mapsalerrors=mapsalerrors; + trim_cal_SAL=cal_SAL; + trim_cal_SAL_err=cal_SAL_err; + + jj=find(isnan(la_ptmp)==1); + trimPRES(jj)=NaN; + trimSAL(jj)=NaN; + trimPTMP(jj)=NaN; + trim_mapped_sal(jj)=NaN; + trim_mapsalerrors(jj)=NaN; + trim_cal_SAL(jj)=NaN; + trim_cal_SAL_err(jj)=NaN; + + if( isempty(use_theta_lt)==0 & isempty(use_theta_gt)==1 ) + jj=find(trimPTMP>use_theta_lt); + trimPRES(jj)=NaN; + trimSAL(jj)=NaN; + trimPTMP(jj)=NaN; + trim_mapped_sal(jj)=NaN; + trim_mapsalerrors(jj)=NaN; + trim_cal_SAL(jj)=NaN; + trim_cal_SAL_err(jj)=NaN; + end - - tplot=[1:2]; - p=length(tplot); - Sint=NaN.*ones(p,n); - Smap=NaN.*ones(p,n); - Smaperr=NaN.*ones(p,n); - Scal=NaN.*ones(p,n); - Scalerr=NaN.*ones(p,n); - Thetalevel_indexes=NaN.*ones(p,n); - - % CC full time series is plotted on specified theta - trimPRES=PRES; % use only manually specified THETA & PRES range --- - trimSAL=SAL; - trimPTMP=PTMP; - trim_mapped_sal=mapped_sal; - trim_mapsalerrors=mapsalerrors; - trim_cal_SAL=cal_SAL; - trim_cal_SAL_err=cal_SAL_err; - - jj=find(isnan(la_ptmp)==1); + if( isempty(use_theta_gt)==0 & isempty(use_theta_lt)==1 ) + jj=find(trimPTMPuse_theta_lt); - trimPRES(jj)=NaN; - trimSAL(jj)=NaN; - trimPTMP(jj)=NaN; - trim_mapped_sal(jj)=NaN; - trim_mapsalerrors(jj)=NaN; - trim_cal_SAL(jj)=NaN; - trim_cal_SAL_err(jj)=NaN; - end - - if( isempty(use_theta_gt)==0 & isempty(use_theta_lt)==1 ) - jj=find(trimPTMPuse_theta_lt) %the middle band is excluded - jj=find(trimPTMPuse_theta_lt); - else - jj=find(trimPTMPuse_theta_lt); - end - trimPRES(jj)=NaN; - trimSAL(jj)=NaN; - trimPTMP(jj)=NaN; - trim_mapped_sal(jj)=NaN; - trim_mapsalerrors(jj)=NaN; - trim_cal_SAL(jj)=NaN; - trim_cal_SAL_err(jj)=NaN; - end - - if( isempty(use_pres_lt)==0 & isempty(use_pres_gt)==1 ) - jj=find(trimPRES>use_pres_lt); - trimPRES(jj)=NaN; - trimSAL(jj)=NaN; - trimPTMP(jj)=NaN; - trim_mapped_sal(jj)=NaN; - trim_mapsalerrors(jj)=NaN; - trim_cal_SAL(jj)=NaN; - trim_cal_SAL_err(jj)=NaN; - end - - if( isempty(use_pres_gt)==0 & isempty(use_pres_lt)==1 ) - jj=find(trimPRESuse_theta_lt) %the middle band is excluded + jj=find(trimPTMPuse_theta_lt); + else + jj=find(trimPTMPuse_theta_lt); end - - if( isempty(use_pres_gt)==0 & isempty(use_pres_lt)==0 ) - if(use_pres_gt>use_pres_lt) %he middle band is excluded - jj=find(trimPRESuse_pres_lt); - else - jj=find(trimPRESuse_pres_lt); - end - trimPRES(jj)=NaN; - trimSAL(jj)=NaN; - trimPTMP(jj)=NaN; - trim_mapped_sal(jj)=NaN; - trim_mapsalerrors(jj)=NaN; - trim_cal_SAL(jj)=NaN; - trim_cal_SAL_err(jj)=NaN; + trimPRES(jj)=NaN; + trimSAL(jj)=NaN; + trimPTMP(jj)=NaN; + trim_mapped_sal(jj)=NaN; + trim_mapsalerrors(jj)=NaN; + trim_cal_SAL(jj)=NaN; + trim_cal_SAL_err(jj)=NaN; + end + + if( isempty(use_pres_lt)==0 & isempty(use_pres_gt)==1 ) + jj=find(trimPRES>use_pres_lt); + trimPRES(jj)=NaN; + trimSAL(jj)=NaN; + trimPTMP(jj)=NaN; + trim_mapped_sal(jj)=NaN; + trim_mapsalerrors(jj)=NaN; + trim_cal_SAL(jj)=NaN; + trim_cal_SAL_err(jj)=NaN; + end + + if( isempty(use_pres_gt)==0 & isempty(use_pres_lt)==1 ) + jj=find(trimPRESuse_pres_lt) %the middle band is excluded + jj=find(trimPRESuse_pres_lt); + else + jj=find(trimPRESuse_pres_lt); end - - for i=1:n - for j=tplot - if(tlevels(j)min(trimPTMP(:,i))) - diffTheta = abs(trimPTMP(:,i)-tlevels(j)); - if isempty(find(~isnan(diffTheta))) - Thetalevel_indexes(j,i) = NaN; - else - Thetalevel_indexes(j,i) = min(find(diffTheta==min(diffTheta))); - end + trimPRES(jj)=NaN; + trimSAL(jj)=NaN; + trimPTMP(jj)=NaN; + trim_mapped_sal(jj)=NaN; + trim_mapsalerrors(jj)=NaN; + trim_cal_SAL(jj)=NaN; + trim_cal_SAL_err(jj)=NaN; + end + + for icycle=1:n + for iplot=tplot + if(tlevels(iplot,iseq)min(trimPTMP(:,icycle))) + diffTheta = abs(trimPTMP(:,icycle)-tlevels(iplot,iseq)); + if isempty(find(~isnan(diffTheta))) + Thetalevel_indexes(iplot,icycle) = NaN; + else + Thetalevel_indexes(iplot,icycle) = min(find(diffTheta==min(diffTheta))); end end end - - for i=tplot % build the S matrix for plotting - for j=1:n - ti=Thetalevel_indexes(i,j); - if ~isnan(ti) - interval=max(ti-1,1):min(ti+1,m); %interval is one above and one below ti - a = trimPTMP(ti,j) - trimPTMP(interval, j); - if( trimPTMP(ti,j)>tlevels(i) ) - gg=find(a>0); - if( ~isempty(gg) ) - b=find(a==min(a(gg))); %find the level with min +ve diff - ki=interval(b); - else - ki=ti; - end - end - if( trimPTMP(ti,j)tlevels(iplot,iseq) ) + gg=find(a>0); + if( ~isempty(gg) ) + b=find(a==min(a(gg))); %find the level with min +ve diff + ki=interval(b); else - Smap(i,j)=trim_mapped_sal(ti,j); % interpolate if possible because that is more accurate than using closest points - Smaperr(i,j)=trim_mapsalerrors(ti,j); % interpolate if possible because that is more accurate than using closest points + ki=ti; end - if( ki~=ti&~isnan(trim_cal_SAL(ti,j))&~isnan(trim_cal_SAL(ki,j))&~isnan(trimPTMP(ti,j))&~isnan(trimPTMP(ki,j)) ) - Scal(i,j) = interp1( [trimPTMP(ti,j), trimPTMP(ki,j)], [trim_cal_SAL(ti,j), trim_cal_SAL(ki,j)], tlevels(i) ); - Scalerr(i,j) = interp1( [trimPTMP(ti,j), trimPTMP(ki,j)], [trim_cal_SAL_err(ti,j), trim_cal_SAL_err(ki,j)], tlevels(i) ); + end + if( trimPTMP(ti,icycle)0))==1 - print('-depsc', strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_6.eps')); - else - print('-depsc', strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_6_split_',num2str(unique_cal(iy)),'.eps')); + end + + figure + set(gcf,'defaultaxeslinewidth',2) + set(gcf,'defaultlinelinewidth',2) + set(gcf,'defaultaxesfontsize',16) + + for k=1:nplot + iplot=tplot(k); + subplot(nplot,1,k) + if(isempty(Sint)==0) + plot(PROFILE_NO,Sint(iplot,:),'b*-'); + hold on + plot(PROFILE_NO,Smap(iplot,:),'r'); + plot(PROFILE_NO,Scal(iplot,:),'g'); + mm=find(isfinite(Scal(iplot,:))==1); ll=PROFILE_NO; ll=ll(mm); kk=Scal(iplot,mm); nn=Scalerr(iplot,mm); + h=fill([ll,fliplr(ll)],[kk+nn,fliplr([kk-nn])],'g'); + set(h,'EdgeColor','g'); + errorbar(PROFILE_NO,Smap(iplot,:),Smaperr(iplot,:),'r-') + plot(PROFILE_NO,Sint(iplot,:),'b*-'); + SMIN = min([Sint(iplot,:),Scal(iplot,:),Smap(iplot,:)]); + SMAX = max([Sint(iplot,:),Scal(iplot,:),Smap(iplot,:)]); + if isfinite(SMIN) & isfinite(SMAX) + axis([0, max(PROFILE_NO)+1, SMIN-.05, SMAX+.05 ]) + end + set(gca,'FontSize',12) + title( strcat(title_floatname, ' salinities with error on \theta= ', num2str(tlevels(iplot,iseq)), '^{\circ}C' ) ); + ylabel('PSS-78') end - end -end + legend('uncal float','mapped salinity','cal float w/1xerr.', 'Location', 'Best'); + xlabel('float profile number'); + + drawnow + set(gcf,'papertype','usletter','paperunits','inches','paperorientation','portrait','paperposition',[.25,.75,8,9.5]); + % CC changes 06/23 figure 6 + if length(unique_cal(unique_cal>0))==1 + print(graph_print_option, strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_6.',graph_format)); + else + print(graph_print_option, strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_6_split_',num2str(unique_cal(iseq)),'.',graph_format)); + end + +end % for iseq % Brian King's plot: salinity anomaly time series on theta levels (figure 7) ------------ @@ -621,7 +651,9 @@ function plot_diagnostics_ow( pn_float_dir, pn_float_name, po_system_configurati fl.useqc = '0'; fl.plot = 1; -fl.yaxes = [2 5 20]; +% DD (2024/09 - 4.5) : automatically set y-axis limits using tlevels. +%fl.yaxes = [2 5 20]; +fl.yaxes = [floor(min(tlevels,[],'all')) ceil(max(tlevels,[],'all')) ceil(max(TEMP))]; d.PSAL = cal_SAL; d.TEMP = TEMP; d.PRES = PRES; @@ -637,101 +669,96 @@ function plot_diagnostics_ow( pn_float_dir, pn_float_name, po_system_configurati drawnow set(gcf,'papertype','usletter','paperunits','inches','paperorientation','portrait','paperposition',[.25,.5,8,10]); -print('-depsc', strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_7.eps')); - +print(graph_print_option, strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_7.',graph_format)); % Paul Robbins' analyse variance plot (figure 8) ------------ % CC changes 06/23 figure 8 to take into account the splitting of the time series -% add l647 to l666 +% DD (2024/09 - 4.1): updated to account for saved thetas instead of re-computing them (see loading step +% here above and calculate_piecewisefit function) + change i,j indexes by +% more explicit names (iplot, iseq, icycle). + unique_cal = unique(lo_float_calseries.calseries); +unique_cal(unique_cal==0)=[]; % as is done in calculate_piecewisefit n_seq = length(unique_cal); -for iy=1:n_seq +for iseq=1:n_seq - calindex = find(lo_float_calseries.calseries==unique_cal(iy)); - k = length(calindex); - if unique_cal(iy)>0 - % choose 10 float theta levels to use in the piecewise linear fit -------- - - unique_SAL = SAL(:, calindex); - unique_cal_SAL = cal_SAL(:, calindex); - unique_PTMP = PTMP(:, calindex); - unique_PRES = PRES(:, calindex); - unique_la_ptmp = la_ptmp(:, calindex); - unique_mapped_sal = mapped_sal(:, calindex); - unique_mapsalerrors = mapsalerrors(:, calindex); - - - [tlevels, plevels, index, var_s_Thetalevels, Thetalevels] = 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); - - figure - set(gcf,'defaultaxeslinewidth',1) - set(gcf,'defaultlinelinewidth',1) - set(gcf,'defaultaxesfontsize',12) - - % plot t-s profile - subplot(222) - %plot(SAL,PTMP,'b-'); - plot(unique_SAL,unique_PTMP,'b-'); %CC changes 06/23 - ylabel('Potential temp (^{\circ}C)') - x = get(gca,'xlim'); - y = get(gca,'ylim'); - xlabel('PSS-78') - title(strcat('OW chosen levels - ', pn_float_name)); - for i=1:10 - hold on - plot(x,[tlevels(i) tlevels(i)] ,'g-'); - end - - % plot s var on t - subplot(221) - plot(var_s_Thetalevels,Thetalevels,'b-') - x = get(gca,'xlim'); - for i=1:10 - hold on - plot(x,[tlevels(i) tlevels(i)] ,'g-'); - end - %title('Salinity Variance on Theta') - title(['Salinity Variance on Theta, cycles ' num2str(PROFILE_NO(calindex(1))) '-' num2str(PROFILE_NO(calindex(end))) ]) %CC changes 06/23 - ylabel('Potential temp (^{\circ}C)') - xlabel('salinity variance') - %set(gca,'ylim',y) + calindex = find(lo_float_calseries.calseries==unique_cal(iseq)); - % plot p-t profile - subplot(223) - %plot(PTMP,-PRES,'b-'); - plot(unique_PTMP,-unique_PRES,'b-'); %CC changes 06/23 - x = get(gca,'xlim'); - xlabel('^{\circ}C') - ylabel('Pressure (dbar)') - title(strcat('OW chosen levels - ', pn_float_name)); - for i=1:10 - hold on - plot(x,[-plevels(i) -plevels(i)] ,'g-'); - end - - % plot p-s profile - subplot(224) - %plot(SAL,-PRES,'b-'); - plot(unique_SAL,-unique_PRES,'b-'); %CC changes 06/23 - x = get(gca,'xlim'); - xlabel('PSS-78') - title(strcat('OW chosen levels - ', pn_float_name)); - for i=1:10 - hold on - plot(x,[-plevels(i) -plevels(i)] ,'g-'); - end - - - drawnow - set(gcf,'papertype','usletter','paperunits','inches','paperorientation','portrait','paperposition',[.5,.25,8,10.25]); - if length(unique_cal(unique_cal>0))==1 - print('-depsc', strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_8.eps')); - else - print('-depsc', strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_8_split_', num2str(unique_cal(iy)), '.eps')); - end + % choose 10 float theta levels to use in the piecewise linear fit -------- + + unique_SAL = SAL(:, calindex); + unique_PTMP = PTMP(:, calindex); + unique_PRES = PRES(:, calindex); + + figure + set(gcf,'defaultaxeslinewidth',1) + set(gcf,'defaultlinelinewidth',1) + set(gcf,'defaultaxesfontsize',12) + + % plot t-s profile + subplot(222) + %plot(SAL,PTMP,'b-'); + plot(unique_SAL,unique_PTMP,'b-'); %CC changes 06/23 + ylabel('Potential temp (^{\circ}C)') + x = get(gca,'xlim'); + y = get(gca,'ylim'); + xlabel('PSS-78') + title(strcat('OW chosen levels - ', pn_float_name)); + for ilevel=1:10 + hold on + plot(x,[tlevels(ilevel,iseq) tlevels(ilevel,iseq)] ,'g-'); end + + % plot s var on t + subplot(221) + plot(var_s_Thetalevels(:,iseq),Thetalevels(:,iseq),'b-') + x = get(gca,'xlim'); + for ilevel=1:10 + hold on + plot(x,[tlevels(ilevel,iseq) tlevels(ilevel,iseq)] ,'g-'); + end + %title('Salinity Variance on Theta') + title(['Salinity Variance on Theta, cycles ' num2str(PROFILE_NO(calindex(1))) '-' num2str(PROFILE_NO(calindex(end))) ]) %CC changes 06/23 + ylabel('Potential temp (^{\circ}C)') + xlabel('salinity variance') + %set(gca,'ylim',y) + + % plot p-t profile + subplot(223) + %plot(PTMP,-PRES,'b-'); + plot(unique_PTMP,-unique_PRES,'b-'); %CC changes 06/23 + x = get(gca,'xlim'); + xlabel('^{\circ}C') + ylabel('Pressure (dbar)') + title(strcat('OW chosen levels - ', pn_float_name)); + for ilevel=1:10 + hold on + plot(x,[-plevels(ilevel,iseq) -plevels(ilevel,iseq)] ,'g-'); + end + + % plot p-s profile + subplot(224) + %plot(SAL,-PRES,'b-'); + plot(unique_SAL,-unique_PRES,'b-'); %CC changes 06/23 + x = get(gca,'xlim'); + xlabel('PSS-78') + title(strcat('OW chosen levels - ', pn_float_name)); + for ilevel=1:10 + hold on + plot(x,[-plevels(ilevel,iseq) -plevels(ilevel,iseq)] ,'g-'); + end + + + drawnow + set(gcf,'papertype','usletter','paperunits','inches','paperorientation','portrait','paperposition',[.5,.25,8,10.25]); + if length(unique_cal(unique_cal>0))==1 + print(graph_print_option, strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_8.',graph_format)); + else + print(graph_print_option, strcat(po_system_configuration.FLOAT_PLOTS_DIRECTORY, pn_float_dir, pn_float_name, '_8_split_', num2str(unique_cal(iseq)), '.',graph_format)); + end + end diff --git a/matlab_codes/retr_region_ow.m b/matlab_codes/retr_region_ow.m index 9391c1c..1e5df1d 100755 --- a/matlab_codes/retr_region_ow.m +++ b/matlab_codes/retr_region_ow.m @@ -1,12 +1,13 @@ -function [ pa_grid_sal, pa_grid_ptmp, pa_grid_pres, pa_grid_lat, pa_grid_long, pa_grid_dates ] = retr_region_ow( pa_wmo_numbers, pa_float_name, po_config_data, index, P, map_p_delta ) ; +function [ pa_grid_sal, pa_grid_ptmp, pa_grid_pres, pa_grid_lat, pa_grid_long, pa_grid_dates , loaded_bbox ] = ... + retr_region_ow( pa_wmo_numbers, pa_float_name, po_config_data, index, P, map_p_delta , loaded_bbox ) % % function [ pa_grid_sal, pa_grid_ptmp, pa_grid_pres, pa_grid_lat, pa_grid_long, pa_grid_dates ] = retr_region_ow( pa_wmo_numbers, pa_float_name, po_config_data, index, P, map_p_delta ) ; % -% This function retrieves historical data from the 25 selected WMO boxes, -% merges the CTD, BOT and Argo files, makes sure longitude is continuous -% around the 0-360 degree mark, and converts the dates of the +% This function retrieves historical data from the WMO boxes corresponding +% to best hist observations, merges the CTD, BOT and Argo files, makes sure +% longitude is continuous around the 0-360 degree mark, and converts the dates of the % historical data from time format2 to year.mo+ (changedates.m). % % pa_wmo_numbers can be NaN: when float profiles are out of range (65N, 65S), @@ -21,10 +22,18 @@ % modified from get_region.m, Dec 2006, Breck Owens. % % Isabelle Gaboury, 26 Sep. 2017: Added check on the dimensions of bottle data. +% +% % Delphine Dobler (DD), August 2024: +% - 3.1: correct for iteration when there is only one box to load +% (consecutive to update in the indexation for performance +% improvement); adapt the creation of output vector to the +% new index input. +% - 3.2: save loaded data, and at each following profile, +% reload only new boxes and delete unused. -zz = find(isnan(P)==0); +zz = (isnan(P)==0); max_p = max(P(zz))+map_p_delta; % max depth to retrieve historical data is deepest Argo point plus map_p_delta pa_index_0 = 0; @@ -36,95 +45,248 @@ pa_grid_long = [ ] ; pa_grid_dates = [ ] ; -[ max_depth, how_many_cols ] = size(pa_grid_pres); +%[ max_depth, how_many_cols ] = size(pa_grid_pres); + +% DD (2024/08-3.2): +% Clean up loaded_bbox for unused boxes for the treated cycle (memory limits) +list_CTD_loaded_boxes=fields(loaded_bbox.CTD); +list_ARGO_loaded_boxes=fields(loaded_bbox.ARGO); +list_BOTTLE_loaded_boxes=fields(loaded_bbox.BOTTLE); + +if size(list_CTD_loaded_boxes,1) > 1 + + list_CTD_loaded_boxes=strrep(list_CTD_loaded_boxes,'box_',''); + list_CTD_loaded_boxes=str2double(list_CTD_loaded_boxes(2:end)); + list_CTD_loaded_boxes_to_delete=list_CTD_loaded_boxes(~ismember(list_CTD_loaded_boxes,pa_wmo_numbers(:,1))); + + if size(list_CTD_loaded_boxes_to_delete,1)> 0 + for i_box_to_delete=1:size(list_CTD_loaded_boxes_to_delete) + loaded_bbox.CTD=rmfield(loaded_bbox.CTD,(['box_' num2str(list_CTD_loaded_boxes_to_delete(i_box_to_delete))])); + end + end + +end + +if size(list_ARGO_loaded_boxes,1) > 1 + + list_ARGO_loaded_boxes=strrep(list_ARGO_loaded_boxes,'box_',''); + list_ARGO_loaded_boxes=str2double(list_ARGO_loaded_boxes(2:end)); + list_ARGO_loaded_boxes_to_delete=list_ARGO_loaded_boxes(~ismember(list_ARGO_loaded_boxes,pa_wmo_numbers(:,1))); + + if size(list_ARGO_loaded_boxes_to_delete)> 0 + for i_box_to_delete=1:size(list_ARGO_loaded_boxes_to_delete) + loaded_bbox.ARGO=rmfield(loaded_bbox.ARGO,(['box_' num2str(list_ARGO_loaded_boxes_to_delete(i_box_to_delete))])); + end + end + +end + +if size(list_BOTTLE_loaded_boxes,1) > 1 + + list_BOTTLE_loaded_boxes=strrep(list_BOTTLE_loaded_boxes,'box_',''); + list_BOTTLE_loaded_boxes=str2double(list_BOTTLE_loaded_boxes(2:end)); + list_BOTTLE_loaded_boxes_to_delete=list_BOTTLE_loaded_boxes(~ismember(list_BOTTLE_loaded_boxes,pa_wmo_numbers(:,1))); + + if size(list_BOTTLE_loaded_boxes_to_delete)> 0 + for i_box_to_delete=1:size(list_BOTTLE_loaded_boxes_to_delete) + loaded_bbox.BOTTLE=rmfield(loaded_bbox.BOTTLE,(['box_' num2str(list_BOTTLE_loaded_boxes_to_delete(i_box_to_delete))])); + end + end + +end -for ln_index = 1:length(pa_wmo_numbers) +% % DD (2024/08-3.1): correct for loop iteration when there is only one box to load. +if size(pa_wmo_numbers,1) == 4 + n_boxes=size(pa_wmo_numbers,2); +else + n_boxes=size(pa_wmo_numbers,1); +end +%for ln_index = 1:length(pa_wmo_numbers) +for ln_index = 1:n_boxes + for ntyp = 2:4 + if( ~isnan(pa_wmo_numbers(ln_index,1)) & pa_wmo_numbers(ln_index,ntyp) ) % check to see if we are supposed to load this data type + if ntyp == 2 % the 2nd column denotes CTD data - lo_box_data = load( strcat( po_config_data.HISTORICAL_DIRECTORY, po_config_data.HISTORICAL_CTD_PREFIX, sprintf( '%4d', pa_wmo_numbers(ln_index,1)))); + + % Test if it is already loaded % change DD (2024/08-3.2) + A=isfield(loaded_bbox.CTD,['box_' num2str(pa_wmo_numbers(ln_index,1))]); + if A == false + lo_box_data = load( strcat( po_config_data.HISTORICAL_DIRECTORY, po_config_data.HISTORICAL_CTD_PREFIX, sprintf( '%4d', pa_wmo_numbers(ln_index,1)))); + loaded_bbox.CTD.(['box_' num2str(pa_wmo_numbers(ln_index,1))])=lo_box_data; + else + lo_box_data=loaded_bbox.CTD.(['box_' num2str(pa_wmo_numbers(ln_index,1))]); + end + elseif ntyp == 3 % the 3rd column denotes historical data - lo_box_data = load( strcat( po_config_data.HISTORICAL_DIRECTORY, po_config_data.HISTORICAL_BOTTLE_PREFIX, sprintf( '%4d', pa_wmo_numbers(ln_index,1)))); - % In some cases the dimensions of the data may not match, causing issues with indexing below - if numel(lo_box_data.lat)==numel(lo_box_data.pres) && size(lo_box_data.lat,2)~=size(lo_box_data.pres,2) - lo_box_data.pres = lo_box_data.pres'; - lo_box_data.ptmp = lo_box_data.ptmp'; - lo_box_data.sal = lo_box_data.sal'; - lo_box_data.temp = lo_box_data.temp'; + + % Test if it is already loaded % change DD (2024/08-3.2) + A=isfield(loaded_bbox.BOTTLE,['box_' num2str(pa_wmo_numbers(ln_index,1))]); + if A == false + lo_box_data = load( strcat( po_config_data.HISTORICAL_DIRECTORY, po_config_data.HISTORICAL_BOTTLE_PREFIX, sprintf( '%4d', pa_wmo_numbers(ln_index,1)))); + % In some cases the dimensions of the data may not match, causing issues with indexing below + if numel(lo_box_data.lat)==numel(lo_box_data.pres) && size(lo_box_data.lat,2)~=size(lo_box_data.pres,2) + lo_box_data.pres = lo_box_data.pres'; + lo_box_data.ptmp = lo_box_data.ptmp'; + lo_box_data.sal = lo_box_data.sal'; + lo_box_data.temp = lo_box_data.temp'; + end + loaded_bbox.BOTTLE.(['box_' num2str(pa_wmo_numbers(ln_index,1))])=lo_box_data; + else + lo_box_data=loaded_bbox.BOTTLE.(['box_' num2str(pa_wmo_numbers(ln_index,1))]); + end elseif ntyp == 4 % the 4th column denotes Argo data - lo_box_data = load( strcat( po_config_data.HISTORICAL_DIRECTORY, po_config_data.HISTORICAL_ARGO_PREFIX, sprintf( '%4d', pa_wmo_numbers(ln_index,1)))); - % exclude Argo float being analysed from the Argo reference data selection, - % must do this step before concatenating the vectors, because "index" comes - % from "get_region_ow.m", which includes this step --------------------- - not_use=[]; - for i=1:length(lo_box_data.lat) - profile=lo_box_data.source{i}; - 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]; - end + + % Test if it is already loaded % change DD (2024/08-3.2) + A=isfield(loaded_bbox.ARGO,['box_' num2str(pa_wmo_numbers(ln_index,1))]); + if A == false + lo_box_data = load( strcat( po_config_data.HISTORICAL_DIRECTORY, po_config_data.HISTORICAL_ARGO_PREFIX, sprintf( '%4d', pa_wmo_numbers(ln_index,1)))); + % exclude Argo float being analysed from the Argo reference data selection, + % must do this step before concatenating the vectors, because "index" comes + % from "get_region_ow.m", which includes this step --------------------- + not_use=[]; + for i=1:length(lo_box_data.lat) + profile=lo_box_data.source{i}; + 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]; + end + end + lo_box_data.lat(not_use)=[]; + lo_box_data.long(not_use)=[]; + lo_box_data.dates(not_use)=[]; + lo_box_data.sal(:,not_use)=[]; + lo_box_data.ptmp(:,not_use)=[]; + lo_box_data.pres(:,not_use)=[]; + + loaded_bbox.ARGO.(['box_' num2str(pa_wmo_numbers(ln_index,1))])=lo_box_data; + else + lo_box_data=loaded_bbox.ARGO.(['box_' num2str(pa_wmo_numbers(ln_index,1))]); end - lo_box_data.lat(not_use)=[]; - lo_box_data.long(not_use)=[]; - lo_box_data.dates(not_use)=[]; - lo_box_data.sal(:,not_use)=[]; - lo_box_data.ptmp(:,not_use)=[]; - lo_box_data.pres(:,not_use)=[]; + %----------------------------------------------------------------------- end - % check index to see if this is a station that should be loaded + + % DD (2024/08-3.1) + % index is a concatenation of all stations inside the ellipse + % whatever the initial box. + % pa_index is a reconstruction of begin and end index of the + % current box. i1 = length(lo_box_data.lat); i2 = 1:i1; pa_index = pa_index_0 + i2; % indices for this set of data pa_index_0 = pa_index_0 + i1; % increment beginning index - - for n=1:i1 % load each station - ok = find(index == pa_index(n)); - if ~isempty(ok) % this entry is in the index list -%aw if n > length(lo_box_data.lat); keyboard; end - - jj = find( isnan(lo_box_data.pres(:,n))==0 ); % only retrieve non-NaN values - pres = lo_box_data.pres(jj,n); - sal = lo_box_data.sal(jj,n); - ptmp = lo_box_data.ptmp(jj,n); - - too_deep = find(pres>max_p); - pres(too_deep)=[]; - sal(too_deep)=[]; - ptmp(too_deep)=[]; - - new_depth = length(pres); - if( new_depth>max_depth & max_depth~=0 ) % patch up number of rows - pa_grid_pres = [ pa_grid_pres; ones( new_depth-max_depth, how_many_cols ) * NaN ]; - pa_grid_ptmp = [ pa_grid_ptmp; ones( new_depth-max_depth, how_many_cols ) * NaN ]; - pa_grid_sal = [ pa_grid_sal; ones( new_depth-max_depth, how_many_cols ) * NaN ]; - elseif( new_depth < max_depth ) - pres = [ pres; ones( max_depth-new_depth, 1 ) * NaN ]; - sal = [ sal ; ones( max_depth-new_depth, 1 ) * NaN ]; - ptmp = [ ptmp; ones( max_depth-new_depth, 1 ) * NaN ]; - end - - pa_grid_sal = [ pa_grid_sal, sal ] ; - pa_grid_ptmp = [ pa_grid_ptmp, ptmp ] ; - pa_grid_pres = [ pa_grid_pres, pres ] ; - - pa_grid_lat = [ pa_grid_lat, lo_box_data.lat(n) ] ; - pa_grid_long = [ pa_grid_long, lo_box_data.long(n) ] ; - pa_grid_dates = [ pa_grid_dates, lo_box_data.dates(n) ] ; - - [ max_depth, how_many_cols ] = size(pa_grid_pres); - + res=ismember(pa_index,index); + i_ok=find(res==1); + n_clim_profiles_to_load=length(i_ok); + %disp(['box_id=',num2str(pa_wmo_numbers(ln_index,1)),... + % ',ntyp=',num2str(ntyp),',n_clim_profiles_to_load=' num2str(n_clim_profiles_to_load)]) + + if n_clim_profiles_to_load > 0 % not really necessary anymore, only for safe-keeping + + % no need to re-scan for each + % station as the index is now saved. The extraction is also + % treated with array operations, instead of looping. + tmp_pa_grid_pres = lo_box_data.pres(:,i_ok) ; + tmp_pa_grid_sal = lo_box_data.sal(:,i_ok) ; + tmp_pa_grid_ptmp = lo_box_data.ptmp(:,i_ok) ; + + too_deep=tmp_pa_grid_pres>max_p; + tmp_pa_grid_pres(too_deep) = NaN; + tmp_pa_grid_sal(too_deep) = NaN; + tmp_pa_grid_ptmp(too_deep) = NaN; + + i_entire_NaN_rows=all(isnan(tmp_pa_grid_pres),2); + tmp_pa_grid_pres(i_entire_NaN_rows,:)=[]; + tmp_pa_grid_sal(i_entire_NaN_rows,:)=[]; + tmp_pa_grid_ptmp(i_entire_NaN_rows,:)=[]; + + [n_pres_level_1,n_obs_1]=size(pa_grid_pres); + [n_pres_level_2,n_obs_2]=size(tmp_pa_grid_pres); + d_pres_level=n_pres_level_2-n_pres_level_1; + if n_pres_level_1 > 0 + if d_pres_level > 0 + nan_array=NaN*ones(d_pres_level,n_obs_1); + pa_grid_pres=cat(2,cat(1,pa_grid_pres,nan_array), tmp_pa_grid_pres); + pa_grid_sal =cat(2,cat(1,pa_grid_sal ,nan_array), tmp_pa_grid_sal ); + pa_grid_ptmp=cat(2,cat(1,pa_grid_ptmp,nan_array), tmp_pa_grid_ptmp); + + else + if d_pres_level < 0 + nan_array=NaN*ones(-d_pres_level,n_obs_2); + pa_grid_pres=cat(2,pa_grid_pres,cat(1,tmp_pa_grid_pres,nan_array)); + pa_grid_sal =cat(2,pa_grid_sal ,cat(1,tmp_pa_grid_sal ,nan_array)); + pa_grid_ptmp=cat(2,pa_grid_ptmp,cat(1,tmp_pa_grid_ptmp,nan_array)); + else + pa_grid_pres=[pa_grid_pres, tmp_pa_grid_pres]; + pa_grid_sal =[pa_grid_sal , tmp_pa_grid_sal ]; + pa_grid_ptmp=[pa_grid_ptmp, tmp_pa_grid_ptmp]; + end + + end + else + pa_grid_pres=[pa_grid_pres, tmp_pa_grid_pres]; + pa_grid_sal =[pa_grid_sal , tmp_pa_grid_sal ]; + pa_grid_ptmp=[pa_grid_ptmp, tmp_pa_grid_ptmp]; end - end %n=1:i1 - end + + pa_grid_lat = [pa_grid_lat , lo_box_data.lat(i_ok) ]; + pa_grid_long = [pa_grid_long , lo_box_data.long(i_ok) ]; + pa_grid_dates = [pa_grid_dates , lo_box_data.dates(i_ok) ]; + + end % if n_clim_profiles_to_load > 0 + % end of DD (2024/08-3.1) + +% %for n=1:i1 % load each station +% +% %ok = find(index == pa_index(n)); +% %if ~isempty(ok) % this entry is in the index list +% %aw if n > length(lo_box_data.lat); keyboard; end +% +% jj = find( isnan(lo_box_data.pres(:,n))==0 ); % only retrieve non-NaN values +% pres = lo_box_data.pres(jj,n); +% sal = lo_box_data.sal(jj,n); +% ptmp = lo_box_data.ptmp(jj,n); +% +% too_deep = find(pres>max_p); +% pres(too_deep)=[]; +% sal(too_deep)=[]; +% ptmp(too_deep)=[]; +% +% new_depth = length(pres); +% if( new_depth>max_depth & max_depth~=0 ) % patch up number of rows +% pa_grid_pres = [ pa_grid_pres; ones( new_depth-max_depth, how_many_cols ) * NaN ]; +% pa_grid_ptmp = [ pa_grid_ptmp; ones( new_depth-max_depth, how_many_cols ) * NaN ]; +% pa_grid_sal = [ pa_grid_sal; ones( new_depth-max_depth, how_many_cols ) * NaN ]; +% elseif( new_depth < max_depth ) +% pres = [ pres; ones( max_depth-new_depth, 1 ) * NaN ]; +% sal = [ sal ; ones( max_depth-new_depth, 1 ) * NaN ]; +% ptmp = [ ptmp; ones( max_depth-new_depth, 1 ) * NaN ]; +% end +% +% pa_grid_sal = [ pa_grid_sal, sal ] ; +% pa_grid_ptmp = [ pa_grid_ptmp, ptmp ] ; +% pa_grid_pres = [ pa_grid_pres, pres ] ; +% +% pa_grid_lat = [ pa_grid_lat, lo_box_data.lat(n) ] ; +% pa_grid_long = [ pa_grid_long, lo_box_data.long(n) ] ; +% pa_grid_dates = [ pa_grid_dates, lo_box_data.dates(n) ] ; +% +% [ max_depth, how_many_cols ] = size(pa_grid_pres); + +% end +% end % for n=1:i1 + + end % if( ~isnan(pa_wmo_numbers(ln_index,1)) & pa_wmo_numbers(ln_index,ntyp) ) end %ntyp -end %ln_index +end % for ln_index = 1:size(pa_wmo_numbers,1) % longitude goes from 0 to 360 degrees diff --git a/matlab_codes/set_calseries.m b/matlab_codes/set_calseries.m index d597de7..78340f7 100755 --- a/matlab_codes/set_calseries.m +++ b/matlab_codes/set_calseries.m @@ -6,8 +6,16 @@ function set_calseries( pn_float_dir, pn_float_name, po_system_configuration ) % % Annie Wong, September 2008 % Breck Owens, October 2006 -% +% +% Delphine Dobler (DD), August 2024: +% 2 - create output directory if they do not exist already + +% DD (2024/08-2) +outdir=[ po_system_configuration.FLOAT_CALIB_DIRECTORY pn_float_dir ]; +if not(isfolder(outdir)) + mkdir(outdir) +end % load data --- diff --git a/matlab_codes/update_salinity_mapping.m b/matlab_codes/update_salinity_mapping.m index 0f837ba..b0e6cd3 100755 --- a/matlab_codes/update_salinity_mapping.m +++ b/matlab_codes/update_salinity_mapping.m @@ -20,6 +20,23 @@ function update_salinity_mapping( pn_float_dir, pn_float_name, po_system_configu % Annie Wong, June 2020: rename "scale_age" to "scale_age_small", and create "scale_age_large" with NaN for old OW mapped profiles. % This is to make the temporal scale variable names consistent with the lat/long scale variable names in OWC. % +% Delphine Dobler (DD), August 2024: +% 1 - Auto-create la_wmo_boxes from climatological repositores when la_wmo_boxes does not exist +% 2 - create output directory if they do not exist already +% 3.1 - Performance: save boxes_id from get_region_ow and +% recompute the index to provide as input to retr_region_ow +% 3.2 - Performance: Inside retr_region_ow, save loaded data, +% and at each following profile, reload only new boxes and delete +% unused. +% 3.3 - Performance: Inside map_data_grid, comment calculation +% of vdataerror as this is never used afterwards. +% + +% DD (2024/08-2) +outdir=[ po_system_configuration.FLOAT_MAPPED_DIRECTORY pn_float_dir ]; +if not(isfolder(outdir)) + mkdir(outdir) +end % load float source data ---------------------------------------- @@ -31,8 +48,12 @@ function update_salinity_mapping( pn_float_dir, pn_float_name, po_system_configu % load mapping parameters ------- - -load( strcat( po_system_configuration.CONFIG_DIRECTORY, po_system_configuration.CONFIG_WMO_BOXES ), 'la_wmo_boxes' ) ; +% DD : 2024/08-1 +la_wmo_boxes_file=strcat( po_system_configuration.CONFIG_DIRECTORY, po_system_configuration.CONFIG_WMO_BOXES); +if ~exist(la_wmo_boxes_file) + create_la_wmo_boxes_file(po_system_configuration) +end +load( la_wmo_boxes_file, 'la_wmo_boxes' ) ; ln_max_casts = str2num( po_system_configuration.CONFIG_MAX_CASTS ) ; map_use_pv = str2num( po_system_configuration.MAP_USE_PV ); @@ -103,10 +124,18 @@ function update_salinity_mapping( pn_float_dir, pn_float_name, po_system_configu end clear a i +% DD (2024/08-3.2): +n_profiles=length( missing_profile_index ); +% loaded_box record initialisation +loaded_bbox.n_profiles=n_profiles; +loaded_bbox.CTD.description="CTD loaded_boxes"; +loaded_bbox.ARGO.description="ARGO loaded_boxes"; +loaded_bbox.BOTTLE.description="BOTTLE loaded_boxes"; + % update mapped data matrix by missing_profile_index -------------------- -for i = 1 : length( missing_profile_index ) +for i = 1 : n_profiles tic disp(['UPDATE_SALINITY_MAPPING: Working on profile ' num2str(i)]) @@ -161,7 +190,10 @@ function update_salinity_mapping( pn_float_dir, pn_float_name, po_system_configu [ la_wmo_numbers ] = find_25boxes( LONG, LAT, la_wmo_boxes ) ; % load in the positions and dates for the surrounding wmo boxes - [ la_grid_lat, la_grid_long, la_grid_dates ] = get_region_ow( la_wmo_numbers, pn_float_name, po_system_configuration ) ; + %[ la_grid_lat, la_grid_long, la_grid_dates ] = get_region_ow( la_wmo_numbers, pn_float_name, po_system_configuration ) ; + % DD (2024/08-3.1) add an output to the function: the box id + [ la_grid_lat, la_grid_long, la_grid_dates, pa_grid_box_id ] = get_region_ow( la_wmo_numbers, pn_float_name, po_system_configuration ) ; + if( la_grid_lat~=999 ) % if no historical data is assigned to the float profile @@ -171,14 +203,14 @@ function update_salinity_mapping( pn_float_dir, pn_float_name, po_system_configu gg=find(la_grid_long>180); la_grid_long1(gg)=la_grid_long(gg)-360; % m_tbase inputs longitudes from 0 to +/- 180 - m_proj('mercator','long', [min(la_grid_long1)-1, max(la_grid_long1)+1], 'lat', [min(la_grid_lat)-1, max(la_grid_lat)+1] ); + m_proj('mercator','long', [min(la_grid_long1)-1, max(la_grid_long1)+1], 'lat', [min(la_grid_lat)-1, max(la_grid_lat)+1] ); [elev,x,y] = m_tbase( [min(la_grid_long1)-1, max(la_grid_long1)+1, min(la_grid_lat)-1, max(la_grid_lat)+1] ); la_grid_Z = -interp2( x,y,elev, la_grid_long1, la_grid_lat, 'linear'); % -ve bathy values % make LONG compatiable with la_grid_long at the 0-360 mark LONG2=LONG; - if(isempty(find(la_grid_long>360))==0) + if(isempty(find(la_grid_long>360))==0) if(LONG>=0&LONG<=20) LONG2=LONG+360; end @@ -187,11 +219,54 @@ function update_salinity_mapping( pn_float_dir, pn_float_name, po_system_configu % find ln_max_casts historical points that are most strongly correlated with the float profile %[ index ] = find_besthist( la_grid_lat, la_grid_long, la_grid_dates, la_grid_Z, LAT, LONG2, DATES, Z, latitude_large, latitude_small, longitude_large, longitude_small, phi_large, phi_small, map_age, map_use_pv, ln_max_casts ); % change config 129 - [ index ] = find_besthist( la_grid_lat, la_grid_long, la_grid_dates, la_grid_Z, LAT, LONG2, DATES, Z, latitude_large, latitude_small, longitude_large, longitude_small, phi_large, phi_small, map_age_small, map_age_large, map_use_pv, ln_max_casts ); %AW June 2020 + [ index_ini ] = find_besthist( la_grid_lat, la_grid_long, la_grid_dates, la_grid_Z, LAT, LONG2, DATES, Z, latitude_large, latitude_small, longitude_large, longitude_small, phi_large, phi_small, map_age_small, map_age_large, map_use_pv, ln_max_casts ); %AW June 2020 clear la_grid_lat la_grid_long la_grid_dates - [ la_bhist_sal, la_bhist_ptmp, la_bhist_pres, la_bhist_lat, la_bhist_long, la_bhist_dates ] = retr_region_ow( la_wmo_numbers, pn_float_name, po_system_configuration, index, PRES, map_p_delta ) ; + % DD (2024/08-3.1) Loading boxes' data is a consuming step in terms of performance. + % The test if the boxes' data should be kept was initially done after the loading + % inside retr_region_ow. Here we recompute the index for the boxes ids to load + % (now known from get_region_ow new output) beforehand to save some execution time: + % 'stable' option is very important in the unique function call. + box_id_2_keep=unique(pa_grid_box_id(index_ini),'stable'); + n_hist_obs=size(index_ini,1); + if size(box_id_2_keep,1)==1 + n_boxes=size(box_id_2_keep,2); + else + n_boxes=size(box_id_2_keep,1); + end + la_wmo_numbers=NaN(n_boxes,4); + n_prev_box=0; + new_global_index = []; + disp(['Fitted observations are located in ' num2str(n_boxes) ' boxes with ' num2str(n_hist_obs) ' historical observations found']) + + for ibox=1:n_boxes + + % extract the new la_wmo_numbers + la_wmo_numbers(ibox,:)=la_wmo_boxes(la_wmo_boxes(:,1)==box_id_2_keep(ibox),:); + % disp(['loaded box ' num2str(la_wmo_numbers(ibox,:))]) + % recompute the index array + full_index_list_for_ibox=find(pa_grid_box_id==box_id_2_keep(ibox)); + index_min=full_index_list_for_ibox(1); + index_max=full_index_list_for_ibox(end); + fitted_index_list_for_ibox=index_ini( (index_ini >=index_min) & (index_ini <=index_max) ); + local_index=fitted_index_list_for_ibox-index_min+1; + new_global_index=[new_global_index; local_index+n_prev_box]; + if(size(full_index_list_for_ibox,1)==1) + n_obs_box=size(full_index_list_for_ibox,2); + else + n_obs_box=size(full_index_list_for_ibox,1); + end + n_prev_box=n_obs_box+n_prev_box; + end + + index=new_global_index; + % end of DD (2024/08-3.1) + + + + [ la_bhist_sal, la_bhist_ptmp, la_bhist_pres, la_bhist_lat, la_bhist_long, la_bhist_dates, loaded_bbox ] = ... + retr_region_ow( la_wmo_numbers, pn_float_name, po_system_configuration, index, PRES, map_p_delta, loaded_bbox ) ; la_bhist_Z = la_grid_Z(index); % include JB's SAF frontal separation criteria if map_use_saf==1 @@ -201,7 +276,7 @@ function update_salinity_mapping( pn_float_dir, pn_float_name, po_system_configu [ la_bhist_sal2, la_bhist_ptmp2, la_bhist_pres2, la_bhist_lat2, la_bhist_long2, la_bhist_dates2, la_bhist_Z2 ] = ... frontalConstraintSAF(la_bhist_sal, la_bhist_ptmp, la_bhist_pres, la_bhist_lat, la_bhist_long, la_bhist_dates, la_bhist_Z, LAT, LONG, PRES, PTMP, SAL, po_system_configuration); - [a2 b2]= size(la_bhist_sal2); + [a2, b2]= size(la_bhist_sal2); if (b2>5) % Use frontal separation only if there are at least 5 profiles la_bhist_sal=la_bhist_sal2; la_bhist_ptmp=la_bhist_ptmp2; @@ -297,7 +372,7 @@ function update_salinity_mapping( pn_float_dir, pn_float_name, po_system_configu % longitude_large, latitude_large, NaN, signal_sal, noise_sal, phi_large, map_use_pv ) ; - [a1,b1,c1,d1]... + [a1,b1,c1,~]... = map_data_grid( la_hist_sal, [ LAT, LONG, DATES, Z ], ... [ la_hist_lat, la_hist_long, la_hist_dates, la_hist_Z ], ... longitude_large, latitude_large, map_age_large, signal_sal, noise_sal, phi_large, map_use_pv ) ; @@ -307,7 +382,7 @@ function update_salinity_mapping( pn_float_dir, pn_float_name, po_system_configu la_residualsal1 = la_hist_sal - c1' ; la_signalresidualsal = signal( la_residualsal1 ) ; - [a2,b2,c2,d2]... + [a2,b2,c2,~]... = map_data_grid( la_residualsal1, [ LAT, LONG, DATES, Z ], ... [ la_hist_lat, la_hist_long, la_hist_dates, la_hist_Z ], ... longitude_small, latitude_small, map_age_small, la_signalresidualsal, noise_sal, phi_small, map_use_pv ) ; %AW June 2020 @@ -337,11 +412,11 @@ function update_salinity_mapping( pn_float_dir, pn_float_name, po_system_configu selected_hist=[la_hist_long(1),la_hist_lat(1),la_profile_no(ln_profile_index)]; end for k = 1 : length(la_hist_long) - [m,n] = size(selected_hist); + [m,~] = size(selected_hist); b = [ la_hist_long(k), la_hist_lat(k) ]; c = selected_hist(:,1:2) - ones(m,1)*b; d = find( abs(c(:,1))<1/60&abs(c(:,2))<1/60 ); %within 1 min, do not save - if( isempty(d)==1 ) + if( isempty(d)==1 ) selected_hist = [ selected_hist; [ la_hist_long(k), la_hist_lat(k), la_profile_no(ln_profile_index) ] ]; end end diff --git a/ow_calibration.m b/ow_calibration.m index 15d13a0..d472903 100755 --- a/ow_calibration.m +++ b/ow_calibration.m @@ -10,6 +10,9 @@ % float_names={'R49000139';'R39033'}; % % these variables have to be set before ow_calibration is called. +% +% Delphine Dobler (DD), August 2024: Add date of end of proccessing for +% performance assessment. % lo_system_configuration = load_configuration( 'ow_config.txt' ); @@ -29,6 +32,8 @@ calculate_piecewisefit( flt_dir, flt_name, lo_system_configuration ); plot_diagnostics_ow( flt_dir, flt_name, lo_system_configuration ); + + disp([datestr(now) ' End of processing for ' flt_name]) end