diff --git a/Code/CTD.json b/Code/CTD.json index 5f281bb..f8e2409 100644 --- a/Code/CTD.json +++ b/Code/CTD.json @@ -200,6 +200,11 @@ "long_name": "Turbidity", "standard_name": "sea_water_turbidity", "coverage_content_type": "physicalMeasurement" + }, + "sn": { + "units": "", + "long_name": "Serial Number", + "coverage_content_type": "auxiliaryInformation" } }, diff --git a/Code/Combo.json b/Code/Combo.json index 3099bfa..2ee24f1 100644 --- a/Code/Combo.json +++ b/Code/Combo.json @@ -546,6 +546,11 @@ "units": "", "long_name": "ODAS Library Version", "coverage_content_type": "auxiliaryInformation" + }, + "sn": { + "units": "", + "long_name": "Serial Number", + "coverage_content_type": "auxiliaryInformation" } }, diff --git a/Code/calc_chi.m b/Code/calc_chi.m new file mode 100644 index 0000000..3b8c40d --- /dev/null +++ b/Code/calc_chi.m @@ -0,0 +1,35 @@ +% +% Calculate chi from temperature gradients given an already calculated epsilon +% +% Aug-2024, Pat Welch, pat@mousebrains.com + +function [dInfo, tbl] = calc_chi(diss, dInfo) +arguments (Input) + diss table % Profile information + dInfo (1,:) table % Summary information about the profile +end % arguments Input +arguments (Output) + dInfo (1,:) table % pInfo with extra fields + tbl table % Tabular form of diss struct +end % arguments Output + +kappa_T = 1.4e-7; % thermal diffusivity [m^2/s] + +tbl = table(); +tbl.t = diss.t; +tbl.depth = diss.depth; + +names = string(diss.Properties.VariableNames); + +for name = names(startsWith(names, "gradT")) + gradT = diss.(name); + gradT2 = gradT.^2; + % gradT2(gradT <= 0) = NaN; + index = extractAfter(name, "gradT"); + tbl.(name) = gradT; + tbl.(sprintf("chi_dT%s_mean", index)) = 2 * kappa_T * diss.N2 .* diss.epsilonMean ./ gradT2; + for j = 1:size(diss.e,2) + tbl.(sprintf("chi_dT%s_e%d", index, j)) = 2 * kappa_T * diss.N2 .* diss.e(:,j) ./ gradT2; + end % for j +end % for name +end % calc_chi \ No newline at end of file diff --git a/Code/chi2binned.m b/Code/chi2binned.m index dd753dd..58fec1a 100644 --- a/Code/chi2binned.m +++ b/Code/chi2binned.m @@ -19,86 +19,73 @@ return; end % if ~row.qProfileOkay -warning("Chi binning not implemented"); - -% fnChi = row.fnChi; -% fnBin = fullfile(pars.chi_binned_root, append(row.name, ".mat")); -% row.fnChiBin = fnBin; -% -% if isnewer(fnBin, fnChi) -% retval = {fnBin, []}; % We want this for combining -% fprintf("%s: %s is newer than %s\n", row.name, row.fnBin, row.fnChi); -% return; -% end % if isnewer -% -% if isempty(a) -% fprintf("Loading %s\n", row.fnChi); -% a = load(row.fnChi); -% end % if isempty -% -% %% Bin the data into depth bins -% -% pInfo = a.info; -% profiles = a.profiles; -% -% fprintf("%s: Binning %d profiles\n", row.name, numel(profiles)); -% -% if pars.profile_direction == "time" % Bin in time -% binSize = seconds(pars.binChi_width); % Bin stepsize in (sec) -% keyName = "t"; -% binFunc = @bin_by_time; -% glueFunc = @glue_lengthwise; -% else % Bin by depth -% binSize = pars.binChi_width; % Bin stepsize (m) -% keyName = pars.binChi_variable; % Variable to bin by -% binFunc = @bin_by_real; -% glueFunc = @glue_widthwise; -% end % if profile_direction -% -% casts = cell(numel(profiles),1); -% for index = 1:numel(profiles) -% profile = profiles{index}; -% nE = size(profile.e, 2); -% prof2 = table(); -% for name = string(profile.Properties.VariableNames) -% sz = size(profile.(name),2); -% if ~ismatrix(profile.(name)), continue; end -% if sz == 1 -% prof2.(name) = profile.(name); -% elseif sz == nE -% for j = 1:sz -% prof2.(append(name, "_", string(j))) = profile.(name)(:,j); -% end % for j; -% end % if sz -% end % for name -% -% casts{index} = binFunc(binSize, keyName, prof2, pars.binChi_method); -% end % for index -% -% qDrop = cellfun(@isempty, casts); % This shouldn't happend -% -% if any(qDrop) -% casts = casts(~qDrop); -% pInfo = pInfo(~qDrop,:); -% end % any qDrop -% -% if isempty(casts) -% row.qProfileOkay = false; -% fprintf("%s: No usable casts found in %s\n", row.name, row.fnProf); -% return; -% end -% -% tbl = glueFunc("bin", casts); -% -% binned = struct ( ... -% "tbl", tbl, ... -% "info", pInfo); -% if isfield(a, "fp07") -% binned.fp07 = a.fp07; -% end % if isfield -% -% my_mk_directory(fnBin); -% save(fnBin, "-struct", "binned", pars.matlab_file_format); -% fprintf("%s: Saving %d profiles to %s\n", row.name, size(binned.info,1), fnBin); -% retval = {fnBin, binned}; +fnChi = row.fnChi; +fnBin = fullfile(pars.chi_binned_root, append(row.name, ".mat")); +row.fnChiBin = fnBin; + +if isnewer(fnBin, fnChi) + retval = {fnBin, []}; % We want this for combining + fprintf("%s: %s is newer than %s\n", row.name, row.fnBin, row.fnChi); + return; +end % if isnewer + +if isempty(a) + fprintf("Loading %s\n", row.fnChi); + a = load(row.fnChi); +end % if isempty + +%% Bin the data into depth bins + +pInfo = a.info; +profiles = a.profiles; + +fprintf("%s: Binning %d profiles\n", row.name, numel(profiles)); + +if pars.profile_direction == "time" % Bin in time + binSize = seconds(pars.binChi_width); % Bin stepsize in (sec) + keyName = "t"; + binFunc = @bin_by_time; + glueFunc = @glue_lengthwise; +else % Bin by depth + binSize = pars.binChi_width; % Bin stepsize (m) + keyName = pars.binChi_variable; % Variable to bin by + binFunc = @bin_by_real; + glueFunc = @glue_widthwise; +end % if profile_direction + +casts = cell(numel(profiles),1); +for index = 1:numel(profiles) + profile = profiles{index}; + prof2 = table(); + for name = string(profile.Properties.VariableNames) + if ~ismatrix(profile.(name)), continue; end + prof2.(name) = profile.(name); + end % for name + + casts{index} = binFunc(binSize, keyName, prof2, pars.binChi_method); +end % for index + +qDrop = cellfun(@isempty, casts); % This shouldn't happend + +if any(qDrop) + casts = casts(~qDrop); + pInfo = pInfo(~qDrop,:); +end % any qDrop + +if isempty(casts) + row.qProfileOkay = false; + fprintf("%s: No usable casts found in %s\n", row.name, row.fnProf); + return; +end + +tbl = glueFunc("bin", casts); + +binned = struct ( ... + "tbl", tbl, ... + "info", pInfo); + +my_mk_directory(fnBin); +save(fnBin, "-struct", "binned", pars.matlab_file_format); +fprintf("%s: Saving %d profiles to %s\n", row.name, size(binned.info,1), fnBin); +retval = {fnBin, binned}; end % chi2binned diff --git a/Code/convert2mat.m b/Code/convert2mat.m index 2fa214e..d7bc497 100644 --- a/Code/convert2mat.m +++ b/Code/convert2mat.m @@ -36,6 +36,14 @@ a = odas_p2mat(char(row.fn), p2args{:}); % extract P file contents my_mk_directory(row.fnMat); save(row.fnMat, "-struct", "a", pars.matlab_file_format); % save into a mat file + + row.sn = setupstr(a.setupfilestr, 'instrument_info', 'sn'); % Get the serial number in the instrument_info stanza + if isempty(row.sn) + row.sn = missing; + else + row.sn = string(row.sn{1}); + end % if isempty + row.qMatOkay = true; fprintf("Took %.2f seconds to convert %s\n", toc(stime), row.name); catch ME diff --git a/Code/ctd2binned.m b/Code/ctd2binned.m index 90a2779..3a5c0ee 100644 --- a/Code/ctd2binned.m +++ b/Code/ctd2binned.m @@ -144,7 +144,7 @@ fnCTD = fullfile(pars.ctd_root, append(row.name, ".mat")); row.fnCTD = fnCTD; -cInfo = row(:,["name", "t0", "tEnd"]); +cInfo = row(:,["name", "t0", "tEnd", "sn"]); cInfo = renamevars(cInfo, "tEnd", "t1"); % For NetCDF time range binned = struct("tbl", tbl, "info", cInfo); @@ -216,4 +216,4 @@ ctd.lon(q) = gps.lon(t) - dLondt .* dt; ctd.lat(q) = gps.lat(t) - dLatdt .* dt; end -end % addGPS \ No newline at end of file +end % addGPS diff --git a/Code/mat2profile.m b/Code/mat2profile.m index d4f03cc..7fa7669 100644 --- a/Code/mat2profile.m +++ b/Code/mat2profile.m @@ -157,7 +157,7 @@ profilesInfo = struct(); profilesInfo.profiles = profiles; -profilesInfo.row = row(1,["name", "date", "fClock" "t0", "t1", "tEnd"]); +profilesInfo.row = row(1,["name", "date", "fClock" "t0", "t1", "tEnd", "sn"]); profilesInfo.pInfo = profileInfo; if qFP07 profilesInfo.fp07Lags = fp07_lags; @@ -181,6 +181,7 @@ tbl = table(); tbl.name = repmat(row.name, nProfiles, 1); tbl.index = (1:nProfiles)'; +tbl.sn = repmat(row.sn, nProfiles,1); tbl.t0 = NaT(nProfiles,1); tbl.t1 = NaT(nProfiles,1); tbl.n_slow = nan(nProfiles,1); diff --git a/Code/osgl_get_netCDF.m b/Code/osgl_get_netCDF.m index 04c2811..30bec60 100644 --- a/Code/osgl_get_netCDF.m +++ b/Code/osgl_get_netCDF.m @@ -15,27 +15,21 @@ % for speed reasons % Dec-2022, Pat Welch, pat@mousebrains.com, add in compliance with CF-1.8 % metadata standards -% +% Jan-2024, Pat Welch, pat@mosuebrains.com, straighten out handling of variable names -function a = osgl_get_netCDF(fn, varargin) +function a = osgl_get_netCDF(fn, names) arguments (Input) fn string {mustBeFile} % Input filename end % arguments Input arguments (Input,Repeating) - varargin % Optional arguments + names (:,1) string % Optional variable names end % arguments Repeating arguments (Output) a struct % Output structure of data read from fn end % arguments output -names = cell(numel(varargin),1); % Pre-allocate - -for i = 1:numel(varargin) - item = string(varargin{i}); - names{i} = item(:); -end % for i - -names = string(names); % Convert to an array of strings from a cell array +names = cellfun(@(x) x(:), names, "UniformOutput", false); % Make sure any arrays are in the same orientation +names = string(vertcat(names{:})); ncid = netcdf.open(fn); % Open the netCDF file cleanUpObj = onCleanup(@() netcdf.close(ncid)); % In case of errors/warnings/... @@ -43,6 +37,7 @@ if isempty(names) [names, varids] = loadNames(ncid); else + names = unique(names); % No duplicates [names, varids] = getVarIDs(ncid, names); end % if @@ -70,6 +65,15 @@ end % osgl_get_netCDF function data = modifyCF(data, ncid, varid) +arguments (Input) + data % Input data column from NetCDF file + ncid double % NetCDF id from netcdf.open + varid double % Variable id within ncid +end % arguments Output +arguments (Output) + data % Possibly modified version of input data +end % arguments Output + qOkay = true(size(data)); % Initially assume all the data is good try @@ -137,19 +141,32 @@ end % modifyCF function data = ctConvert(data, ncid, varid, dtUnits, timeStr, qOkay) +arguments (Input) + data % Input from NetCDF file various formats + ncid double % NetCDF id from netcdf.open + varid double % Variable id within ncid + dtUnits string % Time units from units attribute (The part before "since") + timeStr string % Time reference from units attribute (The part after "since") + qOkay logical % % Which data entries are good +end % arguments Input +arguments (Output) + data % possibly datetime, if converted +end % arguments Output + try - calendar = netcdf.getAtt(ncid, varid, "calendar"); - goodCalendars = ["", "standard", "gregorian", "proleptic_gregorian"]; - if ~ismember(calendar, goodCalendars) - warning("Unsupported calendar, %s, for %s", ... - calendar, netcdf.inqVar(ncid, varid)); - return - end % if + calendar = netcdf.getAtt(ncid, varid, "calendar"); % Should be there, but somem people don't have it catch ME getReport(ME) - % Do nothing + calendar = "proleptic_gregorian"; % Default end +goodCalendars = ["", "standard", "gregorian", "proleptic_gregorian"]; +if ~ismember(calendar, goodCalendars) + warning("Unsupported calendar, %s, for %s", ... + calendar, netcdf.inqVar(ncid, varid)); + return +end % if + switch lower(dtUnits) % The part before since case {"years", "year", "y"} dt = data(qOkay) * 365.242198781; @@ -175,11 +192,16 @@ return end % switch -fmt = "^\s*(\d{4})-(\d{1,2})-(\d{1,2})(\s+\d{1,2}:\d{1,2}:\d{1,2}([.]\d*|)(\s+([A-Za-z/_]+|[+-]?\d{1,2}([:]?\d{1,2}|))|)|)\s*$"; -tfmt = "^\s+(\d{1,2}):(\d{1,2}):(\d{1,2}([.]?\d*|))(\s+([A-Za-z/_]+|[+-]?\d{1,2}([:]?\d{1,2}|))|)$"; +fmt = strjoin([ + "^(\d{4})-(\d{1,2})-(\d{1,2})", % yyyy-mm-dd + "(|(\s+|[Tt])\d{1,2}:\d{1,2}:\d{1,2}(|[.]\d*))", % (\s|T)HH:MM:SS.ssss + "(|\s+([A-Za-z/_]+|[+-]?\d{1,2}(|[:]?\d{1,2})))", % (UTC|[+-]00:00) + "\s*$"], ""); +tfmt = "^[Tt\s]+(\d{1,2}):(\d{1,2}):(\d{1,2}([.]?\d*|))(\s+([A-Za-z/_]+|[+-]?\d{1,2}([:]?\d{1,2}|))|)$"; tokens = regexp(timeStr, fmt, "tokens", "once"); if isempty(tokens) + fmt warning("Unable to parse time string, %s, for %s", ... timeStr, netcdf.inqVar(ncid, varid)); return @@ -228,6 +250,14 @@ end % ctConvert function [names, varids] = loadNames(ncid) +arguments (Input) + ncid double % NetCDF file id from netcdf.open +end % arguments Input +arguments (Output) + names (:,1) string % Variable names in NetCDF file + varids (:,1) double % Variable ids in ncid +end % arguments Output + [~, n] = netcdf.inq(ncid); % Get the number of variables names = cell(n,1); varids = zeros(size(names)) - 1; @@ -246,6 +276,15 @@ end % loadNames function [names, varids] = getVarIDs(ncid, names) +arguments (Input) + ncid double % NetCDF file id from netcdf.open + names (:,1) string % Variable names to get from the NetCDF file +end % arguments Input +arguments (Output) + names (:,1) string % Variable names to get from the NetCDF file + varids (:,1) double % Variable ids in ncid +end % arguments Output + varids = zeros(size(names)); for i = 1:numel(names) varids(i) = netcdf.inqVarID(ncid, names(i)); diff --git a/Code/profile2chi.m b/Code/profile2chi.m index 3791a72..97ff8c8 100644 --- a/Code/profile2chi.m +++ b/Code/profile2chi.m @@ -1,5 +1,6 @@ % -% For each profile, calculate dissipation estimates + +% For each profile, calculate chi estimates % % Nov-2023, Pat Welch, pat@mousebrains.com @@ -7,7 +8,7 @@ arguments (Input) row (1,:) table % row to work on profileInfo struct % Output of mat2profile - dissInfo struct % Output of profile2diss + dissInfo struct % Dissipation for this this profile pars struct % Parameters, defaults from get_info end % arguments Input arguments (Output) @@ -30,44 +31,83 @@ profileInfo = load(row.fnProf); end -if isempty(dissInfo) % we need to load dissInfo +if isempty(chiInfo) % we need to load dissInfo fprintf("Loading %s\n", row.fnDiss); - dissInfo = load(row.fnDiss); + chiInfo = load(row.fnDiss); end -warning("Chi not implemented"); -% return; -% -% stime = tic(); -% -% profiles = profileInfo.profiles; -% pInfo = profileInfo.pInfo; -% nProfiles = numel(profiles); -% -% tbl = cell(nProfiles, 1); -% dInfo = cell(size(tbl)); -% -% for index = 1:nProfiles -% [dInfo{index}, tbl{index}] = calc_diss_shear(profiles{index}, pInfo(index,:), pars); -% end % for index -% -% qEmpty = cellfun(@isempty, tbl); -% if all(qEmpty) % No dissipation estimates -% fprintf("%s: No dissipation estimates found", row.name); -% return; -% end -% -% tbl = tbl(~qEmpty); -% dInfo = dInfo(~qEmpty); -% -% dissInfo = struct(); -% dissInfo.info = vertcat(dInfo{:}); -% dissInfo.profiles = tbl; -% if isfield(profileInfo, "fp07Lags") -% dissInfo.fp07 = profileInfo.fp07Lags; -% end % if isfield -% -% my_mk_directory(fnChi); -% save(fnChi, "-struct", "chiInfo", pars.matlab_file_format); -% fprintf("Took %.2f seconds to create %s\n", toc(stime), fnChi); +stime = tic(); + +profiles = profileInfo.profiles; +pInfo = profileInfo.pInfo; +nProfiles = numel(profiles); + +dProfiles = dissInfo.profiles; +dInfo = dissInfo.info; + +if size(pInfo,1) ~= size(dInfo,1) + warning("Mismatched number of profiles in chi calculation, %d ~= %d", size(pInfo,1), size(dInfo,1)); + return; +end % if size + +tbl = cell(nProfiles, 1); +dInfo = cell(size(tbl)); + +for index = 1:nProfiles + profile = profiles{index}; + fast = profile.fast; + slow = profile.slow; + diss = dProfiles{index}; + fast.dDepth = interp1(diss.depth, diss.depth, fast.depth, "nearest", "extrap"); + fast.grp = findgroups(fast.dDepth); + a = rowfun(@(x) x(1), fast, ... + InputVariables="dDepth", ... + GroupingVariables="grp", ... + OutputVariableNames="depth"); + names = string(fast.Properties.VariableNames); + names = names(startsWith(names, "gradT")); + for name = names + a.(name) = rowfun(@(x) median(x, "omitnan"), fast, ... + InputVariables=name, ... + GroupingVariables="grp", ... + OutputFormat="uniform"); + end % for + + N2 = table(); + [N2.N2, N2.pMid] = gsw_Nsquared(slow.SA, slow.theta, slow.P_slow, slow.lat); + N2.depth = (slow.depth(1:end-1) + slow.depth(2:end)) / 2; + N2.dDepth = interp1(diss.depth, diss.depth, N2.depth, "nearest", "extrap"); + N2.grp = findgroups(N2.dDepth); + + b = rowfun(@(x) x(1), N2, ... + InputVariables="dDepth", ... + GroupingVariables="grp", ... + OutputVariableNames="depth"); + b.N2 = rowfun(@(x) median(x, "omitnan"), N2, ... + InputVariables="N2", ... + GroupingVariables="grp", ... + OutputFormat="uniform"); + + diss = innerjoin(diss, a, Keys="depth", RightVariables=names); + diss = innerjoin(diss, b, Keys="depth", RightVariables="N2"); + + [dInfo{index}, tbl{index}] = calc_chi(diss, pInfo(index,:)); +end % for index + +qEmpty = cellfun(@isempty, tbl); +if all(qEmpty) % No dissipation estimates + fprintf("%s: No dissipation estimates found", row.name); + return; +end + +tbl = tbl(~qEmpty); +dInfo = dInfo(~qEmpty); + +chiInfo = struct(); +chiInfo.info = vertcat(dInfo{:}); +chiInfo.profiles = tbl; + +my_mk_directory(fnChi); +save(fnChi, "-struct", "chiInfo", pars.matlab_file_format); +fprintf("Took %.2f seconds to create %s\n", toc(stime), fnChi); end % profile2diss \ No newline at end of file diff --git a/Examples/example0.m b/Examples/example0.m index 99511bf..31b58c4 100644 --- a/Examples/example0.m +++ b/Examples/example0.m @@ -17,6 +17,7 @@ pars = process_P_files( ... "debug", true, ... + "chi_enable", true, ... % Calculate chi "p_file_root", p_file_root, ... % Where the input .P files are located "p_file_pattern", "SN*/*", ... % Glob pattern appended to p_file_root to locate P files "output_root", output_root, ... % Where to write output to