Skip to content

Commit

Permalink
Merge pull request #1 from ArgoDMQC/develop
Browse files Browse the repository at this point in the history
Develop
  • Loading branch information
cabanesc committed May 24, 2018
2 parents 55c7b33 + 4d58f58 commit 74d45ba
Show file tree
Hide file tree
Showing 14 changed files with 399 additions and 51 deletions.
2 changes: 2 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -1,2 +1,4 @@
ow_config_*
ow_calibration_*
set_calseries_*
~$*
65 changes: 63 additions & 2 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,2 +1,63 @@
# matlabow
The OW Matlab toolbox
# The Matlab OW toolbox

This is a package of MATLAB routines for calibrating profiling float conductivity sensor drift. A description of the algorithms can be found in "An improved calibration method for the drift of the conductivity sensor on autonomous CTD profiling floats by θ-S climatology", by W.B. Owens and A.P.S. Wong, in Deep-Sea Research Part I: Oceanographic Research Papers, 56(3), 450-457, 2009.
Lately, modifications suggested in “Improvement of bias detection in Argo float conductivity sensors and its application in the North Atlantic” , by C. Cabanes, V. Thierry and C. Lagadec, in Deep-Sea Research Part I, 114, 128-136, 2016 have been taken into account.

# How to install the toolbox ?

Either clone the latest version of the git repository:

git clone https://github.com/ArgoDMQC/matlabow.git

or download and unzip the zip file (Clone or download button)

You can also access the different releases here:
https://github.com/ArgoDMQC/matlabow/releases

# How to run the analysis?
Here is a summary of what should be done to run the analysis, please read the ./doc/README.doc file for more details


1. All files are to be used in MATLAB. The full package was tested with MATLAB R2014a. In addition, you will need:
a). The MATLAB Optimization Toolbox;
b). The ITS-90 version of the CSIRO SEAWATER library. The version 3\_3.0 of this library can be found in ./lib/seawater\_330\_its90. Please update if necessary.
c). The M_MAP toolbox. The version 1.4.c of this library can be found in ./lib/m\_map1.4. Please update if necessary.

2. Add the necessary path to your matlab path: addpath('./lib/seawater\_330\_its90';'./lib/m\_map1.4';'./matlab\_codes/')

3. Put your reference data in ./data/climatology/historical\_ctd, /historical\_bot, /argo\_profiles.

REFERENCE DATA can be obtain at ftp.ifremer.fr
cd /coriolis/data/DMQC-ARGO/ (if you need a login/pswd ask codac@ifremer.fr)

Then, create/update your ./data/constants/wmo\_boxes.mat file (more details in ./doc/README.doc, p3)

4. After you have decided where you want to install the package on your computer, edit ow\_config.txt at the following lines so the correct pathways are specified:

* HISTORICAL\_DIRECTORY =
* FLOAT\_SOURCE\_DIRECTORY =
* FLOAT\_MAPPED\_DIRECTORY =
* FLOAT\_CALIB\_DIRECTORY =
* FLOAT\_PLOTS\_DIRECTORY =
* CONFIG\_DIRECTORY =

5. The last section of ow\_config.txt below the heading "Objective Mapping Parameters" is where you set the various parameters (more details in .doc.README.doc, p4-6)

6. If this is the first time you are using this system, then the 4 directories /data/float\_source, /float\_mapped, /float\_calib, and /float\_plots should be empty. Decide how you want to organise your floats, e.g. under different project names or different investigator names. Then make identical subdirectories under each of these 4 directories. For example:

/data/float\_source/project\_xx
/data/float\_mapped/project\_xx
/data/float\_calib/project\_xx
/data/float\_plots/project\_xx

7. Create the float source file (./data/float\_source/project\_xx/$flt\_name$.mat) from the original netcdf files (more details in ./doc/README.doc,p6)

