diff --git a/AutomaticQC/imosDensityInversionSetQC.m b/AutomaticQC/imosDensityInversionSetQC.m new file mode 100644 index 000000000..a4998257c --- /dev/null +++ b/AutomaticQC/imosDensityInversionSetQC.m @@ -0,0 +1,179 @@ +function [sample_data, varChecked, paramsLog] = imosDensityInversionSetQC( sample_data, auto ) +%DENSITYINVERSIONSETQC Flags any PSAL + TEMP + CNDC + DENS value that shows an inversion in density +% or an increase/decrease in density that is greater than a threshold set in +% imosDensityInversionQC.txt. +% +% Density inversion test on profiles data from ARGO which finds and flags +% good any data which value Vn passes the tests: +% - density difference when going down <= threshold +% - density difference when going up >= -threshold +% +% Inputs: +% sample_data - struct containing the data set. +% +% auto - logical, run QC in batch mode +% +% Outputs: +% sample_data - same as input, with QC flags added for variable/dimension +% data. +% +% varChecked - cell array of variables' name which have been checked +% +% paramsLog - string containing details about params' procedure to include in QC log +% +% Author: Guillaume Galibert +% + +% +% Copyright (C) 2017, Australian Ocean Data Network (AODN) and Integrated +% Marine Observing System (IMOS). +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation version 3 of the License. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. + +% You should have received a copy of the GNU General Public License +% along with this program. +% If not, see . +% + +narginchk(1, 2); +if ~isstruct(sample_data), error('sample_data must be a struct'); end + +% auto logical in input to enable running under batch processing +if nargin<5, auto=false; end + +varChecked = {}; +paramsLog = []; + +% this test only applies to profile mode +mode = readProperty('toolbox.mode'); +if ~strcmpi(mode, 'profile') + return; +end + +% this test only applies on PSAL when TEMP and pressure data is available +tempIdx = getVar(sample_data.variables, 'TEMP'); +psalIdx = getVar(sample_data.variables, 'PSAL'); + +[presRel, zName, ~] = getPresRelForGSW(sample_data); +[lat, lon] = getLatLonForGSW(sample_data); + +if ~(psalIdx && tempIdx && ~isempty(zName) && ~isempty(lat) && ~isempty(lon)), return; end + +temp = sample_data.variables{tempIdx}.data; +psal = sample_data.variables{psalIdx}.data; + +% read threshold value from properties file +threshold = str2double(readProperty('threshold', fullfile('AutomaticQC', 'imosDensityInversionSetQC.txt'))); + +% read dataset QC parameters if exist and override previous parameters file +currentQCtest = mfilename; +threshold = readDatasetParameter(sample_data.toolbox_input_file, currentQCtest, 'threshold', 0.03); % default threshold value of 0.03kg/m3 + +qcSet = str2double(readProperty('toolbox.qc_set')); +rawFlag = imosQCFlag('raw', qcSet, 'flag'); +goodFlag = imosQCFlag('good', qcSet, 'flag'); +probGoodFlag = imosQCFlag('probablyGood',qcSet, 'flag'); +probBadFlag = imosQCFlag('probablyBad', qcSet, 'flag'); +badFlag = imosQCFlag('bad', qcSet, 'flag'); + +badFlags = [probBadFlag, badFlag]; + +paramsLog = ['threshold=' num2str(threshold)]; + +% matrix case, we unfold the matrix in one vector for profile study +% purpose +isMatrix = size(temp, 1)>1 & size(temp, 2)>1; +if isMatrix + len1 = size(temp, 1); + len2 = size(temp, 2); + len3 = size(temp, 3); + temp = temp(:); +end + +lenData = length(temp); + +% we don't consider already bad data in the current test +tempFlags = sample_data.variables{tempIdx}.flags; +psalFlags = sample_data.variables{psalIdx}.flags; + +zIdx = getVar(sample_data.variables, zName); +if zIdx == 0 + % we either have DEPTH as a dimension or we have NOMINAL_DEPTH: in both + % cases we don't have any flag information, then it is safe to duplicate TEMP flags + % for example. + presFlags = tempFlags; +else + presFlags = sample_data.variables{zIdx}.flags; +end + +iBadData = ismember(tempFlags, badFlags) | ismember(psalFlags, badFlags) | ismember(presFlags, badFlags); + +% we need at least 2 valid measurements to perform this test +if sum(~iBadData) <= 1, return; end + +temp = temp(~iBadData); +psal = psal(~iBadData); +presRel = presRel(~iBadData); + +lenDataTested = length(temp); + +flags = ones(lenData, 1, 'int8')*rawFlag; +flagsTested = ones(lenDataTested, 1, 'int8')*rawFlag; + +I = true(lenDataTested, 1); +I(end) = false; + +Ip1 = [false; I(1:end-1)]; + +presRelRef = (presRel(I) + presRel(Ip1))/2; + +SA = gsw_SA_from_SP(psal, presRel, lon, lat); +CT = gsw_CT_from_t(SA, temp, presRel); +rhoPotI = gsw_rho(SA(I), CT(I), presRelRef); % potential density referenced to the mid-point pressure +rhoPotIp1 = gsw_rho(SA(Ip1), CT(Ip1), presRelRef); + +deltaUpBottom = [rhoPotI - rhoPotIp1; NaN]; +deltaBottomUp = [NaN; rhoPotIp1 - rhoPotI]; + +threshold = ones(lenDataTested, 1) .* threshold; + +iDensInv = deltaUpBottom >= threshold | deltaBottomUp <= -threshold; + +flagsTested(iDensInv) = badFlag; +flagsTested(~iDensInv) = goodFlag; +flags(~iBadData) = flagsTested; + +if isMatrix + % we fold the vector back into a matrix + flags = reshape(flags, [len1, len2, len3]); +end + +sample_data.variables{tempIdx}.flags = flags; +sample_data.variables{psalIdx}.flags = flags; + +varChecked = {sample_data.variables{tempIdx}.name, ... + sample_data.variables{psalIdx}.name}; + +densIdx = getVar(sample_data.variables, 'DENS'); +if densIdx + sample_data.variables{densIdx}.flags = flags; + varChecked = [varChecked, {sample_data.variables{densIdx}.name}]; +end + +cndcIdx = getVar(sample_data.variables, 'CNDC'); +if cndcIdx + if any(ismember(sample_data.variables{cndcIdx}.flags, [goodFlag, probGoodFlag])) + sample_data.variables{cndcIdx}.flags = flags; + varChecked = [varChecked, {sample_data.variables{cndcIdx}.name}]; + end +end + +% write/update dataset QC parameters +writeDatasetParameter(sample_data.toolbox_input_file, currentQCtest, 'threshold', threshold); \ No newline at end of file diff --git a/AutomaticQC/imosDensityInversionSetQC.txt b/AutomaticQC/imosDensityInversionSetQC.txt new file mode 100644 index 000000000..192204adb --- /dev/null +++ b/AutomaticQC/imosDensityInversionSetQC.txt @@ -0,0 +1 @@ +threshold = 0.03 diff --git a/AutomaticQC/imosVerticalSpikeQC.m b/AutomaticQC/imosVerticalSpikeQC.m index f8c6a7b88..1b3bc469f 100644 --- a/AutomaticQC/imosVerticalSpikeQC.m +++ b/AutomaticQC/imosVerticalSpikeQC.m @@ -59,6 +59,12 @@ paramsLog = []; flags = []; +% this test only applies to profile mode +mode = readProperty('toolbox.mode'); +if ~strcmpi(mode, 'profile') + return; +end + % read all values from imosSpikeQC properties file values = readProperty('*', fullfile('AutomaticQC', 'imosVerticalSpikeQC.txt')); @@ -91,6 +97,7 @@ passFlag = imosQCFlag('good', qcSet, 'flag'); failFlag = imosQCFlag('probablyBad', qcSet, 'flag'); badFlag = imosQCFlag('bad', qcSet, 'flag'); + badFlags = [failFlag, badFlag]; paramsLog = ['threshold=' thresholdExpr{iParam}]; @@ -105,7 +112,7 @@ end % we don't consider already bad data in the current test - iBadData = sample_data.variables{k}.flags == badFlag; + iBadData = ismember(sample_data.variables{k}.flags, badFlags); dataTested = data(~iBadData); if isempty(dataTested), return; end diff --git a/Graph/DepthProfile/highlightDepthProfileGeneric.m b/Graph/DepthProfile/highlightDepthProfileGeneric.m index d9d428e44..5f9a73845 100644 --- a/Graph/DepthProfile/highlightDepthProfileGeneric.m +++ b/Graph/DepthProfile/highlightDepthProfileGeneric.m @@ -19,7 +19,8 @@ % data. If no data points lie within the highlight region, % an empty matrix is returned. % -% Author: Paul McCarthy +% Author: Paul McCarthy +% Contributor: Guillaume Galibert % % @@ -39,4 +40,48 @@ % along with this program. % If not, see . % -highlight = highlightTimeSeriesGeneric(region, data, variable, type); +narginchk(4, 4); + +if ~isnumeric(region) || ~isvector(region) || length(region) ~= 4 + error('region must be a numeric vector of length 4'); +end + +if ~ishandle(data), error('data must be a graphics handle'); end +if ~isstruct(variable), error('variable must be a struct'); end +if ~ischar(type), error('type must be a string'); end + +xdata = get(data(1), 'XData'); % data(1) retrieves the current graphic handles only, in case extra sample is selected +ydata = get(data(1), 'YData'); + +if iscell(xdata) + xdata = cell2mat(xdata)'; + ydata = cell2mat(ydata)'; +end + +if strcmp(type, 'alt') % right click + % figure out indices of all data points outside the X range but within + % the Y range + xidx = (xdata < region(1)) | (xdata > region(3)); + yidx = (ydata >= region(2)) & (ydata <= region(4)); +else + % figure out indices of all data points within the X and Y range + xidx = (xdata >= region(1)) & (xdata <= region(3)); + yidx = (ydata >= region(2)) & (ydata <= region(4)); +end + +% figure out indices of all the points to be highlighted +idx = xidx & yidx; + +if ~any(idx) + % return nothing if no points to plot + highlight = []; +else + % create the highlight + highlight = line(xdata(idx),ydata(idx), ... + 'UserData', idx, ... + 'Parent', gca, ... + 'LineStyle', 'none', ... + 'Marker', 'o', ... + 'MarkerEdgeColor', 'white', ... + 'MarkerFaceColor', 'white'); +end diff --git a/Graph/TimeSeries/highlightTimeSeriesGeneric.m b/Graph/TimeSeries/highlightTimeSeriesGeneric.m index 7dd1dfa43..fa1401990 100644 --- a/Graph/TimeSeries/highlightTimeSeriesGeneric.m +++ b/Graph/TimeSeries/highlightTimeSeriesGeneric.m @@ -55,20 +55,16 @@ ydata = cell2mat(ydata)'; end -% on right click highlight, only highlight -% unflagged data points in the region -if strcmp(type, 'alt') - f = variable.flags; - f = f == 0; - - xdata = xdata(f); - ydata = ydata(f); +if strcmp(type, 'alt') % right click + % figure out indices of all data points within the X range but outside the Y range + xidx = (xdata >= region(1)) & (xdata <= region(3)); + yidx = (ydata < region(2)) | (ydata > region(4)); +else + % figure out indices of all data points within the X and Y range + xidx = (xdata >= region(1)) & (xdata <= region(3)); + yidx = (ydata >= region(2)) & (ydata <= region(4)); end -% figure out indices of all data points within the range -xidx = (xdata >= region(1) & xdata <= region(3)); -yidx = (ydata >= region(2) & ydata <= region(4)); - % figure out indices of all the points to be highlighted idx = xidx & yidx; diff --git a/Graph/TimeSeries/highlightTimeSeriesTimeDepth.m b/Graph/TimeSeries/highlightTimeSeriesTimeDepth.m index 62d55d597..42e3a17d6 100644 --- a/Graph/TimeSeries/highlightTimeSeriesTimeDepth.m +++ b/Graph/TimeSeries/highlightTimeSeriesTimeDepth.m @@ -61,21 +61,18 @@ Ydata = repmat(ydata', [size(xdata), 1]); clear xdata ydata; -% figure out indices of all data points within the range -X = ((Xdata >= region(1)) & (Xdata <= region(3))); -Y = ((Ydata >= region(2)) & (Ydata <= region(4))); +if strcmp(type, 'alt') % right click + % figure out indices of all data points within the X range but outside the Y range + X = (Xdata >= region(1)) & (Xdata <= region(3)); + Y = (Ydata < region(2)) | (Ydata > region(4)); +else + % figure out indices of all data points within the X and Y range + X = (Xdata >= region(1)) & (Xdata <= region(3)); + Y = (Ydata >= region(2)) & (Ydata <= region(4)); +end idx = X & Y; clear X Y; -% on right click, only highlight unflagged points -if strcmp(type, 'alt') - - % see if any points in the region haven't been flagged - idx = (variable.flags == 0) & idx; - if ~any(any(variable.flags)), return; end - -end - if ~any(any(idx)) highlight = []; else diff --git a/Graph/TimeSeries/highlightTimeSeriesTimeFrequencyDirection.m b/Graph/TimeSeries/highlightTimeSeriesTimeFrequencyDirection.m index 42c57631a..6106b0377 100644 --- a/Graph/TimeSeries/highlightTimeSeriesTimeFrequencyDirection.m +++ b/Graph/TimeSeries/highlightTimeSeriesTimeFrequencyDirection.m @@ -50,21 +50,18 @@ Xdata = get(data, 'XData'); Ydata = get(data, 'YData'); -% figure out indices of all data points within the range -X = ((Xdata >= region(1)) & (Xdata <= region(3))); -Y = ((Ydata >= region(2)) & (Ydata <= region(4))); +if strcmp(type, 'alt') % right click + % figure out indices of all data points outside the X and Y range + X = ((Xdata < region(1)) | (Xdata > region(3))); + Y = ((Ydata < region(2)) | (Ydata > region(4))); +else + % figure out indices of all data points within the X and Y range + X = (Xdata >= region(1)) & (Xdata <= region(3)); + Y = (Ydata >= region(2)) & (Ydata <= region(4)); +end idx = X & Y; clear X Y; -% on right click, only highlight unflagged points -if strcmp(type, 'alt') - - % see if any points in the region haven't been flagged - idx = (variable.flags == 0) & idx; - if ~any(any(variable.flags)), return; end - -end - if ~any(any(idx)) highlight = []; else diff --git a/Graph/TimeSeries/parameters.txt b/Graph/TimeSeries/parameters.txt index 7ef4f8acc..e9e70537e 100644 --- a/Graph/TimeSeries/parameters.txt +++ b/Graph/TimeSeries/parameters.txt @@ -35,6 +35,7 @@ NORTEK_QC, TimeDepth VEL1, TimeDepth VEL2, TimeDepth VEL3, TimeDepth +VEL4, TimeDepth VDEN, TimeFrequency VDEV, TimeFrequency VDEP, TimeFrequency diff --git a/Graph/TimeSeries/setTimeSerieColorbarContextMenu.m b/Graph/TimeSeries/setTimeSerieColorbarContextMenu.m index f0104a5a8..291a00ba4 100644 --- a/Graph/TimeSeries/setTimeSerieColorbarContextMenu.m +++ b/Graph/TimeSeries/setTimeSerieColorbarContextMenu.m @@ -51,7 +51,7 @@ end switch upper(var.name(1:4)) - case {'UCUR', 'VCUR', 'WCUR', 'ECUR', 'VEL1', 'VEL2', 'VEL3'} % 0 centred parameters + case {'UCUR', 'VCUR', 'WCUR', 'ECUR', 'VEL1', 'VEL2', 'VEL3', 'VEL4'} % 0 centred parameters colormap(r_b); cbCLimRange('', '', 'auto, 0 centred', var.data); diff --git a/Graph/diagramMooring1DVarAgainstOther.m b/Graph/diagramMooring1DVarAgainstOther.m index 304578b52..103c6ca78 100644 --- a/Graph/diagramMooring1DVarAgainstOther.m +++ b/Graph/diagramMooring1DVarAgainstOther.m @@ -141,7 +141,7 @@ function diagramMooring1DVarAgainstOther(sample_data, varName, yAxisVarName, isQ if iVar && iYAxisVar && ... size(sample_data{iSort(i)}.variables{iVar}.data, 2) == 1 && ... size(sample_data{iSort(i)}.variables{iYAxisVar}.data, 2) == 1 && ... % we're only plotting 1D variables with depth variable but no current - all(~strncmpi(sample_data{iSort(i)}.variables{iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3'}, 4)) + all(~strncmpi(sample_data{iSort(i)}.variables{iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3', 'VEL4'}, 4)) iGood = true(size(sample_data{iSort(i)}.variables{iVar}.data)); if isQC %get time, depth and var QC information diff --git a/Graph/diagramMooring2DVarAgainstOther.m b/Graph/diagramMooring2DVarAgainstOther.m index b4e5cc750..498fb3cbb 100644 --- a/Graph/diagramMooring2DVarAgainstOther.m +++ b/Graph/diagramMooring2DVarAgainstOther.m @@ -184,7 +184,7 @@ function diagramMooring2DVarAgainstOther(sample_data, varName, yAxisVarName, isQ end elseif iVar && iYAxisVar && ... - any(strncmpi(sample_data{iSort(i)}.variables{iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3'}, 4)) && ... + any(strncmpi(sample_data{iSort(i)}.variables{iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3', 'VEL4'}, 4)) && ... size(sample_data{iSort(i)}.variables{iVar}.data, 2) == 1 % we're plotting current metre 1D variables with DEPTH variable. iGood = true(size(sample_data{iSort(i)}.variables{iVar}.data)); if isQC diff --git a/Graph/lineMooring1DVar.m b/Graph/lineMooring1DVar.m index 3834eef52..d94e06fdc 100644 --- a/Graph/lineMooring1DVar.m +++ b/Graph/lineMooring1DVar.m @@ -143,7 +143,7 @@ function lineMooring1DVar(sample_data, varName, isQC, saveToFile, exportDir) iVar = getVar(sample_data{iSort(i)}.(typeVar), varName); if iVar > 0 && size(sample_data{iSort(i)}.(typeVar){iVar}.data, 2) == 1 && ... % we're only plotting 1D variables but no current - all(~strncmpi(sample_data{iSort(i)}.(typeVar){iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3'}, 4)) + all(~strncmpi(sample_data{iSort(i)}.(typeVar){iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3', 'VEL4'}, 4)) if initiateFigure fileName = genIMOSFileName(sample_data{iSort(i)}, 'png'); visible = 'on'; diff --git a/Graph/pcolorMooring2DVar.m b/Graph/pcolorMooring2DVar.m index 3971da4c8..e611292e6 100644 --- a/Graph/pcolorMooring2DVar.m +++ b/Graph/pcolorMooring2DVar.m @@ -108,7 +108,7 @@ function pcolorMooring2DVar(sample_data, varName, isQC, saveToFile, exportDir) instrumentDesc{i} = [strrep(instrumentDesc{i}, '_', ' ') ' (' num2str(metaDepth(i)) 'm' instrumentSN ')']; switch varName(1:4) - case {'UCUR', 'VCUR', 'WCUR', 'VEL1', 'VEL2', 'VEL3'} % 0 centred parameters + case {'UCUR', 'VCUR', 'WCUR', 'VEL1', 'VEL2', 'VEL3', 'VEL4'} % 0 centred parameters cMap = 'r_b'; cType = 'centeredOnZero'; case {'CDIR', 'SSWD'} % directions [0; 360[ diff --git a/Graph/scatterMooring1DVarAgainstDepth.m b/Graph/scatterMooring1DVarAgainstDepth.m index 583f0699c..cbf549ce9 100644 --- a/Graph/scatterMooring1DVarAgainstDepth.m +++ b/Graph/scatterMooring1DVarAgainstDepth.m @@ -126,7 +126,7 @@ function scatterMooring1DVarAgainstDepth(sample_data, varName, isQC, saveToFile, iVar = getVar(sample_data{iSort(i)}.variables, varName); if iVar > 0 && size(sample_data{iSort(i)}.variables{iVar}.data, 2) == 1 && ... % we're only plotting 1D variables with depth variable but no current - all(~strncmpi(sample_data{iSort(i)}.variables{iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3'}, 4)) + all(~strncmpi(sample_data{iSort(i)}.variables{iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3', 'VEL4'}, 4)) iGood = true(size(sample_data{iSort(i)}.variables{iVar}.data)); if isQC %get time, depth and var QC information diff --git a/Graph/scatterMooring2DVarAgainstDepth.m b/Graph/scatterMooring2DVarAgainstDepth.m index 2c23bb054..5e06231f7 100644 --- a/Graph/scatterMooring2DVarAgainstDepth.m +++ b/Graph/scatterMooring2DVarAgainstDepth.m @@ -148,7 +148,7 @@ function scatterMooring2DVarAgainstDepth(sample_data, varName, isQC, saveToFile, end elseif iVar > 0 && ... - any(strncmpi(sample_data{iSort(i)}.variables{iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3'}, 4)) && ... + any(strncmpi(sample_data{iSort(i)}.variables{iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3', 'VEL4'}, 4)) && ... size(sample_data{iSort(i)}.variables{iVar}.data, 2) == 1 % we're plotting current metre 1D variables with DEPTH variable. iGood = true(size(sample_data{iSort(i)}.variables{iVar}.data)); if isQC @@ -182,7 +182,7 @@ function scatterMooring2DVarAgainstDepth(sample_data, varName, isQC, saveToFile, % define cMap, cLim and cType per parameter switch varName(1:4) - case {'UCUR', 'VCUR', 'WCUR', 'ECUR', 'VEL1', 'VEL2', 'VEL3'} % 0 centred parameters + case {'UCUR', 'VCUR', 'WCUR', 'ECUR', 'VEL1', 'VEL2', 'VEL3', 'VEL4'} % 0 centred parameters cMap = 'r_b'; cType = 'centeredOnZero'; CLim = [-max(abs(yLimMin), abs(yLimMax)) max(abs(yLimMax), abs(yLimMin))]; diff --git a/IMOS/imosParameters.txt b/IMOS/imosParameters.txt index f54b7d9fe..0ca6945ce 100644 --- a/IMOS/imosParameters.txt +++ b/IMOS/imosParameters.txt @@ -162,6 +162,7 @@ VDIRT, 1, sea_surface_wave_to_direction, VEL1, 0, sea_water_velocity_from_acoustic_beam_1, m s-1, , , V, 999999.0, -10.0, 10.0, float VEL2, 0, sea_water_velocity_from_acoustic_beam_2, m s-1, , , V, 999999.0, -10.0, 10.0, float VEL3, 0, sea_water_velocity_from_acoustic_beam_3, m s-1, , , V, 999999.0, -10.0, 10.0, float +VEL4, 0, sea_water_velocity_from_acoustic_beam_4, m s-1, , , V, 999999.0, -10.0, 10.0, float VOLT, 0, voltage, V, , , E, 999999.0, -100.0, 100.0, float VSF412, 0, volume_scattering_function_of_radiative_flux_in_sea_water_for_wavelength_412nm, m-1 sr-1, , , E, 999999.0, , , float VSF440, 0, volume_scattering_function_of_radiative_flux_in_sea_water_for_wavelength_440nm, m-1 sr-1, , , E, 999999.0, , , float diff --git a/Java/ddb.jar b/Java/ddb.jar index 0265c8faf..ddf76c603 100644 Binary files a/Java/ddb.jar and b/Java/ddb.jar differ diff --git a/Java/src/org/imos/ddb/JDBCDDB.java b/Java/src/org/imos/ddb/JDBCDDB.java index eb26b3c60..1a49729bf 100644 --- a/Java/src/org/imos/ddb/JDBCDDB.java +++ b/Java/src/org/imos/ddb/JDBCDDB.java @@ -100,6 +100,8 @@ public ArrayList executeQuery( throws Exception { Connection conn = null; + Statement stmt = null; + ResultSet rs = null; ArrayList results = null; try { @@ -126,8 +128,6 @@ public ArrayList executeQuery( } //execute the query - Statement stmt = null; - ResultSet rs = null; try { stmt = conn.createStatement(); rs = stmt.executeQuery(query); @@ -177,7 +177,13 @@ public ArrayList executeQuery( } //always close db connection - finally {try {conn.close();} catch (Exception e) {}} + finally { + try { + rs.close(); + stmt.close(); + conn.close(); + } + catch (Exception e) {}} return results; } diff --git a/Parser/readAD2CPBinary.m b/Parser/readAD2CPBinary.m index 7e4356879..3b693748b 100644 --- a/Parser/readAD2CPBinary.m +++ b/Parser/readAD2CPBinary.m @@ -77,10 +77,12 @@ [sect, len] = readSection(filename, data, dIdx, cpuEndianness); if ~isempty(sect) + % grouping by Id, Version and Size is our only option, there is no + % information on whether record type is XXX or Alt_XXX like in .mat files. if isfield(sect.Data, 'Version') - curField = ['Id' sprintf('%s', sect.Header.Id) 'Version' sprintf('%d', sect.Data.Version)]; + curField = ['Id' sprintf('%s', sect.Header.Id) '_Version' sprintf('%d', sect.Data.Version) '_Size' sprintf('%u', sect.Header.DataSize)]; else - curField = ['Id' sprintf('%s', sect.Header.Id)]; + curField = ['Id' sprintf('%s', sect.Header.Id) '_Size' sprintf('%u', sect.Header.DataSize)]; end theFieldNames = fieldnames(sect); nField = length(theFieldNames); @@ -150,20 +152,46 @@ % read the section in id = dec2hex(data(idx+2)); switch id - case '15' + case '15' % Burst Data Record [sect, ~, off] = readBurstAverage (data, idx, cpuEndianness); % 0x15 - case '16' + + case '16' % Average Data Record [sect, ~, off] = readBurstAverage (data, idx, cpuEndianness); % 0x16 - case '17' + + case '17' % Bottom Track Data Record [sect, ~, off] = readBottomTrack (data, idx, cpuEndianness); % 0x17 - case '18' + + case '18' % Interleaved Burst Data Record (beam 5) [sect, ~, off] = readBurstAverage (data, idx, cpuEndianness); % 0x18 - case 'A0' + + case '1A' % Burst Altimeter Raw Record + [sect, ~, off] = readBurstAverage (data, idx, cpuEndianness); % 0x1A + + case '1B' % DVL Bottom Track Record + [sect, ~, off] = readBottomTrack (data, idx, cpuEndianness); % 0x1B + + case '1C' % Echo Sounder Record + disp('Echo Sounder Record not supported'); + off = len; + + case '1D' % DVL Water Track Record + disp('DVL Water Track Record not supported'); + off = len; + + case '1E' % Altimeter Record + [sect, ~, off] = readBurstAverage (data, idx, cpuEndianness); % 0x1E + + case '1F' % Avg Altimeter Raw Record + [sect, ~, off] = readBurstAverage (data, idx, cpuEndianness); % 0x1F + + case 'A0' % String Data Record, eg. GPS NMEA data, comment from the FWRITE command [sect, ~, off] = readString (data, idx, cpuEndianness); % 0xA0 + otherwise disp('Unknown section type'); disp(['ID : hex ' id ' at offset ' num2str(idx+2)]); off = len; + end if isempty(sect), return; end @@ -285,8 +313,11 @@ iEndCell = 10; iStartBeam = 13; iEndBeam = 16; +iStartCoord = 11; +iEndCoord = 12; sect.nCells = bin2dec(sect.Beams_CoordSys_Cells(end-iEndCell+1:end-iStartCell+1)); sect.nBeams = bin2dec(sect.Beams_CoordSys_Cells(end-iEndBeam+1:end-iStartBeam+1)); +sect.coordSys = bin2dec(sect.Beams_CoordSys_Cells(end-iEndCoord+1:end-iStartCoord+1)); sect.CellSize = bytecast(data(idx+30:idx+31), 'L', 'uint16', cpuEndianness); % 1 mm sect.Blanking = bytecast(data(idx+32:idx+33), 'L', 'uint16', cpuEndianness); % 1 mm @@ -343,8 +374,11 @@ iEndCell = 10; iStartBeam = 13; iEndBeam = 16; +iStartCoord = 11; +iEndCoord = 12; sect.nCells = bin2dec(sect.Beams_CoordSys_Cells(end-iEndCell+1:end-iStartCell+1)); sect.nBeams = bin2dec(sect.Beams_CoordSys_Cells(end-iEndBeam+1:end-iStartBeam+1)); +sect.coordSys = bin2dec(sect.Beams_CoordSys_Cells(end-iEndCoord+1:end-iStartCoord+1)); sect.CellSize = bytecast(data(idx+36:idx+37), 'L', 'uint16', cpuEndianness); % 1 mm sect.Blanking = bytecast(data(idx+38:idx+39), 'L', 'uint16', cpuEndianness); % 1 mm @@ -356,7 +390,7 @@ sect.AccRawX = bytecast(data(idx+50:idx+51), 'L', 'int16', cpuEndianness); % accelerometer Raw, X axis value in last measurement interval (16384 = 1.0) sect.AccRawY = bytecast(data(idx+52:idx+53), 'L', 'int16', cpuEndianness); sect.AccRawZ = bytecast(data(idx+54:idx+55), 'L', 'int16', cpuEndianness); -sect.AmbiguityVelocity = bytecast(data(idx+56:idx+57), 'L', 'uint16', cpuEndianness); % 0.1mm/s ; corrected for sound velocity +sect.AmbiguityVel = bytecast(data(idx+56:idx+57), 'L', 'uint16', cpuEndianness); % 0.1mm/s ; corrected for sound velocity sect.DatasetDesc = dec2bin(bytecast(data(idx+58:idx+59), 'L', 'uint16', cpuEndianness), 16); sect.TransmitEnergy = bytecast(data(idx+60:idx+61), 'L', 'uint16', cpuEndianness); sect.VelocityScaling = bytecast(data(idx+62), 'L', 'int8', cpuEndianness); @@ -416,7 +450,7 @@ iEndCoord = 12; sect.nCells = bin2dec(sect.Beams_CoordSys_Cells(end-iEndCell+1:end-iStartCell+1)); sect.nBeams = bin2dec(sect.Beams_CoordSys_Cells(end-iEndBeam+1:end-iStartBeam+1)); -sect.coordSys = sect.Beams_CoordSys_Cells(end-iEndCoord+1:end-iStartCoord+1); +sect.coordSys = bin2dec(sect.Beams_CoordSys_Cells(end-iEndCoord+1:end-iStartCoord+1)); sect.CellSize = bytecast(data(idx+32:idx+33), 'L', 'uint16', cpuEndianness); % 1 mm sect.Blanking = bytecast(data(idx+34:idx+35), 'L', 'uint16', cpuEndianness); % 1 mm @@ -429,7 +463,7 @@ sect.AccRawX = bytecast(data(idx+46:idx+47), 'L', 'int16', cpuEndianness); % accelerometer Raw, X axis value in last measurement interval (16384 = 1.0) sect.AccRawY = bytecast(data(idx+48:idx+49), 'L', 'int16', cpuEndianness); sect.AccRawZ = bytecast(data(idx+50:idx+51), 'L', 'int16', cpuEndianness); -sect.AmbiguityVelocity = bytecast(data(idx+52:idx+53), 'L', 'uint16', cpuEndianness); % 10^(Velocity scaling) m/s ; corrected for sound velocity +sect.AmbiguityVel = bytecast(data(idx+52:idx+53), 'L', 'uint16', cpuEndianness); % 10^(Velocity scaling) m/s ; corrected for sound velocity sect.DatasetDesc = dec2bin(bytecast(data(idx+54:idx+55), 'L', 'uint16', cpuEndianness), 16); sect.TransmitEnergy = bytecast(data(idx+56:idx+57), 'L', 'uint16', cpuEndianness); sect.VelocityScaling = bytecast(data(idx+58), 'L', 'int8', cpuEndianness); @@ -521,7 +555,7 @@ sect.Data.AccRawX = bytecast(data(idx+46:idx+47), 'L', 'int16', cpuEndianness); % accelerometer Raw, X axis value in last measurement interval (16384 = 1.0) sect.Data.AccRawY = bytecast(data(idx+48:idx+49), 'L', 'int16', cpuEndianness); sect.Data.AccRawZ = bytecast(data(idx+50:idx+51), 'L', 'int16', cpuEndianness); -sect.Data.AmbiguityVelocity = bytecast(data(idx+52:idx+53), 'L', 'uint16', cpuEndianness); % 10^(Velocity scaling) m/s ; corrected for sound velocity +sect.Data.AmbiguityVel = bytecast(data(idx+52:idx+53), 'L', 'uint16', cpuEndianness); % 10^(Velocity scaling) m/s ; corrected for sound velocity sect.Data.DatasetDesc = dec2bin(bytecast(data(idx+54:idx+55), 'L', 'uint16', cpuEndianness), 16); sect.Data.TransmitEnergy = bytecast(data(idx+56:idx+57), 'L', 'uint16', cpuEndianness); sect.Data.VelocityScaling = bytecast(data(idx+58), 'L', 'int8', cpuEndianness); diff --git a/Parser/signatureParse.m b/Parser/signatureParse.m index f535b360d..c42d62bc7 100644 --- a/Parser/signatureParse.m +++ b/Parser/signatureParse.m @@ -37,7 +37,7 @@ % only one file supported filename = filename{1}; -[path, name, ext] = fileparts(filename); +[~, ~, ext] = fileparts(filename); isMagBias = false; magDec = 0; @@ -47,87 +47,183 @@ % read in all of the structures in the raw file structures = readAD2CPBinary(filename); + acquisitionMode = fieldnames(structures); + instrumentModel = ''; - if isfield(structures, 'IdA0') + % look for the string data record A0 + iA0 = strncmpi('IdA0', acquisitionMode, 4); + if any(iA0) % this is a string data record - if structures.IdA0.Data.Id == 16 + if structures.(acquisitionMode{iA0}).Data.Id == 16 % looking for instrument model - instrumentModel = regexp(structures.IdA0.Data.String, 'Signature[0-9]*', 'match', 'once'); + instrumentModel = regexp(structures.(acquisitionMode{iA0}).Data.String, 'Signature[0-9]*', 'match', 'once'); % check for magnetic declination - stringCell = textscan(structures.IdA0.Data.String, '%s', 'Delimiter', ','); + stringCell = textscan(structures.(acquisitionMode{iA0}).Data.String, '%s', 'Delimiter', ','); iMagDec = strncmp('DECL=', stringCell{1}, 5); if any(iMagDec) magDec = stringCell{1}{iMagDec}; magDec = textscan(magDec, 'DECL=%f'); magDec = magDec{1}; end + + structures = rmfield(structures, acquisitionMode{iA0}); + acquisitionMode(iA0) = []; end end - if isfield(structures, 'Id15Version3') - % this is a burst data record - dataRecordType = 'Id15Version3'; - elseif isfield(structures, 'Id16Version3') - % this is an average data record - dataRecordType = 'Id16Version3'; + % look for burst/average data record version 3 + i15 = strncmpi('Id15_Version3', acquisitionMode, 13); + i16 = strncmpi('Id16_Version3', acquisitionMode, 13); + iSupported = i15 | i16; + if any(~iSupported) + dataRecordNamesNotSupported = acquisitionMode(~iSupported); + acquisitionMode(~iSupported) = []; + for i=1:length(dataRecordNamesNotSupported) + structures = rmfield(structures, dataRecordNamesNotSupported{i}); + disp(['Warning : data record ' dataRecordNamesNotSupported{i} ' is not supported.']); + end + end + if ~any(iSupported) + error('Not a single burst/average data record version 3 format found (only format supported).'); else - % string/bottom track/interleaved data record is not supported - error('Only burst/average data record version 3 format is supported'); + serialNumber = num2str(structures.(acquisitionMode{1}).Data(1).SerialNumber); end - nSamples = length(structures.(dataRecordType).Header); - nCells = structures.(dataRecordType).Data(1).nCells; - - serialNumber = num2str(unique(vertcat(structures.(dataRecordType).Data.SerialNumber))); - - cellSize = unique(vertcat(structures.(dataRecordType).Data.CellSize))*0.001; % m - blankDist = unique(vertcat(structures.(dataRecordType).Data.Blanking))*0.001; % m - if length(cellSize) > 1, error('Multiple cell sizes/blanking distance not supported'); end - - % retrieve sample data - Time = vertcat(structures.(dataRecordType).Data.Time); - speedOfSound = vertcat(structures.(dataRecordType).Data.SpeedOfSound)*0.1; % m/s - Temperature = vertcat(structures.(dataRecordType).Data.Temperature)*0.01; % degree Celsius - Pressure = vertcat(structures.(dataRecordType).Data.Pressure)*0.001; % dBar - Heading = vertcat(structures.(dataRecordType).Data.Heading)*0.01; % deg - Pitch = vertcat(structures.(dataRecordType).Data.Pitch)*0.01; % deg - Roll = vertcat(structures.(dataRecordType).Data.Roll)*0.01; % deg - Battery = vertcat(structures.(dataRecordType).Data.BatteryVoltage)*0.1; % Volt - Status = vertcat(structures.(dataRecordType).Data.Status); - velocity = cat(3, structures.(dataRecordType).Data.VelocityData); - velocityScaling = repmat(10.^vertcat(structures.(dataRecordType).Data.VelocityScaling), 1, nCells); - Velocity_E = squeeze(velocity(1,:,:))'.*velocityScaling; % m/s - Velocity_N = squeeze(velocity(2,:,:))'.*velocityScaling; - Velocity_U = squeeze(velocity(3,:,:))'.*velocityScaling; - Velocity_U2 = squeeze(velocity(4,:,:))'.*velocityScaling; - amplitude = cat(3, structures.(dataRecordType).Data.AmplitudeData); - Backscatter1 = squeeze(amplitude(1,:,:))'; % count - Backscatter2 = squeeze(amplitude(2,:,:))'; - Backscatter3 = squeeze(amplitude(3,:,:))'; - Backscatter4 = squeeze(amplitude(4,:,:))'; - correlation = cat(3, structures.(dataRecordType).Data.CorrelationData); - Correlation1 = squeeze(correlation(1,:,:))'; % percent [0 - 100] - Correlation2 = squeeze(correlation(2,:,:))'; - Correlation3 = squeeze(correlation(3,:,:))'; - Correlation4 = squeeze(correlation(4,:,:))'; - clear structures velocity velocityScaling amplitude correlation; + data = struct; + nDataset = length(acquisitionMode); + for i=1:nDataset + data.(acquisitionMode{i}).nSamples = length(structures.(acquisitionMode{i}).Header); + data.(acquisitionMode{i}).nCells = unique(vertcat(structures.(acquisitionMode{i}).Data.nCells)); + data.(acquisitionMode{i}).nBeams = unique(vertcat(structures.(acquisitionMode{i}).Data.nBeams)); + data.(acquisitionMode{i}).coordSys = unique(vertcat(structures.(acquisitionMode{i}).Data.coordSys)); + data.(acquisitionMode{i}).cellSize = unique(vertcat(structures.(acquisitionMode{i}).Data.CellSize))*0.001; % m + data.(acquisitionMode{i}).blankDist = unique(vertcat(structures.(acquisitionMode{i}).Data.Blanking))*0.001; % m + + if length(data.(acquisitionMode{i}).nCells) > 1, error('Multiple nCells values not supported'); end + if length(data.(acquisitionMode{i}).nBeams) > 1, error('Multiple nBeams values not supported'); end + if length(data.(acquisitionMode{i}).coordSys) > 1, error('Multiple coordSys values not supported'); end + if length(data.(acquisitionMode{i}).cellSize) > 1, error('Multiple cellSize values not supported'); end + if length(data.(acquisitionMode{i}).blankDist) > 1, error('Multiple blankDist values not supported'); end + + % retrieve sample data + data.(acquisitionMode{i}).Time = vertcat(structures.(acquisitionMode{i}).Data.Time); + data.(acquisitionMode{i}).speedOfSound = vertcat(structures.(acquisitionMode{i}).Data.SpeedOfSound)*0.1; % m/s + data.(acquisitionMode{i}).Temperature = vertcat(structures.(acquisitionMode{i}).Data.Temperature)*0.01; % degree Celsius + data.(acquisitionMode{i}).Pressure = vertcat(structures.(acquisitionMode{i}).Data.Pressure)*0.001; % dBar + data.(acquisitionMode{i}).Heading = vertcat(structures.(acquisitionMode{i}).Data.Heading)*0.01; % deg + data.(acquisitionMode{i}).Pitch = vertcat(structures.(acquisitionMode{i}).Data.Pitch)*0.01; % deg + data.(acquisitionMode{i}).Roll = vertcat(structures.(acquisitionMode{i}).Data.Roll)*0.01; % deg + data.(acquisitionMode{i}).Battery = vertcat(structures.(acquisitionMode{i}).Data.BatteryVoltage)*0.1; % Volt + data.(acquisitionMode{i}).Status = vertcat(structures.(acquisitionMode{i}).Data.Status); + data.(acquisitionMode{i}).Error = vertcat(structures.(acquisitionMode{i}).Data.Error); % error codes for each cell of a velocity profile inferred from the beams. 0=good; otherwise error. See http://www.nortek-as.com/en/knowledge-center/forum/waves/20001875?b_start=0#769595815 + data.(acquisitionMode{i}).AmbiguityVel = vertcat(structures.(acquisitionMode{i}).Data.AmbiguityVel); + data.(acquisitionMode{i}).TransmitEnergy = vertcat(structures.(acquisitionMode{i}).Data.TransmitEnergy); + data.(acquisitionMode{i}).NominalCorrelation = vertcat(structures.(acquisitionMode{i}).Data.NominalCorrelation); + + % support velocity data in ENU or beam coordinates + data.(acquisitionMode{i}).isVelocityData = false; + if isfield(structures.(acquisitionMode{i}).Data, 'VelocityData') + switch data.(acquisitionMode{i}).coordSys + case 0 % 0 is ENU, 1 is XYZ and 2 is BEAM + velocity = cat(3, structures.(acquisitionMode{i}).Data.VelocityData); + velocityScaling = repmat(10.^vertcat(structures.(acquisitionMode{i}).Data.VelocityScaling), 1, data.(acquisitionMode{i}).nCells); + data.(acquisitionMode{i}).Velocity_E = squeeze(velocity(1,:,:))'.*velocityScaling; % m/s + data.(acquisitionMode{i}).Velocity_N = squeeze(velocity(2,:,:))'.*velocityScaling; + data.(acquisitionMode{i}).Velocity_U = squeeze(velocity(3,:,:))'.*velocityScaling; + data.(acquisitionMode{i}).Velocity_U2 = squeeze(velocity(4,:,:))'.*velocityScaling; + clear velocity velocityScaling; + + data.(acquisitionMode{i}).isVelocityData = true; + data.(acquisitionMode{i}).coordSys = 'ENU'; + + case 2 + velocity = cat(3, structures.(acquisitionMode{i}).Data.VelocityData); + velocityScaling = repmat(10.^vertcat(structures.(acquisitionMode{i}).Data.VelocityScaling), 1, data.(acquisitionMode{i}).nCells); + data.(acquisitionMode{i}).Velocity_1 = squeeze(velocity(1,:,:))'.*velocityScaling; % m/s + data.(acquisitionMode{i}).Velocity_2 = squeeze(velocity(2,:,:))'.*velocityScaling; + data.(acquisitionMode{i}).Velocity_3 = squeeze(velocity(3,:,:))'.*velocityScaling; + data.(acquisitionMode{i}).Velocity_4 = squeeze(velocity(4,:,:))'.*velocityScaling; + clear velocity velocityScaling; + + data.(acquisitionMode{i}).isVelocityData = true; + data.(acquisitionMode{i}).coordSys = 'BEAM'; + + end + end + + if data.(acquisitionMode{i}).isVelocityData + amplitude = cat(3, structures.(acquisitionMode{i}).Data.AmplitudeData); + data.(acquisitionMode{i}).Backscatter1 = squeeze(amplitude(1,:,:))'; % count + data.(acquisitionMode{i}).Backscatter2 = squeeze(amplitude(2,:,:))'; + data.(acquisitionMode{i}).Backscatter3 = squeeze(amplitude(3,:,:))'; + data.(acquisitionMode{i}).Backscatter4 = squeeze(amplitude(4,:,:))'; + clear amplitude; + + correlation = cat(3, structures.(acquisitionMode{i}).Data.CorrelationData); + data.(acquisitionMode{i}).Correlation1 = squeeze(correlation(1,:,:))'; % percent [0 - 100] + data.(acquisitionMode{i}).Correlation2 = squeeze(correlation(2,:,:))'; + data.(acquisitionMode{i}).Correlation3 = squeeze(correlation(3,:,:))'; + data.(acquisitionMode{i}).Correlation4 = squeeze(correlation(4,:,:))'; + clear correlation; + end + + data.(acquisitionMode{i}).isAltimeterData = false; + if isfield(structures.(acquisitionMode{i}).Data, 'AltimeterDistance') + data.(acquisitionMode{i}).AltimeterDistanceLE = vertcat(structures.(acquisitionMode{i}).Data.AltimeterDistance); + data.(acquisitionMode{i}).AltimeterQualityLE = vertcat(structures.(acquisitionMode{i}).Data.AltimeterQuality); +% data.(acquisitionMode{i}).AltimeterStatus = vertcat(structures.(acquisitionMode{i}).Data.AltimeterStatus); % flags for when pitch or roll is > 5 or > 10. We already have pitch and roll data... + + data.(acquisitionMode{i}).isAltimeterData = true; + end + + data.(acquisitionMode{i}).isASTData = false; + if isfield(structures.(acquisitionMode{i}).Data, 'ASTDistance') + data.(acquisitionMode{i}).AltimeterDistanceAST = vertcat(structures.(acquisitionMode{i}).Data.ASTDistance); + data.(acquisitionMode{i}).AltimeterQualityAST = vertcat(structures.(acquisitionMode{i}).Data.ASTQuality); + data.(acquisitionMode{i}).AltimeterTimeOffsetAST = vertcat(structures.(acquisitionMode{i}).Data.ASTOffset100uSec); + data.(acquisitionMode{i}).AltimeterPressure = vertcat(structures.(acquisitionMode{i}).Data.ASTPressure); + + data.(acquisitionMode{i}).isASTData = true; + end + end + clear structures; case '.mat' + [~, ~, cpuEndianness] = computer; + % read in all of the structures in the file structures = load(filename); dataFields = fieldnames(structures.Data); % look for potentially more data in other .mat files [path, name, ext] = fileparts(filename); - iMatFile = 1; + + isAvgd = false; + if any(strfind(name, '_avgd')) + isAvgd = true; + end + + if strcmpi(name(end-1:end), '_1') + % we assume the first file is name_1.mat so the next one would + % be name_2.mat + name(end-1:end) = []; + iMatFile = 2; + else + % we assume the first file is name.mat so the next one would be + % name_1.mat + iMatFile = 1; + end nextFilename = fullfile(path, [name '_' num2str(iMatFile) ext]); + while exist(nextFilename, 'file') nextStructures = load(nextFilename); % for each field of the structure we append new data to the % existing ones for i=1:length(dataFields) - structures.Data.(dataFields{i}) = [structures.Data.(dataFields{i}); nextStructures.Data.(dataFields{i})]; + if isfield(nextStructures.Data, dataFields{i}) + structures.Data.(dataFields{i}) = [structures.Data.(dataFields{i}); nextStructures.Data.(dataFields{i})]; + end end iMatFile = iMatFile + 1; @@ -135,88 +231,360 @@ clear nextStructures; end - % investigate which mode has been set for acquisition - if strcmpi(structures.Config.Plan_BurstEnabled, 'TRUE') - acquisitionMode = 'Burst'; - end - - if strcmpi(structures.Config.Plan_AverageEnabled, 'TRUE') - acquisitionMode = 'Average'; - end - - % we need to sort the data by time just in case - [Time, iSort] = sort(structures.Data.([acquisitionMode '_Time'])); - for i=1:length(dataFields) - structures.Data.(dataFields{i}) = structures.Data.(dataFields{i})(iSort, :); - end - - nSamples = length(Time); - nCells = double(structures.Config.([acquisitionMode '_NCells'])); - serialNumber = num2str(structures.Config.SerialNo); instrumentModel = structures.Config.InstrumentName; - - cellSize = structures.Config.([acquisitionMode '_CellSize']); - blankDist = structures.Config.([acquisitionMode '_BlankingDistance']); - + % check for magnetic declination magDec = structures.Config.Declination; - [~, ~, cpuEndianness] = computer; - - Status = dec2bin(bytecast(structures.Data.([acquisitionMode '_Status']), 'L', 'uint32', cpuEndianness), 32); - % Error = structures.Data.Average_Error; % error codes for each cell of a velocity profile inferred from the beams. 0=good; otherwise error. See http://www.nortek-as.com/en/knowledge-center/forum/waves/20001875?b_start=0#769595815 - % AmbiguityVel = structures.Data.Average_AmbiguityVel; - - if isfield(structures.Data, [acquisitionMode '_Velocity_ENU']) - Velocity_E = squeeze(structures.Data.([acquisitionMode '_ENU'])(:, 1, :)); - Velocity_N = squeeze(structures.Data.([acquisitionMode '_ENU'])(:, 2, :)); - Velocity_U = squeeze(structures.Data.([acquisitionMode '_ENU'])(:, 3, :)); - Velocity_U2 = squeeze(structures.Data.([acquisitionMode '_ENU'])(:, 4, :)); - else - Velocity_E = structures.Data.([acquisitionMode '_VelEast']); - Velocity_N = structures.Data.([acquisitionMode '_VelNorth']); - Velocity_U = structures.Data.([acquisitionMode '_VelUp']); - Velocity_U2 = structures.Data.([acquisitionMode '_VelUp2']); + % investigate which mode(s) has been set for acquisition + acquisitionMode = {}; + configFields = fieldnames(structures.Config); + if strcmpi(structures.Config.Plan_BurstEnabled, 'True') && ~isAvgd + acquisitionMode{end+1} = 'Burst_'; end - if isfield(structures.Data, [acquisitionMode '_Amplitude_Beam']) - Backscatter1 = squeeze(structures.Data.([acquisitionMode '_Amplitude_Beam'])(:, 1, :))*2; % looks like the .mat format is giving raw counts * 0.5 by default - Backscatter2 = squeeze(structures.Data.([acquisitionMode '_Amplitude_Beam'])(:, 2, :))*2; - Backscatter3 = squeeze(structures.Data.([acquisitionMode '_Amplitude_Beam'])(:, 3, :))*2; - Backscatter4 = squeeze(structures.Data.([acquisitionMode '_Amplitude_Beam'])(:, 4, :))*2; - else - Backscatter1 = structures.Data.([acquisitionMode '_AmpBeam1'])*2; - Backscatter2 = structures.Data.([acquisitionMode '_AmpBeam2'])*2; - Backscatter3 = structures.Data.([acquisitionMode '_AmpBeam3'])*2; - Backscatter4 = structures.Data.([acquisitionMode '_AmpBeam4'])*2; + if strcmpi(structures.Config.Plan_AverageEnabled, 'True') + acquisitionMode{end+1} = 'Average_'; end - if isfield(structures.Data, [acquisitionMode '_Correlation_Beam']) - Correlation1 = squeeze(structures.Data.([acquisitionMode '_Correlation_Beam'])(:, 1, :)); - Correlation2 = squeeze(structures.Data.([acquisitionMode '_Correlation_Beam'])(:, 2, :)); - Correlation3 = squeeze(structures.Data.([acquisitionMode '_Correlation_Beam'])(:, 3, :)); - Correlation4 = squeeze(structures.Data.([acquisitionMode '_Correlation_Beam'])(:, 4, :)); - else - Correlation1 = structures.Data.([acquisitionMode '_CorBeam1']); - Correlation2 = structures.Data.([acquisitionMode '_CorBeam2']); - Correlation3 = structures.Data.([acquisitionMode '_CorBeam3']); - Correlation4 = structures.Data.([acquisitionMode '_CorBeam4']); + if strcmpi(structures.Config.Alt_Plan_BurstEnabled, 'True') && ~isAvgd + % we need to compare structures.Config.Plan_xxx VS structures.Config.Alt_Plan_xxx + % to see if we can put together the two datasets + configPlanFields = configFields(strncmpi('Plan_', configFields, 5) & ~strncmpi('Plan_Average', configFields, 12)); + configAltPlanFields = configFields(strncmpi('Alt_Plan_', configFields, 9) & ~strncmpi('Alt_Plan_Average', configFields, 16)); + + nPlanFields = length(configPlanFields); + samePlan = false(nPlanFields, 1); + for i=1:nPlanFields + if ischar(structures.Config.(configPlanFields{i})) + samePlan(i) = strcmp(structures.Config.(configPlanFields{i}), structures.Config.(configAltPlanFields{i})); + else + samePlan(i) = structures.Config.(configPlanFields{i}) == structures.Config.(configAltPlanFields{i}); + end + end + + configBurstFields = configFields(strncmpi('Burst_', configFields, 6)); + configAltBurstFields = configFields(strncmpi('Alt_Burst_', configFields, 10)); + + nBurstFields = length(configBurstFields); + sameBurst = false(nBurstFields, 1); + for i=1:nBurstFields + if ischar(structures.Config.(configBurstFields{i})) + sameBurst(i) = strcmp(structures.Config.(configBurstFields{i}), structures.Config.(configAltBurstFields{i})); + else + sameBurst(i) = all(all(structures.Config.(configBurstFields{i}) == structures.Config.(configAltBurstFields{i}))); + end + end + + if all(samePlan) && all(sameBurst) + % can be added to former dataset + [~, iBurst] = find(strcmp(acquisitionMode, 'Burst_')); + acquisitionMode{end+1, iBurst} = 'Alt_Burst_'; + else + % distinct dataset + acquisitionMode{1, end+1} = 'Alt_Burst_'; + end + end + if strcmpi(structures.Config.Alt_Plan_AverageEnabled, 'True') + % we need to compare structures.Config.Plan_xxx VS structures.Config.Alt_Plan_xxx + % to see if we can put together the two datasets + configPlanFields = configFields(strncmpi('Plan_', configFields, 5) & ~strncmpi('Plan_Burst', configFields, 10)); + configAltPlanFields = configFields(strncmpi('Alt_Plan_', configFields, 9) & ~strncmpi('Alt_Plan_Burst', configFields, 14)); + + nPlanFields = length(configPlanFields); + samePlan = false(nPlanFields, 1); + for i=1:nPlanFields + if ischar(structures.Config.(configPlanFields{i})) + samePlan(i) = strcmp(structures.Config.(configPlanFields{i}), structures.Config.(configAltPlanFields{i})); + else + samePlan(i) = structures.Config.(configPlanFields{i}) == structures.Config.(configAltPlanFields{i}); + end + end + + configAvgFields = configFields(strncmpi('Average_', configFields, 8)); + configAltAvgFields = configFields(strncmpi('Alt_Average_', configFields, 12)); + + nAvgFields = length(configAvgFields); + sameAvg = false(nAvgFields, 1); + for i=1:nAvgFields + if ischar(structures.Config.(configAvgFields{i})) + sameAvg(i) = strcmp(structures.Config.(configAvgFields{i}), structures.Config.(configAltAvgFields{i})); + else + sameAvg(i) = all(all(structures.Config.(configAvgFields{i}) == structures.Config.(configAltAvgFields{i}))); + end + end + + if all(samePlan) && all(sameAvg) + % can be added to former dataset + [~, iAverage] = find(strcmp(acquisitionMode, 'Average_')); + acquisitionMode{end+1, iAverage} = 'Alt_Average_'; + else + % distinct dataset + acquisitionMode{1, end+1} = 'Alt_Average_'; + end + end +% if strcmpi(structures.Config.Burst_Altimeter, 'True') +% acquisitionMode{end+1} = 'BurstRawAltimeter_'; +% end +% if strcmpi(structures.Config.Alt_Burst_Altimeter, 'True') +% acquisitionMode{end+1} = 'Alt_BurstRawAltimeter_'; +% end + + data = struct; + nDataset = size(acquisitionMode, 2); + for i=1:nDataset + data.(acquisitionMode{1, i}).Time = structures.Data.([acquisitionMode{1, i} 'Time']); + + data.(acquisitionMode{1, i}).instrument_sample_interval = 1; + + data.(acquisitionMode{1, i}).nSamples = length(data.(acquisitionMode{1, i}).Time); + data.(acquisitionMode{1, i}).nCells = double(structures.Config.([acquisitionMode{1, i} 'NCells'])); + data.(acquisitionMode{1, i}).nBeams = double(structures.Config.([acquisitionMode{1, i} 'NBeams'])); + + data.(acquisitionMode{1, i}).cellSize = structures.Config.([acquisitionMode{1, i} 'CellSize']); + data.(acquisitionMode{1, i}).blankDist = structures.Config.([acquisitionMode{1, i} 'BlankingDistance']); + + data.(acquisitionMode{1, i}).Status = dec2bin(bytecast(structures.Data.([acquisitionMode{1, i} 'Status']), 'L', 'uint32', cpuEndianness), 32); + data.(acquisitionMode{1, i}).Error = structures.Data.([acquisitionMode{1, i} 'Error']); % error codes for each cell of a velocity profile inferred from the beams. 0=good; otherwise error. See http://www.nortek-as.com/en/knowledge-center/forum/waves/20001875?b_start=0#769595815 + data.(acquisitionMode{1, i}).AmbiguityVel = structures.Data.([acquisitionMode{1, i} 'AmbiguityVel']); + data.(acquisitionMode{1, i}).TransmitEnergy = structures.Data.([acquisitionMode{1, i} 'TransmitEnergy']); + data.(acquisitionMode{1, i}).NominalCorrelation = structures.Data.([acquisitionMode{1, i} 'NominalCorrelation']); + + % support velocity data in ENU or beam coordinates + data.(acquisitionMode{1, i}).isVelocityData = false; + if isfield(structures.Data, [acquisitionMode{1, i} 'Velocity_ENU']) + data.(acquisitionMode{1, i}).Velocity_E = squeeze(structures.Data.([acquisitionMode{1, i} 'ENU'])(:, 1, :)); + data.(acquisitionMode{1, i}).Velocity_N = squeeze(structures.Data.([acquisitionMode{1, i} 'ENU'])(:, 2, :)); + data.(acquisitionMode{1, i}).Velocity_U = squeeze(structures.Data.([acquisitionMode{1, i} 'ENU'])(:, 3, :)); + data.(acquisitionMode{1, i}).Velocity_U2 = squeeze(structures.Data.([acquisitionMode{1, i} 'ENU'])(:, 4, :)); + + data.(acquisitionMode{1, i}).isVelocityData = true; + data.(acquisitionMode{1, i}).coordSys = 'ENU'; + elseif isfield(structures.Data, [acquisitionMode{1, i} 'VelEast']) + data.(acquisitionMode{1, i}).Velocity_E = structures.Data.([acquisitionMode{1, i} 'VelEast']); + data.(acquisitionMode{1, i}).Velocity_N = structures.Data.([acquisitionMode{1, i} 'VelNorth']); + if isfield(structures.Data, [acquisitionMode{1, i} 'VelUp']) + data.(acquisitionMode{1, i}).Velocity_U = structures.Data.([acquisitionMode{1, i} 'VelUp']); + else + data.(acquisitionMode{1, i}).Velocity_U = structures.Data.([acquisitionMode{1, i} 'VelUp1']); + end + data.(acquisitionMode{1, i}).Velocity_U2 = structures.Data.([acquisitionMode{1, i} 'VelUp2']); + + data.(acquisitionMode{1, i}).isVelocityData = true; + data.(acquisitionMode{1, i}).coordSys = 'ENU'; + elseif isfield(structures.Data, [acquisitionMode{1, i} 'VelBeam1']) + data.(acquisitionMode{1, i}).Velocity_1 = structures.Data.([acquisitionMode{1, i} 'VelBeam1']); + data.(acquisitionMode{1, i}).Velocity_2 = structures.Data.([acquisitionMode{1, i} 'VelBeam2']); + data.(acquisitionMode{1, i}).Velocity_3 = structures.Data.([acquisitionMode{1, i} 'VelBeam3']); + data.(acquisitionMode{1, i}).Velocity_4 = structures.Data.([acquisitionMode{1, i} 'VelBeam4']); + + data.(acquisitionMode{1, i}).isVelocityData = true; + data.(acquisitionMode{1, i}).coordSys = 'BEAM'; + end + if data.(acquisitionMode{1, i}).isVelocityData + data.(acquisitionMode{1, i}).Beam2xyz = structures.Config.([acquisitionMode{1, i} 'Beam2xyz']); + if isfield(structures.Data, [acquisitionMode{1, i} 'Amplitude_Beam']) + data.(acquisitionMode{1, i}).Backscatter1 = squeeze(structures.Data.([acquisitionMode{1, i} 'Amplitude_Beam'])(:, 1, :))*2; % looks like the .mat format is giving dB by default with dB = raw counts * 0.5 + data.(acquisitionMode{1, i}).Backscatter2 = squeeze(structures.Data.([acquisitionMode{1, i} 'Amplitude_Beam'])(:, 2, :))*2; + data.(acquisitionMode{1, i}).Backscatter3 = squeeze(structures.Data.([acquisitionMode{1, i} 'Amplitude_Beam'])(:, 3, :))*2; + data.(acquisitionMode{1, i}).Backscatter4 = squeeze(structures.Data.([acquisitionMode{1, i} 'Amplitude_Beam'])(:, 4, :))*2; + elseif isfield(structures.Data, [acquisitionMode{1, i} 'AmpBeam1']) + data.(acquisitionMode{1, i}).Backscatter1 = structures.Data.([acquisitionMode{1, i} 'AmpBeam1'])*2; + data.(acquisitionMode{1, i}).Backscatter2 = structures.Data.([acquisitionMode{1, i} 'AmpBeam2'])*2; + data.(acquisitionMode{1, i}).Backscatter3 = structures.Data.([acquisitionMode{1, i} 'AmpBeam3'])*2; + data.(acquisitionMode{1, i}).Backscatter4 = structures.Data.([acquisitionMode{1, i} 'AmpBeam4'])*2; + end + if isfield(structures.Data, [acquisitionMode{1, i} 'Correlation_Beam']) + data.(acquisitionMode{1, i}).Correlation1 = squeeze(structures.Data.([acquisitionMode{1, i} 'Correlation_Beam'])(:, 1, :)); + data.(acquisitionMode{1, i}).Correlation2 = squeeze(structures.Data.([acquisitionMode{1, i} 'Correlation_Beam'])(:, 2, :)); + data.(acquisitionMode{1, i}).Correlation3 = squeeze(structures.Data.([acquisitionMode{1, i} 'Correlation_Beam'])(:, 3, :)); + data.(acquisitionMode{1, i}).Correlation4 = squeeze(structures.Data.([acquisitionMode{1, i} 'Correlation_Beam'])(:, 4, :)); + elseif isfield(structures.Data, [acquisitionMode{1, i} 'CorBeam1']) + data.(acquisitionMode{1, i}).Correlation1 = structures.Data.([acquisitionMode{1, i} 'CorBeam1']); + data.(acquisitionMode{1, i}).Correlation2 = structures.Data.([acquisitionMode{1, i} 'CorBeam2']); + data.(acquisitionMode{1, i}).Correlation3 = structures.Data.([acquisitionMode{1, i} 'CorBeam3']); + data.(acquisitionMode{1, i}).Correlation4 = structures.Data.([acquisitionMode{1, i} 'CorBeam4']); + end + end + + data.(acquisitionMode{1, i}).isAltimeterData = false; + if isfield(structures.Data, [acquisitionMode{1, i} 'AltimeterDistanceLE']) + data.(acquisitionMode{1, i}).AltimeterDistanceLE = structures.Data.([acquisitionMode{1, i} 'AltimeterDistanceLE']); + data.(acquisitionMode{1, i}).AltimeterQualityLE = structures.Data.([acquisitionMode{1, i} 'AltimeterQualityLE']); +% data.(acquisitionMode{1, i}).AltimeterStatus = structures.Data.([acquisitionMode{1, i} 'AltimeterStatus']); % flags for when pitch or roll is > 5 or > 10. We already have pitch and roll data... + + data.(acquisitionMode{1, i}).isAltimeterData = true; + end + + data.(acquisitionMode{1, i}).isASTData = false; + if isfield(structures.Data, [acquisitionMode{1, i} 'AltimeterDistanceAST']) + data.(acquisitionMode{1, i}).AltimeterDistanceAST = structures.Data.([acquisitionMode{1, i} 'AltimeterDistanceAST']); + data.(acquisitionMode{1, i}).AltimeterQualityAST = structures.Data.([acquisitionMode{1, i} 'AltimeterQualityAST']); + data.(acquisitionMode{1, i}).AltimeterTimeOffsetAST = structures.Data.([acquisitionMode{1, i} 'AltimeterTimeOffsetAST']); + data.(acquisitionMode{1, i}).AltimeterPressure = structures.Data.([acquisitionMode{1, i} 'AltimeterPressure']); + + data.(acquisitionMode{1, i}).isASTData = true; + end + + data.(acquisitionMode{1, i}).Battery = structures.Data.([acquisitionMode{1, i} 'Battery']); + data.(acquisitionMode{1, i}).Heading = structures.Data.([acquisitionMode{1, i} 'Heading']); + data.(acquisitionMode{1, i}).Pitch = structures.Data.([acquisitionMode{1, i} 'Pitch']); + data.(acquisitionMode{1, i}).Roll = structures.Data.([acquisitionMode{1, i} 'Roll']); + data.(acquisitionMode{1, i}).Temperature = structures.Data.([acquisitionMode{1, i} 'Temperature']); + data.(acquisitionMode{1, i}).speedOfSound = structures.Data.([acquisitionMode{1, i} 'Soundspeed']); + data.(acquisitionMode{1, i}).Pressure = structures.Data.([acquisitionMode{1, i} 'Pressure']); + + if ~isempty(acquisitionMode{2, i}) + % we're adding the similar dataset to the existing one + data.(acquisitionMode{1, i}).Time = [data.(acquisitionMode{1, i}).Time; ... + structures.Data.([acquisitionMode{2, i} 'Time'])]; + + data.(acquisitionMode{1, i}).nSamples = length(data.(acquisitionMode{1, i}).Time); + + data.(acquisitionMode{1, i}).Status = [data.(acquisitionMode{1, i}).Status; ... + dec2bin(bytecast(structures.Data.([acquisitionMode{2, i} 'Status']), 'L', 'uint32', cpuEndianness), 32)]; + data.(acquisitionMode{1, i}).Error = [data.(acquisitionMode{1, i}).Error; ... + structures.Data.([acquisitionMode{2, i} 'Error'])]; % error codes for each cell of a velocity profile inferred from the beams. 0=good; otherwise error. See http://www.nortek-as.com/en/knowledge-center/forum/waves/20001875?b_start=0#769595815 + data.(acquisitionMode{1, i}).AmbiguityVel = [data.(acquisitionMode{1, i}).AmbiguityVel; ... + structures.Data.([acquisitionMode{2, i} 'AmbiguityVel'])]; + data.(acquisitionMode{1, i}).TransmitEnergy = [data.(acquisitionMode{1, i}).TransmitEnergy; ... + structures.Data.([acquisitionMode{2, i} 'TransmitEnergy'])]; + data.(acquisitionMode{1, i}).NominalCorrelation = [data.(acquisitionMode{1, i}).NominalCorrelation; ... + structures.Data.([acquisitionMode{2, i} 'NominalCorrelation'])]; + + % only support velocity data in ENU coordinates + if isfield(structures.Data, [acquisitionMode{2, i} 'Velocity_ENU']) + data.(acquisitionMode{1, i}).Velocity_E = [data.(acquisitionMode{1, i}).Velocity_E; ... + squeeze(structures.Data.([acquisitionMode{2, i} 'ENU'])(:, 1, :))]; + data.(acquisitionMode{1, i}).Velocity_N = [data.(acquisitionMode{1, i}).Velocity_N; ... + squeeze(structures.Data.([acquisitionMode{2, i} 'ENU'])(:, 2, :))]; + data.(acquisitionMode{1, i}).Velocity_U = [data.(acquisitionMode{1, i}).Velocity_U; ... + squeeze(structures.Data.([acquisitionMode{2, i} 'ENU'])(:, 3, :))]; + data.(acquisitionMode{1, i}).Velocity_U2 = [data.(acquisitionMode{1, i}).Velocity_U2; ... + squeeze(structures.Data.([acquisitionMode{2, i} 'ENU'])(:, 4, :))]; + + elseif isfield(structures.Data, [acquisitionMode{2, i} 'VelEast']) + data.(acquisitionMode{1, i}).Velocity_E = [data.(acquisitionMode{1, i}).Velocity_E; ... + structures.Data.([acquisitionMode{2, i} 'VelEast'])]; + data.(acquisitionMode{1, i}).Velocity_N = [data.(acquisitionMode{1, i}).Velocity_N; ... + structures.Data.([acquisitionMode{2, i} 'VelNorth'])]; + if isfield(structures.Data, [acquisitionMode{2, i} 'VelUp']) + data.(acquisitionMode{1, i}).Velocity_U = [data.(acquisitionMode{1, i}).Velocity_U; ... + structures.Data.([acquisitionMode{2, i} 'VelUp'])]; + else + data.(acquisitionMode{1, i}).Velocity_U = [data.(acquisitionMode{1, i}).Velocity_U; ... + structures.Data.([acquisitionMode{2, i} 'VelUp1'])]; + end + data.(acquisitionMode{1, i}).Velocity_U2 = [data.(acquisitionMode{1, i}).Velocity_U2; ... + structures.Data.([acquisitionMode{2, i} 'VelUp2'])]; + + elseif isfield(structures.Data, [acquisitionMode{1, i} 'VelBeam1']) + data.(acquisitionMode{1, i}).Velocity_1 = [data.(acquisitionMode{1, i}).Velocity_1; ... + structures.Data.([acquisitionMode{2, i} 'VelBeam1'])]; + data.(acquisitionMode{1, i}).Velocity_2 = [data.(acquisitionMode{1, i}).Velocity_2; ... + structures.Data.([acquisitionMode{2, i} 'VelBeam2'])]; + data.(acquisitionMode{1, i}).Velocity_3 = [data.(acquisitionMode{1, i}).Velocity_3; ... + structures.Data.([acquisitionMode{2, i} 'VelBeam3'])]; + data.(acquisitionMode{1, i}).Velocity_4 = [data.(acquisitionMode{1, i}).Velocity_4; ... + structures.Data.([acquisitionMode{2, i} 'VelBeam4'])]; + + end + if data.(acquisitionMode{1, i}).isVelocityData + if isfield(structures.Data, [acquisitionMode{2, i} 'Amplitude_Beam']) + data.(acquisitionMode{1, i}).Backscatter1 = [data.(acquisitionMode{1, i}).Backscatter1; ... + squeeze(structures.Data.([acquisitionMode{2, i} 'Amplitude_Beam'])(:, 1, :))*2]; % looks like the .mat format is giving dB by default with dB = raw counts * 0.5 + data.(acquisitionMode{1, i}).Backscatter2 = [data.(acquisitionMode{1, i}).Backscatter2; ... + squeeze(structures.Data.([acquisitionMode{2, i} 'Amplitude_Beam'])(:, 2, :))*2]; + data.(acquisitionMode{1, i}).Backscatter3 = [data.(acquisitionMode{1, i}).Backscatter3; ... + squeeze(structures.Data.([acquisitionMode{2, i} 'Amplitude_Beam'])(:, 3, :))*2]; + data.(acquisitionMode{1, i}).Backscatter4 = [data.(acquisitionMode{1, i}).Backscatter4; ... + squeeze(structures.Data.([acquisitionMode{2, i} 'Amplitude_Beam'])(:, 4, :))*2]; + elseif isfield(structures.Data, [acquisitionMode{2, i} 'AmpBeam1']) + data.(acquisitionMode{1, i}).Backscatter1 = [data.(acquisitionMode{1, i}).Backscatter1; ... + structures.Data.([acquisitionMode{2, i} 'AmpBeam1'])*2]; + data.(acquisitionMode{1, i}).Backscatter2 = [data.(acquisitionMode{1, i}).Backscatter2; ... + structures.Data.([acquisitionMode{2, i} 'AmpBeam2'])*2]; + data.(acquisitionMode{1, i}).Backscatter3 = [data.(acquisitionMode{1, i}).Backscatter3; ... + structures.Data.([acquisitionMode{2, i} 'AmpBeam3'])*2]; + data.(acquisitionMode{1, i}).Backscatter4 = [data.(acquisitionMode{1, i}).Backscatter4; ... + structures.Data.([acquisitionMode{2, i} 'AmpBeam4'])*2]; + end + if isfield(structures.Data, [acquisitionMode{2, i} 'Correlation_Beam']) + data.(acquisitionMode{1, i}).Correlation1 = [data.(acquisitionMode{1, i}).Correlation1; ... + squeeze(structures.Data.([acquisitionMode{2, i} 'Correlation_Beam'])(:, 1, :))]; + data.(acquisitionMode{1, i}).Correlation2 = [data.(acquisitionMode{1, i}).Correlation2; ... + squeeze(structures.Data.([acquisitionMode{2, i} 'Correlation_Beam'])(:, 2, :))]; + data.(acquisitionMode{1, i}).Correlation3 = [data.(acquisitionMode{1, i}).Correlation3; ... + squeeze(structures.Data.([acquisitionMode{2, i} 'Correlation_Beam'])(:, 3, :))]; + data.(acquisitionMode{1, i}).Correlation4 = [data.(acquisitionMode{1, i}).Correlation4; ... + squeeze(structures.Data.([acquisitionMode{2, i} 'Correlation_Beam'])(:, 4, :))]; + elseif isfield(structures.Data, [acquisitionMode{2, i} 'CorBeam1']) + data.(acquisitionMode{1, i}).Correlation1 = [data.(acquisitionMode{1, i}).Correlation1; ... + structures.Data.([acquisitionMode{2, i} 'CorBeam1'])]; + data.(acquisitionMode{1, i}).Correlation2 = [data.(acquisitionMode{1, i}).Correlation2; ... + structures.Data.([acquisitionMode{2, i} 'CorBeam2'])]; + data.(acquisitionMode{1, i}).Correlation3 = [data.(acquisitionMode{1, i}).Correlation3; ... + structures.Data.([acquisitionMode{2, i} 'CorBeam3'])]; + data.(acquisitionMode{1, i}).Correlation4 = [data.(acquisitionMode{1, i}).Correlation4; ... + structures.Data.([acquisitionMode{2, i} 'CorBeam4'])]; + end + end + + if isfield(structures.Data, [acquisitionMode{2, i} 'AltimeterDistanceLE']) + data.(acquisitionMode{1, i}).AltimeterDistanceLE = [data.(acquisitionMode{1, i}).AltimeterDistanceLE; ... + structures.Data.([acquisitionMode{2, i} 'AltimeterDistanceLE'])]; + data.(acquisitionMode{1, i}).AltimeterQualityLE = [data.(acquisitionMode{1, i}).AltimeterQualityLE; ... + structures.Data.([acquisitionMode{2, i} 'AltimeterQualityLE'])]; +% data.(acquisitionMode{i}).AltimeterStatus = [data.(acquisitionMode{1, i}).AltimeterStatus; ... +% structures.Data.([acquisitionMode{2, i} 'AltimeterStatus'])]; % flags for when pitch or roll is > 5 or > 10. We already have pitch and roll data... + end + + if isfield(structures.Data, [acquisitionMode{2, i} 'AltimeterDistanceAST']) + data.(acquisitionMode{1, i}).AltimeterDistanceAST = [data.(acquisitionMode{1, i}).AltimeterDistanceAST; ... + structures.Data.([acquisitionMode{2, i} 'AltimeterDistanceAST'])]; + data.(acquisitionMode{1, i}).AltimeterQualityAST = [data.(acquisitionMode{1, i}).AltimeterQualityAST; ... + structures.Data.([acquisitionMode{2, i} 'AltimeterQualityAST'])]; + data.(acquisitionMode{1, i}).AltimeterTimeOffsetAST = [data.(acquisitionMode{1, i}).AltimeterTimeOffsetAST; ... + structures.Data.([acquisitionMode{2, i} 'AltimeterTimeOffsetAST'])]; + data.(acquisitionMode{1, i}).AltimeterPressure = [data.(acquisitionMode{1, i}).AltimeterPressure; ... + structures.Data.([acquisitionMode{2, i} 'AltimeterPressure'])]; + end + + data.(acquisitionMode{1, i}).Battery = [data.(acquisitionMode{1, i}).Battery; ... + structures.Data.([acquisitionMode{2, i} 'Battery'])]; + data.(acquisitionMode{1, i}).Heading = [data.(acquisitionMode{1, i}).Heading; ... + structures.Data.([acquisitionMode{2, i} 'Heading'])]; + data.(acquisitionMode{1, i}).Pitch = [data.(acquisitionMode{1, i}).Pitch; ... + structures.Data.([acquisitionMode{2, i} 'Pitch'])]; + data.(acquisitionMode{1, i}).Roll = [data.(acquisitionMode{1, i}).Roll; ... + structures.Data.([acquisitionMode{2, i} 'Roll'])]; + data.(acquisitionMode{1, i}).Temperature = [data.(acquisitionMode{1, i}).Temperature; ... + structures.Data.([acquisitionMode{2, i} 'Temperature'])]; + data.(acquisitionMode{1, i}).speedOfSound = [data.(acquisitionMode{1, i}).speedOfSound; ... + structures.Data.([acquisitionMode{2, i} 'Soundspeed'])]; + data.(acquisitionMode{1, i}).Pressure = [data.(acquisitionMode{1, i}).Pressure; ... + structures.Data.([acquisitionMode{2, i} 'Pressure'])]; + + % we now need to sort the data chronologically + [~, iSort] = sort(data.(acquisitionMode{1, i}).Time); + nSample = length(data.(acquisitionMode{1, i}).Time); + dataFields = fieldnames(data.(acquisitionMode{1, i})); + for j=1:length(dataFields) + if size(data.(acquisitionMode{1, i}).(dataFields{j}), 1) == nSample + data.(acquisitionMode{1, i}).(dataFields{j}) = data.(acquisitionMode{1, i}).(dataFields{j})(iSort, :); + end + end + + end end - Battery = structures.Data.([acquisitionMode '_Battery']); - Heading = structures.Data.([acquisitionMode '_Heading']); - Pitch = structures.Data.([acquisitionMode '_Pitch']); - Roll = structures.Data.([acquisitionMode '_Roll']); - Temperature = structures.Data.([acquisitionMode '_Temperature']); - speedOfSound = structures.Data.([acquisitionMode '_Soundspeed']); - Pressure = structures.Data.([acquisitionMode '_Pressure']); clear structures; + + acquisitionMode = acquisitionMode(1, :); otherwise error('Data format not supported'); end -sample_data = struct; - if magDec ~= 0 isMagBias = true; magBiasComment = ['A compass correction of ' num2str(magDec) ... @@ -224,131 +592,213 @@ '(usually to account for magnetic declination).']; end -sample_data.toolbox_input_file = filename; -sample_data.meta.featureType = ''; % strictly this dataset cannot be described as timeSeriesProfile since it also includes timeSeries data like TEMP -sample_data.meta.binSize = cellSize; -sample_data.meta.instrument_make = 'Nortek'; -if isempty(instrumentModel) - sample_data.meta.instrument_model = 'Signature'; -else - sample_data.meta.instrument_model = instrumentModel; +sample_data = cell(1, nDataset); +for i=1:nDataset + sample_data{i}.toolbox_input_file = filename; + sample_data{i}.meta.featureType = ''; % strictly this dataset cannot be described as timeSeriesProfile since it also includes timeSeries data like TEMP + sample_data{i}.meta.binSize = data.(acquisitionMode{i}).cellSize; + sample_data{i}.meta.nBeams = data.(acquisitionMode{i}).nBeams; + sample_data{i}.meta.instrument_make = 'Nortek'; + if isempty(instrumentModel) + sample_data{i}.meta.instrument_model = 'Signature'; + else + sample_data{i}.meta.instrument_model = instrumentModel; + end + sample_data{i}.meta.instrument_serial_no = serialNumber; -end -sample_data.meta.instrument_serial_no = serialNumber; -sample_data.meta.instrument_sample_interval = median(diff(Time*24*3600)); -switch sample_data.meta.instrument_model - case 'Signature250' - sample_data.meta.beam_angle = 20; - otherwise % Signature500 and Signature1000 - sample_data.meta.beam_angle = 25; -end - -% generate distance values -distance = nan(nCells, 1); -distance(:) = (blankDist): ... - (cellSize): ... - (blankDist + (nCells-1) * cellSize); - -% distance between the ADCP's transducers and the middle of each cell -% See http://www.nortek-bv.nl/en/knowledge-center/forum/current-profilers-and-current-meters/579860330 -distance = (distance + cellSize); - -% add dimensions with their data mapped -iStartOrientation = 26; -iEndOrientation = 28; -adcpOrientations = bin2dec(Status(:, end-iEndOrientation+1:end-iStartOrientation+1)); -adcpOrientation = mode(adcpOrientations); % hopefully the most frequent value reflects the orientation when deployed -% we assume adcpOrientation == 4 by default "ZUP" -if adcpOrientation == 5 - % case of a downward looking ADCP -> negative values - distance = -distance; -end -iBadOriented = any(adcpOrientations ~= repmat(adcpOrientation, nSamples, 1), 2); % we'll only keep velocity data collected when ADCP is oriented as expected -Velocity_N(iBadOriented, :) = NaN; -Velocity_E(iBadOriented, :) = NaN; -Velocity_U(iBadOriented, :) = NaN; -Velocity_U2(iBadOriented, :) = NaN; -Backscatter1(iBadOriented, :) = NaN; -Backscatter2(iBadOriented, :) = NaN; -Backscatter3(iBadOriented, :) = NaN; -Backscatter4(iBadOriented, :) = NaN; -Correlation1(iBadOriented, :) = NaN; -Correlation2(iBadOriented, :) = NaN; -Correlation3(iBadOriented, :) = NaN; -Correlation4(iBadOriented, :) = NaN; -dims = { - 'TIME', Time, ''; ... - 'DIST_ALONG_BEAMS', distance, 'Nortek instrument data is not vertically bin-mapped (no tilt correction applied). Cells are lying parallel to the beams, at heights above sensor that vary with tilt.' - }; -clear Time distance; - -nDims = size(dims, 1); -sample_data.dimensions = cell(nDims, 1); -for i=1:nDims - sample_data.dimensions{i}.name = dims{i, 1}; - sample_data.dimensions{i}.typeCastFunc = str2func(netcdf3ToMatlabType(imosParameters(dims{i, 1}, 'type'))); - sample_data.dimensions{i}.data = sample_data.dimensions{i}.typeCastFunc(dims{i, 2}); - sample_data.dimensions{i}.comment = dims{i, 3}; -end -clear dims; - -% add variables with their dimensions and data mapped. -% we assume no correction for magnetic declination has been applied -if isMagBias - magExt = ''; -else - magExt = '_MAG'; -end -iDimVel = nDims; -iDimDiag = nDims; -vars = { - 'TIMESERIES', [], 1;... - 'LATITUDE', [], NaN; ... - 'LONGITUDE', [], NaN; ... - 'NOMINAL_DEPTH', [], NaN; ... - ['VCUR' magExt], [1 iDimVel], Velocity_N; ... % V - ['UCUR' magExt], [1 iDimVel], Velocity_E; ... % U - 'WCUR', [1 iDimVel], Velocity_U; ... - 'WCUR_2', [1 iDimVel], Velocity_U2; ... - 'ABSIC1', [1 iDimDiag], Backscatter1; ... - 'ABSIC2', [1 iDimDiag], Backscatter2; ... - 'ABSIC3', [1 iDimDiag], Backscatter3; ... - 'ABSIC4', [1 iDimDiag], Backscatter4; ... - 'CMAG1', [1 iDimDiag], Correlation1; ... - 'CMAG2', [1 iDimDiag], Correlation2; ... - 'CMAG3', [1 iDimDiag], Correlation3; ... - 'CMAG4', [1 iDimDiag], Correlation4; ... - 'TEMP', 1, Temperature; ... - 'PRES_REL', 1, Pressure; ... - 'SSPD', 1, speedOfSound; ... - 'VOLT', 1, Battery; ... - 'PITCH', 1, Pitch; ... - 'ROLL', 1, Roll; ... - ['HEADING' magExt], 1, Heading - }; -clear Velocity_N Velocity_E Velocity_U Velocity_U2 ... - Backscatter1 Backscatter2 Backscatter3 Backscatter4 ... - Correlation1 Correlation2 Correlation3 Correlation4 ... - Temperature Pressure speedOfSound Battery Pitch Roll Heading Status; - -nVars = size(vars, 1); -sample_data.variables = cell(nVars, 1); -for i=1:nVars - sample_data.variables{i}.name = vars{i, 1}; - sample_data.variables{i}.typeCastFunc = str2func(netcdf3ToMatlabType(imosParameters(vars{i, 1}, 'type'))); - sample_data.variables{i}.dimensions = vars{i, 2}; - if ~isempty(vars{i, 2}) % we don't want this for scalar variables - if length(sample_data.variables{i}.dimensions) == 2 - sample_data.variables{i}.coordinates = ['TIME LATITUDE LONGITUDE ' sample_data.dimensions{sample_data.variables{i}.dimensions(2)}.name]; - else - sample_data.variables{i}.coordinates = 'TIME LATITUDE LONGITUDE NOMINAL_DEPTH'; + diffTimeInSec = diff(data.(acquisitionMode{i}).Time*24*3600); + + % look for most frequents modes + Ms = unique(round(diffTimeInSec*100)/100); % unique() sorts in order of most frequent + + sample_data{i}.meta.instrument_sample_interval = Ms(1); + + % we look for the second burst (less likely to be chopped) + % and we hope it is valid + dt = [0; diffTimeInSec]; + iBurst = find(dt >= Ms(2) - 1/100, 2, 'first'); + if length(iBurst) == 2 + sample_data{i}.meta.instrument_burst_interval = round((data.(acquisitionMode{i}).Time(iBurst(2)) - ... + data.(acquisitionMode{i}).Time(iBurst(1)))*24*3600 * 100)/100; + + sample_data{i}.meta.instrument_burst_duration = round((data.(acquisitionMode{i}).Time(iBurst(2)-1) - ... + data.(acquisitionMode{i}).Time(iBurst(1)))*24*3600 * 100)/100; + end + + switch sample_data{i}.meta.instrument_model + case 'Signature250' + sample_data{i}.meta.beam_angle = 20; + otherwise % Signature500 and Signature1000 + sample_data{i}.meta.beam_angle = 25; + end +% if data.(acquisitionMode{i}).isVelocityData +% sample_data{i}.meta.beam_to_xyz_transform = data.(acquisitionMode{i}).Beam2xyz; % need to find out how to compute this matrix for .ad2cp format. See https://pypkg.com/pypi/nortek/f/nortek/files.py +% end + + % generate distance values + distance = nan(data.(acquisitionMode{i}).nCells, 1); + distance(:) = (data.(acquisitionMode{i}).blankDist): ... + (data.(acquisitionMode{i}).cellSize): ... + (data.(acquisitionMode{i}).blankDist + (data.(acquisitionMode{i}).nCells-1) * data.(acquisitionMode{i}).cellSize); + + % distance between the ADCP's transducers and the middle of each cell + % See http://www.nortek-bv.nl/en/knowledge-center/forum/current-profilers-and-current-meters/579860330 + distance = (distance + data.(acquisitionMode{i}).cellSize); + + % add dimensions with their data mapped + iStartOrientation = 26; + iEndOrientation = 28; + adcpOrientations = bin2dec(data.(acquisitionMode{i}).Status(:, end-iEndOrientation+1:end-iStartOrientation+1)); + adcpOrientation = mode(adcpOrientations); % hopefully the most frequent value reflects the orientation when deployed + % we assume adcpOrientation == 4 by default "ZUP" + if adcpOrientation == 5 + % case of a downward looking ADCP -> negative values + distance = -distance; + end + + dims = {'TIME', data.(acquisitionMode{i}).Time, ''}; + + if data.(acquisitionMode{i}).isVelocityData + % we'll only keep velocity data collected when ADCP is oriented as expected + iBadOriented = any(adcpOrientations ~= repmat(adcpOrientation, data.(acquisitionMode{i}).nSamples, 1), 2); + + switch data.(acquisitionMode{i}).coordSys + case 'ENU' + data.(acquisitionMode{i}).Velocity_N (iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Velocity_E (iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Velocity_U (iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Velocity_U2 (iBadOriented, :) = NaN; + + case 'BEAM' + data.(acquisitionMode{i}).Velocity_1 (iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Velocity_2 (iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Velocity_3 (iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Velocity_4 (iBadOriented, :) = NaN; + end + data.(acquisitionMode{i}).Backscatter1(iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Backscatter2(iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Backscatter3(iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Backscatter4(iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Correlation1(iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Correlation2(iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Correlation3(iBadOriented, :) = NaN; + data.(acquisitionMode{i}).Correlation4(iBadOriented, :) = NaN; + + dims = [dims; ... + {'DIST_ALONG_BEAMS', distance, 'Nortek instrument data is not vertically bin-mapped (no tilt correction applied). Cells are lying parallel to the beams, at heights above sensor that vary with tilt.'}]; end - sample_data.variables{i}.data = sample_data.variables{i}.typeCastFunc(vars{i, 3}); - if any(strcmpi(vars{i, 1}, {'VCUR', 'UCUR', 'HEADING'})) - sample_data.variables{i}.compass_correction_applied = magDec; - sample_data.variables{i}.comment = magBiasComment; + nDims = size(dims, 1); + sample_data{i}.dimensions = cell(nDims, 1); + for j=1:nDims + sample_data{i}.dimensions{j}.name = dims{j, 1}; + sample_data{i}.dimensions{j}.typeCastFunc = str2func(netcdf3ToMatlabType(imosParameters(dims{j, 1}, 'type'))); + sample_data{i}.dimensions{j}.data = sample_data{i}.dimensions{j}.typeCastFunc(dims{j, 2}); + sample_data{i}.dimensions{j}.comment = dims{j, 3}; end -end -clear vars; + clear dims; + + % add variables with their dimensions and data mapped. + % we assume no correction for magnetic declination has been applied + if isMagBias + magExt = ''; + else + magExt = '_MAG'; + end + iDimVel = nDims; + iDimDiag = nDims; + vars = { + 'TIMESERIES', [], 1;... + 'LATITUDE', [], NaN; ... + 'LONGITUDE', [], NaN; ... + 'NOMINAL_DEPTH', [], NaN + }; + + if data.(acquisitionMode{i}).isVelocityData + switch data.(acquisitionMode{i}).coordSys + case 'ENU' + vars = [vars; { + ['VCUR' magExt], [1 iDimVel], data.(acquisitionMode{i}).Velocity_N; ... % V + ['UCUR' magExt], [1 iDimVel], data.(acquisitionMode{i}).Velocity_E; ... % U + 'WCUR', [1 iDimVel], data.(acquisitionMode{i}).Velocity_U; ... + 'WCUR_2', [1 iDimVel], data.(acquisitionMode{i}).Velocity_U2 + }]; + + case 'BEAM' + vars = [vars; { + 'VEL1', [1 iDimVel], data.(acquisitionMode{i}).Velocity_1; ... + 'VEL2', [1 iDimVel], data.(acquisitionMode{i}).Velocity_2; ... + 'VEL3', [1 iDimVel], data.(acquisitionMode{i}).Velocity_3; ... + 'VEL4', [1 iDimVel], data.(acquisitionMode{i}).Velocity_4 + }]; + + end + vars = [vars; { + 'ABSIC1', [1 iDimDiag], data.(acquisitionMode{i}).Backscatter1; ... + 'ABSIC2', [1 iDimDiag], data.(acquisitionMode{i}).Backscatter2; ... + 'ABSIC3', [1 iDimDiag], data.(acquisitionMode{i}).Backscatter3; ... + 'ABSIC4', [1 iDimDiag], data.(acquisitionMode{i}).Backscatter4; ... + 'ABSI1', [1 iDimDiag], data.(acquisitionMode{i}).Backscatter1*0.5; ... % 1 count = 0.5 dB according to manual + 'ABSI2', [1 iDimDiag], data.(acquisitionMode{i}).Backscatter2*0.5; ... + 'ABSI3', [1 iDimDiag], data.(acquisitionMode{i}).Backscatter3*0.5; ... + 'ABSI4', [1 iDimDiag], data.(acquisitionMode{i}).Backscatter4*0.5; ... + 'CMAG1', [1 iDimDiag], data.(acquisitionMode{i}).Correlation1; ... + 'CMAG2', [1 iDimDiag], data.(acquisitionMode{i}).Correlation2; ... + 'CMAG3', [1 iDimDiag], data.(acquisitionMode{i}).Correlation3; ... + 'CMAG4', [1 iDimDiag], data.(acquisitionMode{i}).Correlation4 + }]; + end + + if data.(acquisitionMode{i}).isAltimeterData + vars = [vars; { + 'LE_DIST', 1, data.(acquisitionMode{i}).AltimeterDistanceLE; ... + 'LE_QUALITY', 1, data.(acquisitionMode{i}).AltimeterQualityLE + }]; + end + + if data.(acquisitionMode{i}).isASTData + vars = [vars; { + 'AST_DIST', 1, data.(acquisitionMode{i}).AltimeterDistanceAST; ... + 'AST_QUALITY', 1, data.(acquisitionMode{i}).AltimeterQualityAST; ... + 'AST_TIME_OFFSET', 1, data.(acquisitionMode{i}).AltimeterTimeOffsetAST; ... + 'ALTIMETER_PRES', 1, data.(acquisitionMode{i}).AltimeterPressure + }]; + end + + vars = [vars; { + 'TEMP', 1, data.(acquisitionMode{i}).Temperature; ... + 'PRES_REL', 1, data.(acquisitionMode{i}).Pressure; ... + 'SSPD', 1, data.(acquisitionMode{i}).speedOfSound; ... + 'VOLT', 1, data.(acquisitionMode{i}).Battery; ... + 'PITCH', 1, data.(acquisitionMode{i}).Pitch; ... + 'ROLL', 1, data.(acquisitionMode{i}).Roll; ... + ['HEADING' magExt], 1, data.(acquisitionMode{i}).Heading; ... + 'ERROR', 1, data.(acquisitionMode{i}).Error; ... + 'AMBIG_VEL', 1, data.(acquisitionMode{i}).AmbiguityVel; ... + 'TRANSMIT_E', 1, data.(acquisitionMode{i}).TransmitEnergy; ... + 'NOMINAL_CORR', 1, data.(acquisitionMode{i}).NominalCorrelation + }]; + + nVars = size(vars, 1); + sample_data{i}.variables = cell(nVars, 1); + for k=1:nVars + sample_data{i}.variables{k}.name = vars{k, 1}; + sample_data{i}.variables{k}.typeCastFunc = str2func(netcdf3ToMatlabType(imosParameters(vars{k, 1}, 'type'))); + sample_data{i}.variables{k}.dimensions = vars{k, 2}; + if ~isempty(vars{k, 2}) % we don't want this for scalar variables + if length(sample_data{i}.variables{k}.dimensions) == 2 + sample_data{i}.variables{k}.coordinates = ['TIME LATITUDE LONGITUDE ' sample_data{i}.dimensions{sample_data{i}.variables{k}.dimensions(2)}.name]; + else + sample_data{i}.variables{k}.coordinates = 'TIME LATITUDE LONGITUDE NOMINAL_DEPTH'; + end + end + sample_data{i}.variables{k}.data = sample_data{i}.variables{k}.typeCastFunc(vars{k, 3}); + + if any(strcmpi(vars{k, 1}, {'VCUR', 'UCUR', 'HEADING'})) + sample_data{i}.variables{k}.compass_correction_applied = magDec; + sample_data{i}.variables{k}.comment = magBiasComment; + end + end + clear vars; +end \ No newline at end of file diff --git a/Preprocessing/oxygenPP.m b/Preprocessing/oxygenPP.m index 839f1bb92..334c0fe06 100644 --- a/Preprocessing/oxygenPP.m +++ b/Preprocessing/oxygenPP.m @@ -86,7 +86,7 @@ tempIdx = getVar(sam.variables, 'TEMP'); psalIdx = getVar(sam.variables, 'PSAL'); - [presRel, presName] = getPresRelForGSW(sam); + [presRel, zName, zComment] = getPresRelForGSW(sam); doxIdx = getVar(sam.variables, 'DOX'); doxyIdx = getVar(sam.variables, 'DOXY'); @@ -95,7 +95,7 @@ doxsIdx = getVar(sam.variables, 'DOXS'); % cancel if psal, temp, pres/pres_rel or depth/nominal depth and any DO parameter not present in data set - if ~(psalIdx && tempIdx && ~isempty(presName) && (doxIdx || doxyIdx || dox1Idx || dox2Idx || doxsIdx)), continue; end + if ~(psalIdx && tempIdx && ~isempty(zName) && (doxIdx || doxyIdx || dox1Idx || dox2Idx || doxsIdx)), continue; end % data set already contains dissolved oxygen DOX1/DOX2 or DOXS if (dox1Idx && dox2Idx && doxsIdx), continue; end @@ -115,7 +115,7 @@ CT = gsw_CT_from_pt(SA, potTemp); potDens = gsw_rho(SA, CT, 0); potDensComment = ['Potential density at atmospheric ' ... - 'pressure was derived from TEMP, PSAL, ' presName ', LATITUDE ' ... + 'pressure was derived from TEMP, PSAL, ' zName ' ' zComment ', LATITUDE ' ... 'and LONGITUDE using gsw_rho, gsw_SA_from_SP, gsw_CT_from_pt and gsw_pt0_from_t ' ... 'from the Gibbs-SeaWater toolbox (TEOS-10) v3.06.']; @@ -129,8 +129,8 @@ % calculate oxygen solubility at atmospheric pressure oxsolSurf = gsw_O2sol_SP_pt(psal, potTemp); % sea bird calculates OXSOL using psal and local temperature, not potential temperature - oxsolSurfComment = ['OXSOL_SURFACE derived from TEMP, PSAL, ' presName ... - ', LATITUDE and LONGITUDE using gsw_O2sol_SP_pt, gsw_pt0_from_t ' ... + oxsolSurfComment = ['OXSOL_SURFACE derived from TEMP, PSAL, ' zName ... + ' ' zComment ', LATITUDE and LONGITUDE using gsw_O2sol_SP_pt, gsw_pt0_from_t ' ... 'and gsw_SA_from_SP from the Gibbs-SeaWater toolbox (TEOS-10) v3.06.']; % add oxygen solubility at atmospheric pressure as new variable in data set diff --git a/Preprocessing/salinityPP.m b/Preprocessing/salinityPP.m index e179bca6f..eb9c4b6ef 100644 --- a/Preprocessing/salinityPP.m +++ b/Preprocessing/salinityPP.m @@ -61,10 +61,10 @@ cndcIdx = getVar(sam.variables, 'CNDC'); tempIdx = getVar(sam.variables, 'TEMP'); - [presRel, presName] = getPresRelForGSW(sam); + [presRel, zName, zComment] = getPresRelForGSW(sam); % cndc, temp, and pres/pres_rel or depth/nominal depth not present in data set - if ~(cndcIdx && tempIdx && ~isempty(presName)), continue; end + if ~(cndcIdx && tempIdx && ~isempty(zName)), continue; end cndc = sam.variables{cndcIdx}.data; temp = sam.variables{tempIdx}.data; @@ -77,7 +77,7 @@ psal = gsw_SP_from_R(R, temp, presRel); dimensions = sam.variables{tempIdx}.dimensions; - salinityComment = ['salinityPP.m: derived from CNDC, TEMP and ' presName ' using the Gibbs-SeaWater toolbox (TEOS-10) v3.06']; + salinityComment = ['salinityPP.m: derived from CNDC, TEMP and ' zName ' ' zComment ' using the Gibbs-SeaWater toolbox (TEOS-10) v3.06']; if isfield(sam.variables{tempIdx}, 'coordinates') coordinates = sam.variables{tempIdx}.coordinates; diff --git a/Seawater/TEOS10/gsw_CT_from_t.m b/Seawater/TEOS10/gsw_CT_from_t.m new file mode 100644 index 000000000..29d0a980d --- /dev/null +++ b/Seawater/TEOS10/gsw_CT_from_t.m @@ -0,0 +1,96 @@ +function CT = gsw_CT_from_t(SA,t,p) + +% gsw_CT_from_t Conservative Temperature from in-situ temperature +%========================================================================== +% +% USAGE: +% CT = gsw_CT_from_t(SA,t,p) +% +% DESCRIPTION: +% Calculates Conservative Temperature of seawater from in-situ +% temperature. +% +% INPUT: +% SA = Absolute Salinity [ g/kg ] +% t = in-situ temperature (ITS-90) [ deg C ] +% p = sea pressure [ dbar ] +% ( i.e. absolute pressure - 10.1325 dbar ) +% +% SA & t need to have the same dimensions. +% p may have dimensions 1x1 or Mx1 or 1xN or MxN, where SA & t are MxN. +% +% OUTPUT: +% CT = Conservative Temperature (ITS-90) [ deg C ] +% +% AUTHOR: +% David Jackett, Trevor McDougall and Paul Barker [ help@teos-10.org ] +% +% VERSION NUMBER: 3.05 (27th January 2015) +% +% REFERENCES: +% IOC, SCOR and IAPSO, 2010: The international thermodynamic equation of +% seawater - 2010: Calculation and use of thermodynamic properties. +% Intergovernmental Oceanographic Commission, Manuals and Guides No. 56, +% UNESCO (English), 196 pp. Available from http://www.TEOS-10.org +% See section 3.3 of this TEOS-10 Manual. +% +% The software is available from http://www.TEOS-10.org +% +%========================================================================== + +%-------------------------------------------------------------------------- +% Check variables and resize if necessary +%-------------------------------------------------------------------------- + +if ~(nargin==3) + error('gsw_CT_from_t: Requires three inputs') +end %if + +[ms,ns] = size(SA); +[mt,nt] = size(t); +[mp,np] = size(p); + +if (mt ~= ms | nt ~= ns) + error('gsw_CT_from_t: SA and t must have same dimensions') +end + +if (mp == 1) & (np == 1) % p scalar - fill to size of SA + p = p*ones(size(SA)); +elseif (ns == np) & (mp == 1) % p is row vector, + p = p(ones(1,ms), :); % copy down each column. +elseif (ms == mp) & (np == 1) % p is column vector, + p = p(:,ones(1,ns)); % copy across each row. +elseif (ns == mp) & (np == 1) % p is a transposed row vector, + p = p.'; % transposed then + p = p(ones(1,ms),:); % copy down each column. +elseif (ms == mp) & (ns == np) + % ok +else + error('gsw_CT_from_t: Inputs array dimensions arguments do not agree') +end %if + +%Find values that are out of range, set them to NaN. +t(p < 100 & (t > 80 | t < -12)) = NaN; +t(p >= 100 & (t > 40 | t < -12)) = NaN; + +if ms == 1 + SA = SA.'; + t = t.'; + p = p.'; + transposed = 1; +else + transposed = 0; +end + +%-------------------------------------------------------------------------- +% Start of the calculation +%-------------------------------------------------------------------------- + +pt0 = gsw_pt0_from_t(SA,t,p); +CT = gsw_CT_from_pt(SA,pt0); + +if transposed + CT = CT.'; +end + +end diff --git a/Util/getPresRelForGSW.m b/Util/getPresRelForGSW.m index 65cbec41c..59383f8ab 100644 --- a/Util/getPresRelForGSW.m +++ b/Util/getPresRelForGSW.m @@ -1,4 +1,4 @@ -function [ presRel, presName ] = getPresRelForGSW( sam ) +function [ presRel, zName, zComment ] = getPresRelForGSW( sam ) %GETPRESRELFORGSW retrieves values of pressure due to sea water in sam for % use in the Gibbs-SeaWater toolbox (TEOS-10). % @@ -17,8 +17,8 @@ % Outputs: % presRel - the pressure due to sea water data retrieved from sam for % use in GSW. -% presName - the name of the variable in sam and the method used to -% produce presRel. +% zName - the name of the variable in sam. +% zComment - the method used to produce presRel. % % Author: Guillaume Galibert % @@ -45,11 +45,34 @@ if ~isstruct(sam), error('sam must be a struct'); end if isempty(sam), return; end +qcSet = str2double(readProperty('toolbox.qc_set')); +rawFlag = imosQCFlag('raw', qcSet, 'flag'); +goodFlag = imosQCFlag('good', qcSet, 'flag'); +probGoodFlag = imosQCFlag('probablyGood', qcSet, 'flag'); +goodFlags = [rawFlag, goodFlag, probGoodFlag]; + presRel = NaN; -presName = ''; +zName = ''; +zComment = ''; presIdx = getVar(sam.variables, 'PRES'); +if presIdx + if isfield(sam.variables{presIdx}, 'flags') + % we don't want to consider PRES if it's all rubbish + if ~any(ismember(sam.variables{presIdx}.flags, goodFlags)) + presIdx = 0; + end + end +end presRelIdx = getVar(sam.variables, 'PRES_REL'); +if presRelIdx + if isfield(sam.variables{presRelIdx}, 'flags') + % we don't want to consider PRES_REL if it's all rubbish + if ~any(ismember(sam.variables{presRelIdx}.flags, goodFlags)) + presRelIdx = 0; + end + end +end isPresVar = logical(presIdx || presRelIdx); isDepthInfo = false; @@ -74,13 +97,14 @@ if isPresVar if presRelIdx > 0 presRel = sam.variables{presRelIdx}.data; - presName = 'PRES_REL'; + zName = 'PRES_REL'; else % update from a relative pressure like SeaBird computes % it in its processed files, substracting a constant value % 10.1325 dbar for nominal atmospheric pressure presRel = sam.variables{presIdx}.data - gsw_P0/10^4; - presName = 'PRES substracting a constant value 10.1325 dbar for nominal atmospheric pressure'; + zName = 'PRES'; + zComment = 'substracting a constant value 10.1325 dbar for nominal atmospheric pressure'; end else % when no pressure variable exists, we use depth information either @@ -89,7 +113,7 @@ if depthIdx > 0 % with depth data depth = sam.(depthType){depthIdx}.data; - presName = 'DEPTH'; + zName = 'DEPTH'; else % with nominal depth information @@ -111,7 +135,7 @@ end depth = sam.instrument_nominal_depth * ones(size(sam.dimensions{dimIdx}.data)); - presName = 'instrument_nominal_depth'; + zName = 'instrument_nominal_depth'; end % any depth values <= -5 are discarded (reminder, depth is @@ -134,7 +158,7 @@ else % without latitude information, we assume 1dbar ~= 1m presRel = depth; - presName = [presName ' (assuming 1 m ~ 1 dbar)']; + zComment = '(assuming 1 m ~ 1 dbar)'; end end diff --git a/imosToolbox.m b/imosToolbox.m index 48467b6ad..b28563b5c 100644 --- a/imosToolbox.m +++ b/imosToolbox.m @@ -37,7 +37,7 @@ function imosToolbox(auto, varargin) % % Set current toolbox version -toolboxVersion = ['2.5.35 - ' computer]; +toolboxVersion = ['2.5.36 - ' computer]; if nargin == 0, auto = 'manual'; end diff --git a/imosToolbox_Linux64.bin b/imosToolbox_Linux64.bin index 8ddc50716..b72a64e7b 100755 Binary files a/imosToolbox_Linux64.bin and b/imosToolbox_Linux64.bin differ diff --git a/imosToolbox_Win32.exe b/imosToolbox_Win32.exe index 2173b8ce7..c7febe767 100644 Binary files a/imosToolbox_Win32.exe and b/imosToolbox_Win32.exe differ diff --git a/imosToolbox_Win64.exe b/imosToolbox_Win64.exe index 69a2c7a35..52a722f86 100644 Binary files a/imosToolbox_Win64.exe and b/imosToolbox_Win64.exe differ