%% Read in table properly fpath = 'C:\Users\bwerb\Documents\CUGNROMS\IOOS glider data 2020'; files = dir(fullfile(fpath, '*.csv')); % filter to CSVs, or use '*.*' for all % i = 5; for i = 5:15 fprintf('file %s', string(i)); fname = files(i).name; if contains(fname, '_ce_') opts = detectImportOptions(fname); opts = setvaropts(opts, opts.VariableNames{1}, 'Type', 'char'); % unnamed index opts = setvaropts(opts, opts.VariableNames{2}, 'Type', 'char'); % time (UTC) opts = setvaropts(opts, opts.VariableNames{3}, 'Type', 'double'); % depth (m) opts = setvaropts(opts, opts.VariableNames{4}, 'Type', 'double'); % latitude opts = setvaropts(opts, opts.VariableNames{5}, 'Type', 'double'); % longitude opts = setvaropts(opts, opts.VariableNames{6}, 'Type', 'double'); % temperature opts = setvaropts(opts, opts.VariableNames{7}, 'Type', 'double'); % salinity opts = setvaropts(opts, opts.VariableNames{8}, 'Type', 'double'); % dissolved_oxygen opts = setvaropts(opts, opts.VariableNames{9}, 'Type', 'double'); % density % opts.SelectedVariableNames = opts.VariableNames(2:end); opts.VariableNames = {'idx', 'time_UTC_','depth','latn','lone','tempc','psal','do','sigma'}; T_RAW = readtable(fname, opts); T_RAW.time_UTC_ = datetime(T_RAW.time_UTC_, 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss''Z''', 'TimeZone', 'UTC'); T_RAW.mon_day_yr = T_RAW.time_UTC_; T_RAW.mon_day_yr.Format = 'MM/dd/uuuu'; T_RAW.hh_mm = T_RAW.time_UTC_; T_RAW.hh_mm.Format = 'HH:mm'; T = table(); T.WMO_ID = repmat({fname(9:11)}, height(T_RAW), 1); T.Prof_num = NaN(height(T_RAW), 1); T.mon_day_yr = T_RAW.mon_day_yr; T.hh_mm = T_RAW.hh_mm; T.Lon_E = T_RAW.lone; T.Lat_N = T_RAW.latn; T.Position_QF = NaN(height(T_RAW), 1); T.Pressure_dbar = gsw_p_from_z(-T_RAW.depth, T_RAW.latn); T.Pressure_dbar_QF = NaN(height(T_RAW), 1); T.Temperature_C = T_RAW.tempc; T.Temperature_C_QF = NaN(height(T_RAW), 1); T.Salinity_pss = T_RAW.psal; T.Salinity_pss_QF = NaN(height(T_RAW), 1); T.Sigma_theta_kg_m3 = T_RAW.sigma; T.Sigma_theta_kg_m3_QF = NaN(height(T_RAW), 1); T.Depth_m = T_RAW.depth; T.Depth_m_QC = NaN(height(T_RAW), 1); T.Oxygen_umol_kg = T_RAW.do; T.Oxygen_umol_kg_QF = NaN(height(T_RAW), 1); T.OxygenSat = NaN(height(T_RAW), 1); T.OxygenSat_QF = NaN(height(T_RAW), 1); T.Nitrate_umol_kg = NaN(height(T_RAW), 1); T.Nitrate_umol_kg_QF = NaN(height(T_RAW), 1); T.pHinsitu_Total = NaN(height(T_RAW), 1); T.pHinsitu_Total_QF = NaN(height(T_RAW), 1); T.TALK_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.DIC_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.pH_ALGORITHM_insitu_total = NaN(height(T_RAW), 1); T.NO3_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.PO4_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.SILICATE_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.DIC_pHTA_umol_kg = NaN(height(T_RAW), 1); T.PCO2_pHTA_uatm = NaN(height(T_RAW), 1); elseif contains(fname, 'gp') % gp-specific readtable setup opts = detectImportOptions(fname); opts = setvaropts(opts, opts.VariableNames{1}, 'Type', 'char'); % unnamed index opts = setvaropts(opts, opts.VariableNames{2}, 'Type', 'char'); % time (UTC) opts = setvaropts(opts, opts.VariableNames{3}, 'Type', 'double'); % depth (m) opts = setvaropts(opts, opts.VariableNames{4}, 'Type', 'double'); % latitude opts = setvaropts(opts, opts.VariableNames{5}, 'Type', 'double'); % longitude opts = setvaropts(opts, opts.VariableNames{6}, 'Type', 'double'); % temperature opts = setvaropts(opts, opts.VariableNames{7}, 'Type', 'double'); % salinity opts = setvaropts(opts, opts.VariableNames{8}, 'Type', 'double'); % dissolved_oxygen opts = setvaropts(opts, opts.VariableNames{9}, 'Type', 'double'); % density % opts.SelectedVariableNames = opts.VariableNames(2:end); opts.VariableNames = {'idx', 'time_UTC_','depth','latn','lone','tempc','psal','do','sigma'}; T_RAW = readtable(fname, opts); T_RAW.time_UTC_ = datetime(T_RAW.time_UTC_, 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss''Z''', 'TimeZone', 'UTC'); T_RAW.mon_day_yr = T_RAW.time_UTC_; T_RAW.mon_day_yr.Format = 'MM/dd/uuuu'; T_RAW.hh_mm = T_RAW.time_UTC_; T_RAW.hh_mm.Format = 'HH:mm'; T = table(); T.WMO_ID = repmat({fname(9:11)}, height(T_RAW), 1); T.Prof_num = NaN(height(T_RAW), 1); T.mon_day_yr = T_RAW.mon_day_yr; T.hh_mm = T_RAW.hh_mm; T.Lon_E = T_RAW.lone; T.Lat_N = T_RAW.latn; T.Position_QF = NaN(height(T_RAW), 1); T.Pressure_dbar = gsw_p_from_z(-T_RAW.depth, T_RAW.latn); T.Pressure_dbar_QF = NaN(height(T_RAW), 1); T.Temperature_C = T_RAW.tempc; T.Temperature_C_QF = NaN(height(T_RAW), 1); T.Salinity_pss = T_RAW.psal; T.Salinity_pss_QF = NaN(height(T_RAW), 1); T.Sigma_theta_kg_m3 = T_RAW.sigma; T.Sigma_theta_kg_m3_QF = NaN(height(T_RAW), 1); T.Depth_m = T_RAW.depth; T.Depth_m_QC = NaN(height(T_RAW), 1); T.Oxygen_umol_kg = T_RAW.do; T.Oxygen_umol_kg_QF = NaN(height(T_RAW), 1); T.OxygenSat = NaN(height(T_RAW), 1); T.OxygenSat_QF = NaN(height(T_RAW), 1); T.Nitrate_umol_kg = NaN(height(T_RAW), 1); T.Nitrate_umol_kg_QF = NaN(height(T_RAW), 1); T.pHinsitu_Total = NaN(height(T_RAW), 1); T.pHinsitu_Total_QF = NaN(height(T_RAW), 1); T.TALK_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.DIC_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.pH_ALGORITHM_insitu_total = NaN(height(T_RAW), 1); T.NO3_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.PO4_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.SILICATE_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.DIC_pHTA_umol_kg = NaN(height(T_RAW), 1); T.PCO2_pHTA_uatm = NaN(height(T_RAW), 1); elseif contains(fname, '_osu') % osu-specific readtable setup opts = detectImportOptions(fname); opts = setvaropts(opts, opts.VariableNames{1}, 'Type', 'char'); % unnamed index opts = setvaropts(opts, opts.VariableNames{2}, 'Type', 'char'); % time (UTC) opts = setvaropts(opts, opts.VariableNames{3}, 'Type', 'double'); % depth (m) opts = setvaropts(opts, opts.VariableNames{4}, 'Type', 'double'); % latitude opts = setvaropts(opts, opts.VariableNames{5}, 'Type', 'double'); % longitude opts = setvaropts(opts, opts.VariableNames{6}, 'Type', 'double'); % temperature opts = setvaropts(opts, opts.VariableNames{7}, 'Type', 'double'); % salinity opts = setvaropts(opts, opts.VariableNames{8}, 'Type', 'double'); % dissolved_oxygen opts = setvaropts(opts, opts.VariableNames{9}, 'Type', 'double'); % density % opts.SelectedVariableNames = opts.VariableNames(2:end); opts.VariableNames = {'idx', 'time_UTC_','depth','latn','lone','tempc','psal','do','sigma'}; T_RAW = readtable(fname, opts); T_RAW.time_UTC_ = datetime(T_RAW.time_UTC_, 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss''Z''', 'TimeZone', 'UTC'); T_RAW.mon_day_yr = T_RAW.time_UTC_; T_RAW.mon_day_yr.Format = 'MM/dd/uuuu'; T_RAW.hh_mm = T_RAW.time_UTC_; T_RAW.hh_mm.Format = 'HH:mm'; T = table(); T.WMO_ID = repmat({fname(9:11)}, height(T_RAW), 1); T.Prof_num = NaN(height(T_RAW), 1); T.mon_day_yr = T_RAW.mon_day_yr; T.hh_mm = T_RAW.hh_mm; T.Lon_E = T_RAW.lone; T.Lat_N = T_RAW.latn; T.Position_QF = NaN(height(T_RAW), 1); T.Pressure_dbar = gsw_p_from_z(-T_RAW.depth, T_RAW.latn); T.Pressure_dbar_QF = NaN(height(T_RAW), 1); T.Temperature_C = T_RAW.tempc; T.Temperature_C_QF = NaN(height(T_RAW), 1); T.Salinity_pss = T_RAW.psal; T.Salinity_pss_QF = NaN(height(T_RAW), 1); T.Sigma_theta_kg_m3 = T_RAW.sigma; T.Sigma_theta_kg_m3_QF = NaN(height(T_RAW), 1); T.Depth_m = T_RAW.depth; T.Depth_m_QC = NaN(height(T_RAW), 1); T.Oxygen_umol_kg = T_RAW.do; T.Oxygen_umol_kg_QF = NaN(height(T_RAW), 1); T.OxygenSat = NaN(height(T_RAW), 1); T.OxygenSat_QF = NaN(height(T_RAW), 1); T.Nitrate_umol_kg = NaN(height(T_RAW), 1); T.Nitrate_umol_kg_QF = NaN(height(T_RAW), 1); T.pHinsitu_Total = NaN(height(T_RAW), 1); T.pHinsitu_Total_QF = NaN(height(T_RAW), 1); T.TALK_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.DIC_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.pH_ALGORITHM_insitu_total = NaN(height(T_RAW), 1); T.NO3_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.PO4_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.SILICATE_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.DIC_pHTA_umol_kg = NaN(height(T_RAW), 1); T.PCO2_pHTA_uatm = NaN(height(T_RAW), 1); elseif contains(fname, '_UW') % A. Quality control values % https://iop.apl.washington.edu/Seaglider_Quality_Control_Manual.html % These are the available quality control names and their numeric equivalents. They are taken from Argo, 2010 with the addition of QC_UNSAMPLED. Not all values are currently used. % % QC_NO_CHANGE 0 - No QC was performed % QC_GOOD 1 - Value is ok % QC_PROBABLY_GOOD 2 - Value is likely good % QC_PROBABLY_BAD 3 - Potentially correctable % QC_BAD 4 - Untrustworthy and uncorrectable % QC_CHANGED 5 - Explicit manual change % QC_UNSAMPLED 6 - Explicitly not sampled (vs. expected but QC_MISSING) % QC_INTERPOLATED 8 - Interpolated value % QC_MISSING 9 - Value missing; instrument timed out opts = detectImportOptions(fname); opts = setvaropts(opts, opts.VariableNames{1}, 'Type', 'char'); % unnamed index opts = setvaropts(opts, opts.VariableNames{2}, 'Type', 'char'); % time (UTC) opts = setvaropts(opts, opts.VariableNames{3}, 'Type', 'double'); % depth (m) opts = setvaropts(opts, opts.VariableNames{4}, 'Type', 'double'); % latitude opts = setvaropts(opts, opts.VariableNames{5}, 'Type', 'double'); % longitude opts = setvaropts(opts, opts.VariableNames{6}, 'Type', 'double'); % temperature opts = setvaropts(opts, opts.VariableNames{7}, 'Type', 'double'); % salinity opts = setvaropts(opts, opts.VariableNames{8}, 'Type', 'double'); % dissolved_oxygen opts = setvaropts(opts, opts.VariableNames{9}, 'Type', 'int8'); % temp_QC opts = setvaropts(opts, opts.VariableNames{9}, 'Type', 'int8'); % psal_QC opts = setvaropts(opts, opts.VariableNames{9}, 'Type', 'int8'); % do_QC opts = setvaropts(opts, opts.VariableNames{9}, 'Type', 'double'); % sigma % opts.SelectedVariableNames = opts.VariableNames(2:end); opts.VariableNames = {'idx', 'time_UTC_','depth','latn','lone','tempc','psal','do','temp_QC','psal_QC','do_QC','sigma'}; T_RAW = readtable(fname, opts); T_RAW.time_UTC_ = datetime(T_RAW.time_UTC_, 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss''Z''', 'TimeZone', 'UTC'); T_RAW.mon_day_yr = T_RAW.time_UTC_; T_RAW.mon_day_yr.Format = 'MM/dd/uuuu'; T_RAW.hh_mm = T_RAW.time_UTC_; T_RAW.hh_mm.Format = 'HH:mm'; QC_idx = find(T_RAW.temp_QC == 1 & T_RAW.psal_QC == 1 & T_RAW.do_QC == 1); T_RAW = T_RAW(QC_idx,:); T = table(); T.WMO_ID = repmat({fname(8:10)}, height(T_RAW), 1); T.Prof_num = NaN(height(T_RAW), 1); T.mon_day_yr = T_RAW.mon_day_yr; T.hh_mm = T_RAW.hh_mm; T.Lon_E = T_RAW.lone; T.Lat_N = T_RAW.latn; T.Position_QF = NaN(height(T_RAW), 1); T.Pressure_dbar = gsw_p_from_z(-T_RAW.depth, T_RAW.latn); T.Pressure_dbar_QF = NaN(height(T_RAW), 1); T.Temperature_C = T_RAW.tempc; T.Temperature_C_QF = T_RAW.temp_QC; T.Salinity_pss = T_RAW.psal; T.Salinity_pss_QF = T_RAW.psal_QC; T.Sigma_theta_kg_m3 = T_RAW.sigma; T.Sigma_theta_kg_m3_QF = NaN(height(T_RAW), 1); T.Depth_m = T_RAW.depth; T.Depth_m_QC = NaN(height(T_RAW), 1); T.Oxygen_umol_kg = T_RAW.do; T.Oxygen_umol_kg_QF = T_RAW.do_QC; T.OxygenSat = NaN(height(T_RAW), 1); T.OxygenSat_QF = NaN(height(T_RAW), 1); T.Nitrate_umol_kg = NaN(height(T_RAW), 1); T.Nitrate_umol_kg_QF = NaN(height(T_RAW), 1); T.pHinsitu_Total = NaN(height(T_RAW), 1); T.pHinsitu_Total_QF = NaN(height(T_RAW), 1); T.TALK_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.DIC_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.pH_ALGORITHM_insitu_total = NaN(height(T_RAW), 1); T.NO3_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.PO4_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.SILICATE_ALGORITHM_umol_kg = NaN(height(T_RAW), 1); T.DIC_pHTA_umol_kg = NaN(height(T_RAW), 1); T.PCO2_pHTA_uatm = NaN(height(T_RAW), 1); else error('Unrecognized file format: %s', fname) end clear T_RAW; %% Compute ESPER sdn = datenum(T.mon_day_yr + timeofday(T.hh_mm)); DesireVar = [1, 2, 3, 4, 5, 6]; %TA, DIC, pH, phosphate, nitrate, silicate OutCoords = horzcat(T.Lon_E, T.Lat_N, T.Depth_m); %lon, lat, depth PredictorTypes = [1 2 6]; % PSAL, TEMP, DOXY_ADJ, MeasEsper = [T.Salinity_pss, T.Temperature_C, T.Oxygen_umol_kg]; refyear = year(sdn)+month(sdn)/12; % reference year for OA adjustment Equations = 7; % S, T, O2 % ESPER-LIR for TA, DIC, pH, and Nitrate disp('Starting ESPER_LIR...'); [EspLir,~,~] = ESPER_LIR(DesireVar, OutCoords, MeasEsper,... PredictorTypes, 'Equations', Equations, 'EstDates', refyear); disp('Finished ESPER_LIR'); % save variables to 's' structure s.ta_esplir = EspLir.TA; s.dic_esplir = EspLir.DIC; s.pH_esplir = EspLir.pH; s.no3_esplir = EspLir.nitrate; s.po4_esplir = EspLir.phosphate; s.sioh4_esplir = EspLir.silicate; % Esper NN disp('Starting ESPER_NN...'); [EspNN,~] = ESPER_NN(DesireVar, OutCoords, MeasEsper,... PredictorTypes, 'Equations', Equations, 'EstDates', refyear); disp('Finished ESPER_NN'); % save variables to 's' structure s.ta_espnn = EspNN.TA; s.dic_espnn = EspNN.DIC; s.pH_espnn = EspNN.pH; s.no3_espnn = EspNN.nitrate; s.po4_espnn = EspNN.phosphate; s.sioh4_espnn = EspNN.silicate; clear DesireVar OutCoords PredictorTypes MeasEsper refyear Equations %% CANB % canb = CANYONB(sdn, T.Lat_N, T.Lon_E, T.Pressure_dbar, T.Temperature_C, T.Salinity_pss, T.Oxygen_umol_kg, {'NO3', 'AT', 'CT', 'pH', 'SiOH4', 'PO4'}); % make sure same weights % % s.pH_canb = canb.pH + (canb.pH*0.0404 - 0.3168); % additional correction to make it consistent to 'spec pH' % s.ta_canb = canb.AT; % s.dic_canb = canb.CT; % s.no3_canb = canb.NO3; % s.po4_canb = canb.PO4; % s.sioh4_canb = canb.SiOH4; %% CANB fast % --- find rows with no NaNs across all CANYONB inputs --- measCanb = [T.Pressure_dbar, T.Temperature_C, T.Salinity_pss, T.Oxygen_umol_kg]; validIdx = find(all(~isnan(measCanb), 2)); nTotal = height(T); nValid = length(validIdx); fprintf('Running CANYONB on %d of %d rows\n', nValid, nTotal); % --- run CANYONB on valid rows only --- canb = CANYONB(... sdn(validIdx), ... T.Lat_N(validIdx), ... T.Lon_E(validIdx), ... T.Pressure_dbar(validIdx), ... T.Temperature_C(validIdx), ... T.Salinity_pss(validIdx), ... T.Oxygen_umol_kg(validIdx), ... {'NO3', 'AT', 'CT', 'pH', 'SiOH4', 'PO4'}); %% % --- reindex outputs to full table size --- pH_canb_raw = expandToFull(canb.pH, validIdx, nTotal); s.pH_canb = pH_canb_raw + (pH_canb_raw * 0.0404 - 0.3168); s.ta_canb = expandToFull(canb.AT, validIdx, nTotal); s.dic_canb = expandToFull(canb.CT, validIdx, nTotal); s.no3_canb = expandToFull(canb.NO3, validIdx, nTotal); s.po4_canb = expandToFull(canb.PO4, validIdx, nTotal); s.sioh4_canb = expandToFull(canb.SiOH4, validIdx, nTotal); %% Add to table: compute the average between ESPERLIR,ESPERNN,CANB T.TALK_ALGORITHM_umol_kg = mean(cat(3,s.ta_esplir, s.ta_espnn, s.ta_canb), 3,'omitnan'); T.DIC_ALGORITHM_umol_kg = mean(cat(3,s.dic_esplir, s.dic_espnn, s.dic_canb),3,'omitnan'); T.pH_ALGORITHM_insitu_total = mean(cat(3,s.pH_esplir, s.pH_espnn, s.pH_canb), 3,'omitnan'); T.NO3_ALGORITHM_umol_kg = mean(cat(3,s.no3_esplir, s.no3_espnn, s.no3_canb),3,'omitnan'); T.PO4_ALGORITHM_umol_kg = mean(cat(3,s.po4_esplir, s.po4_espnn, s.po4_canb),3,'omitnan'); T.SILICATE_ALGORITHM_umol_kg = mean(cat(3,s.sioh4_esplir, s.sioh4_espnn, s.sioh4_canb),3,'omitnan'); %% Now use CO2SYS with TA_ALGORITHM, pH, nutrients to estimate DIC and PCO2. % need everything in vectors ta = T.TALK_ALGORITHM_umol_kg; pHin = T.pHinsitu_Total; tc = T.Temperature_C; psal = T.Salinity_pss; pres = T.Pressure_dbar; sil = T.SILICATE_ALGORITHM_umol_kg; po4 = T.PO4_ALGORITHM_umol_kg; trex = CO2SYSv3(ta, pHin, 1, 3, psal, tc, 25, pres, 0, sil, po4, 0, 0, 1, 10, 1, 2, 2); % extract variables T.DIC_pHTA_umol_kg = trex(:,2); T.PCO2_pHTA_uatm = trex(:,4); clear ta pHin tc psal pres sil po4 trex %% -10000000000, -999 -> nan T = standardizeMissing(T, [-10000000000, -999]);%% save files %% SAVE fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS','IOOS glider data 2020',... 'IOOS_DERIVED_PARAMS',[fname(1:end-4), '_ROMS.csv']); % Define metadata lines metadata = { '// IOOS GLIDER WITH DERIVED PARAMETERS' '// Date: 2026-03-27' '// Estimated parameters are average of CANYONB, ESPER-LIR, and ESPER-NN' '// Uncertainty: TALK_ALGORITHM_umol_kg = 15 umol/kg, DIC_ALGORITHM_umol_kg = 20 umol/kg, pH_ALGORITHM_insitu_total = 0.05, NO3_ALGORITHM_umol_kg = 3 umol/kg, pHinsitu_Total = 0.01, Nitrate_umol_kg = 1 umol/kg, DIC_pHTA_umol_kg = 10 umol/kg (DIC calculated from measured pH + algorithm TA), Oxygen_umol_kg = 1%' }; % Open file and write metadata fid = fopen(fsave, 'w'); for i = 1:length(metadata) fprintf(fid, '%s\n', metadata{i}); end fclose(fid); writetable(T,fsave,"WriteMode","append",WriteVariableNames=true); clear fsave; end %% % --- helper (if not already defined) --- function out = expandToFull(vals, idx, n) out = nan(n, 1); out(idx) = vals; end