8. Open MATLAB in the top directory. List all the float files in a cell array "float\_names", with the corresponding subdirectories in another cell array "float_dirs". For example,
float_dirs = { 'project\_xx/'; 'project\_xx/'; 'jones/'; 'jones/' };
float_names = { 'float0001'; 'float0002'; 'myfloat\_a'; 'myfloat\_b' }.

Tips: If the files are not saved under a subdirectory and are only saved under ./float\_source/, specify float\_dirs = { ''; ''; ''; '' }, etc.

Run ow\_calibration.m.


Binary file added doc/README.doc
Binary file not shown.
Binary file removed doc/README.pdf
Binary file not shown.
142 changes: 142 additions & 0 deletions matlab_codes/build_cov.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,142 @@

function cov = build_cov( ptmp,coord_float,po_system_configuration);

% This function was created from build_ptmp_cov which the description is written below:
%
% Description of build_ptmp_cov
% ____
% This function builds the covariance matrix (a square matrix) that has nxn number
% of tiles, each tile is of mxm size.
%
% The vertical covariance matrix is the building tile. It has 1 at its diagonal,
% then decreases exponentially in the off-diagonals, which represents the vertical
% covariance between water masses.
%
% This version assums that each profile is independent from the other ones.
%
% Annie Wong, May 2001
% Breck Owens, November 2007
%____
%
% Cecile Cabanes, June 2013
% Compared to build_ptmp_cov, build_cov( ptmp,coord_float,po_system_configuration) take into account the horizontal covariance between different mapped profiles.
% This lateral covariance take into account the fact that a mapped profile on a Argo position is build from a set of historical profiles that is not very different from the set used to build a mapped profile at the next or previous Argo profile position.
% This lateral covariance between two mapped profiles is constructed using a Gaussian function and the large spatial scales
%
% Cecile Cabanes, 2017
% use the small spatial scales instead of the large spatial scales to build
% the lateral covariance: found to be the best compromise between
% informative errors and large enough NDF for AIC criterium, at least for
% the scales defined for the North Atlantic bassin

% Set up the theta boundaries for water masses

ptboundaries = [ 30; 24; 18; 12; 8; 4; 2.5; 1; -2 ];
ptscale_down = [ 6; 6; 6; 4; 4; 1.5; 1.5; 1; 1 ];
ptscale_up = [ 6; 6; 6; 6; 4; 4; 1.5; 1.5; 1 ];


% set up the building tile = vertical covariance matrix
%
% upper triangle of the matrix = covariance of each ptlevel with every ptlevel below it,
% looking down the water column from the diagonal
% lower triangle of the matrix = covariance of each ptlevel with every ptlevel above it,
% looking up the water column from the diagonal

[m,n]=size(ptmp);

cov=zeros(m*n,m); % set up the output covariance matrix, with 1's along the diagonal


for k=1:n % for each profile
k1 = (k-1)*m;
for l=1:1 % for each profile
l1 = (l-1)*m;
for i=1:m % for all levels
for j=1:m
if(i<j) % upper triangle, look down the water column for vertical scale
Ltheta = interp1( ptboundaries, ptscale_down, ptmp(i,k), 'linear' );
cov(i+k1,j+l1) = exp( - ( ptmp(j,k) - ptmp(i,k) ).^2/ Ltheta.^2 );
elseif(i>j) % lower triangle, look up the water column for vertical scale
Ltheta = interp1( ptboundaries, ptscale_up, ptmp(i,k), 'linear' );
cov(i+k1,j+l1) = exp( - ( ptmp(j,k) - ptmp(i,k) ).^2/ Ltheta.^2 );
elseif(i==j)
cov(i+k1,j+l1) =1;
end
if isnan(cov(i+k1,j+l1))
cov(i+k1,j+l1)=1;
end
end
end
end
end


%
% set up the horizontal covariance matrix
scale_lon = str2num(po_system_configuration.MAPSCALE_LONGITUDE_SMALL);
scale_lat = str2num(po_system_configuration.MAPSCALE_LATITUDE_SMALL);
scale_phi = str2num(po_system_configuration.MAPSCALE_PHI_SMALL);
map_use_pv = str2num(po_system_configuration.MAP_USE_PV);

for k=1:n
b(k,:)=covarxy_pv(coord_float(k,:),coord_float,scale_lon,scale_lat,scale_phi,map_use_pv);
end
b=b(:,1:n);

covh=b;

%keyboard
% build the final covariance matrix, taking into account horizontal and vertical covariance

covn=repmat(cov,[1,n]);
for k=1:n % for each profile
kd = (k-1)*m+1;
kf = k*m;
for l=1:n % for each profile
ld = (l-1)*m+1;
lf = (l)*m;
covn(kd:kf,ld:lf)=covh(k,l)*covn(kd:kf,ld:lf);
end
end

cov=covn;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% function [a] = covarxy_pv(x1,x2,gs,ts,phi,pv)

function [a] = covarxy_pv(x1,x2,gs,ts, phi, pv)

m=size(x1);
n=size(x2);
a=zeros(m(1),n(1));
one=ones(n(1),1);


% derive planetary vorticity at each data point ----

Zx1=x1(:,3); % depth at x1 datapoint
Zx2=x2(:,3); % depth at x2 datapoint

PV_x1=(2*7.292*10^-5.*sin(x1(:,1).*pi/180))./Zx1;
PV_x2=(2*7.292*10^-5.*sin(x2(:,1).*pi/180))./Zx2;


% calculate covariance term, use pv as optional ----

for i=1:m(1)
tmp=((x1(i,1)-x2(:,1))./ts).^2 + ...
((x1(i,2)-x2(:,2))./gs).^2;
if(pv &PV_x1~=0&PV_x2~=0) % if PV is wanted and the point is not on the equator ---
tmp = tmp + ...
( (PV_x1(i)-PV_x2)./sqrt( PV_x1(i).^2+PV_x2.^2 )./phi ).^2;
end


a(i,:)=exp(-tmp');
end

return
43 changes: 36 additions & 7 deletions matlab_codes/calculate_piecewisefit.m
Original file line number Diff line number Diff line change
@@ -1,10 +1,13 @@

function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configuration )

%
% Annie Wong, October 2008
% Breck Owens, November 2007


%
% Cecile Cabanes, June 2013 : calculate off-diagonal terms for error estimate: add horizontal covariance to track changes: see "change config 129"

%pn_float_dir='testfloats/';
%pn_float_name='robbins4900178';
Expand All @@ -31,6 +34,29 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur
la_ptmp = lo_float_mapped_data.la_ptmp; % float potential temperature where mapping is done



%+++++ change config 129----FROM HERE

% retrieve coordinate (XYZ) of the float position (coord_float) that is used in build_cov.m
if ~isempty(lo_float_mapped_data.selected_hist)

% Calculate elevation at the float position

if(LONG>180) % m_tbase inputs longitudes from 0 to +/- 180
LONG1=LONG-360;
else
LONG1=LONG;
end
m_proj('mercator','long', [min(LONG1)-1, max(LONG1)+1], 'lat', [min(LAT)-1, max(LAT)+1] );
[elev,x,y] = m_tbase( [min(LONG1)-1, max(LONG1)+1, min(LAT)-1, max(LAT)+1] );
Z = -interp2( x,y,elev, LONG1, LAT, 'linear'); % -ve bathy values

coord_float=[LONG',LAT',Z'];

end
%+++++ change config 129-----TO HERE


% load calibration settings -----------------

lo_float_calseries = load( strcat( po_system_configuration.FLOAT_CALIB_DIRECTORY, pn_float_dir, po_system_configuration.FLOAT_CALSERIES_PREFIX , pn_float_name, po_system_configuration.FLOAT_MAPPED_POSTFIX ) );
Expand Down Expand Up @@ -170,7 +196,10 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur

% calculate off-diagonal terms for error estimate --------

covariance = build_ptmp_cov(ten_PTMP); % build the data covariance matrix
%covariance = build_ptmp_cov(ten_PTMP); % build the data covariance matrix % vertical covariance only

covariance = build_cov(ten_PTMP,coord_float,po_system_configuration); % change config 129 : vertical and horizontal covariances


% for debugging purposes to speed up calculations, use next line for first time calculation
% and then comment out the call to build_ptmp_cov and load the covariance matrix
Expand All @@ -186,7 +215,7 @@ 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(i));
else
breaks_in = breaks(i,:);
breaks_in = breaks_in(find(isfinite(breaks_in)));
Expand All @@ -197,7 +226,7 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur
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(i));
end
end

Expand All @@ -208,13 +237,13 @@ function calculate_piecewisefit( pn_float_dir, pn_float_name, po_system_configur
unique_COND = sw_c3515*sw_cndr( unique_SAL, unique_PTMP, 0);
cal_COND(:,calindex) = ( ones(m,1)*pcond_factor(calindex) ).*unique_COND;
cal_SAL(:,calindex) = sw_salt( cal_COND(:,calindex)/sw_c3515, unique_PTMP, 0);

%pcond_factor(calindex)
% estimate the error in salinity ---------------------------------

%pcond_factor_err(calindex)
cal_COND_err(:,calindex) = ( ones(m,1)*pcond_factor_err(calindex) ).*unique_COND;
cal_SAL1(:,calindex) = sw_salt( (cal_COND(:,calindex)+cal_COND_err(:,calindex))/sw_c3515, unique_PTMP, 0 );
cal_SAL_err(:,calindex) = abs(cal_SAL(:,calindex)-cal_SAL1(:,calindex));

% keyboard
% estimate the error in salinity for station by station fit ----

sta_COND(:,calindex) = ( ones(m,1)*sta_mean(calindex) ).*unique_COND;
Expand Down
57 changes: 43 additions & 14 deletions matlab_codes/changedates.m
Original file line number Diff line number Diff line change
Expand Up @@ -6,26 +6,55 @@
%
% A. Wong, 29 May 2001
%
% C.Cabanes June 2013 : speed up this function

function [dates]=changedates(dates_format2);


% dates_format2 to single column

if size(dates_format2,2)>1
dates_format2=dates_format2';
end

dates=NaN.*ones(length(dates_format2),1); %organise dates in a single column

for i=1:length(dates_format2)
if(isnan(dates_format2(i))==0)
junk=int2str(dates_format2(i));
yr=str2num(junk(:,1:4));
mo=str2num(junk(:,5:6));
day=str2num(junk(:,7:8));
hr=str2num(junk(:,9:10));
min=str2num(junk(:,11:12));
if(mo<1|mo>12|day<1|day>31)
dates(i)=yr;
else
dates(i)=yr+cal2dec(mo,day,hr,min)./365;
end
end
indnotnan=find(~isnan(dates_format2'));

junk=int2str(dates_format2(indnotnan));

if length(dates_format2)>0

yr=str2num(junk(:,1:4));
mo=str2num(junk(:,5:6));
day=str2num(junk(:,7:8));
hr=str2num(junk(:,9:10));
min=str2num(junk(:,11:12));

indselect=indnotnan(mo<1|mo>12|day<1|day>31);

dates(indselect)=yr(indselect);

isok=(mo>=1&mo<=12&day>=1&day<=31);
indselect=indnotnan(isok);
dates(indselect)=yr(isok)+cal2dec(mo(isok),day(isok),hr(isok),min(isok))./365;

end
%
% for i=1:length(dates_format2)
%
% if(isnan(dates_format2(i))==0)
% junk=int2str(dates_format2(i));
% yr=str2num(junk(:,1:4));
% mo=str2num(junk(:,5:6));
% day=str2num(junk(:,7:8));
% hr=str2num(junk(:,9:10));
% min=str2num(junk(:,11:12));
% if(mo<1|mo>12|day<1|day>31)
% dates(i)=yr;
% else
% dates(i)=yr+cal2dec(mo,day,hr,min)./365;
% end
% end
% end

Loading

0 comments on commit 74d45ba

Please sign in to comment.