url = 'https://www.aoml.noaa.gov/ocd/gcc/GOMECC2/GOMECC2_station_data_version4.csv';
filename = 'GOMECC2_station_data_version4.csv';  % you can specify a full path here if needed
websave(filename, url);
Goemec2 = readtable(filename);

%

load("\\atlas.shore.mbari.org\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\Spray_Data\25420901\25420901sat.mat")

%

tic
% Change depth to be positive
s.depth = -s.depth;
%% This section is for calculating ESPER for a full deployment
% Prepare vars
pHin_ESPER = NaN(size(s.pHin)); % Preallocate array for estimated pH
nprofiles = length(s.sdn);

for profile = 1:nprofiles
    try
        DesiredVariables = 3; % pH insitu (total scale) (from function)
        depth = s.depth(:,profile);
        nsamps = length(depth);
        lon = repmat(s.lon(profile),nsamps,1); % Size matches depth column
        lat = repmat(s.lat(profile),nsamps,1); % Size matches depth column
        OutputCoordinates = [lon, lat, depth]; % Lon, Lat, Depth
        
        
        PredictorMeasurements = [s.psal(:,profile), s.tc(:,profile), s.doxy(:,profile)]; % Salinity, Temperature, Dissolved Oxygen
        PredictorTypes = [1 2 6]; % Salinity, Temperature, Dissolved Oxygen (from function)
        eqs = 7; % Salinity, Temperature, Dissolved Oxygen (equation from function)
        
        % Convert serial date number to decimal year
        [year, month, day, hour, minute, second] = datevec(s.sdn(profile));
        start_of_year = datenum(year, 1, 1);
        end_of_year = datenum(year + 1, 1, 1);
        % Calculate the decimal year
        dec_yr = year + (s.sdn(profile) - start_of_year) ./ (end_of_year - start_of_year);
        dec_yr = repmat(dec_yr,nsamps,1); % Size matches depth column
        
        [Estimates,~,~]= ESPER_LIR(DesiredVariables,OutputCoordinates,PredictorMeasurements,PredictorTypes, ... 
                    'Equations',eqs,'EstDates', dec_yr); % ESPER model to predict pH
        
        pHin_ESPER(:,profile) = Estimates.pH;
        disp(['Calculating Profile ', num2str(profile), ' / ', num2str(nprofiles)])
    catch
        disp(['Skipping Profile ', num2str(profile), ' / ', num2str(nprofiles)])
        continue;
    end
end
toc
MissionID = '25402901';
% Save variable to .mat file
fname = ['pHin_Esper_', MissionID,'.mat'];
save(fname,'pHin_ESPER','-mat') % save variable

%

% need everything in vectors
ta = s.ta_canb(:);
pHin_esper = pHin_ESPER(:);
tc = s.tc(:);
psal = s.psal(:);
pres = s.pres(:);
sil = s.sil_canb(:);
po4 = s.po4_canb(:); 

% CO2SYS(PAR1,PAR2,PAR1TYPE,PAR2TYPE,SAL,TEMPIN,TEMPOUT,PRESIN,PRESOUT,SI,PO4,NH4,H2S,pHSCALEIN,K1K2CONSTANTS,KSO4CONSTANT,KFCONSTANT,BORON)
trex = CO2SYSv3(ta, pHin_esper, 1, 3, psal, tc, 25, pres, 0, sil, po4, 0, 0, 1, 10, 1, 2, 2);
% extract variables
s.pH25atm_ESPER = reshape(trex(:,20), size(s.pres)); % changed from 18 to 20 following CO2SYSv3 update
s.pHin_ESPER = pHin_ESPER;

%

save([MissionID, '_ESPER_CANB_42525.mat'],'s','-mat') % save variable

%
% Investigating GOEMEC vs GliderFET, ESPER, and CanyonB below 500 m
%

load("pHin_Esper_25402901.mat");
load("25402901_ESPER_CANB_42525.mat");
filename = 'GOMECC2_station_data_version4.csv';  % you can specify a full path here if needed
Goemec2 = readtable(filename);

%

figure
geobasemap('satellite')
geoscatter(Goemec2.LATITUDE,Goemec2.LONGITUDE,Goemec2.STNNBR + 1,'cyan','filled','o')
hold on;
geoscatter(s.lat,s.lon,'red','filled','o')
legend('GOEMEC 2', 'Spray2',Location='northwest')

%
% Nice crossovers on line 5. Here testing ESPER skill to predict pHtot25atm at 1000m. It is pretty spot on. Should use that instead of canyonB to correct pH in gulf stream. Our k0 is way off. Need to fix that in post!
%

idx_cast50 = Goemec2.STNNBR >= 52 & Goemec2.STNNBR <= 52 & Goemec2.PH_TOT ~= -999;
pHtot_cast50 = Goemec2.PH_TOT(idx_cast50); % First intersect with spray2 (pH in total scale @ 25 C)
depth_cast50 = Goemec2.DEPTH(idx_cast50);

idx_spray = 15;
SpraypH25atm = s.pH25atm(:,idx_spray);
SpraypH25atm(SpraypH25atm == -999) = NaN;
SprayDepth = -s.depth(:,idx_spray);
SprayDepth(SprayDepth == -999) = NaN;

pH25atm_ESPER = s.pH25atm_ESPER(:,idx_spray);
pH25atm_ESPER(pH25atm_ESPER == -999) = NaN;

offset = .72;
% offset = 0;

figure
scatter(pHtot_cast50,depth_cast50,'cyan','filled','o')
hold on;
scatter(SpraypH25atm+offset,SprayDepth,'red','filled','o')
scatter(pH25atm_ESPER,SprayDepth,'yellow','filled','o')


set(gca,'Ydir','reverse')
grid on;
legend('GOEMEC','SPRAY2 + 0.72 offset','pH25atm_ESPER',Location='southwest')

figure
geobasemap('satellite')
geoscatter(Goemec2.LATITUDE(idx_cast50),Goemec2.LONGITUDE(idx_cast50),'cyan','filled','o')
hold on;
geoscatter(s.lat(idx_spray),s.lon(idx_spray),'red','filled','o')
legend('GOEMEC 2', 'Spray2',Location='northwest')

%
% Lets do some hovmollers for offshore distance vs depth and pH25atm. crossovers between 47 and 56
%

% Convert COEMEC2 pH25atm to pH insitu for comparison with CanB, ESPER, and
% Spray
% ALKALI: Alkalinity – UMOL/KG
% PH_TOT: pH in total scale @ 25 C
% SILCAT: Silicate – UMOL/KG
% PHSPHT: Phosphate – UMOL/KG
% TCO2 and pH
trex = CO2SYSv3(Goemec2.ALKALI, Goemec2.PH_TOT, 1, 3, Goemec2.CTDSAL, 25, Goemec2.CTDTMP, 0, Goemec2.CTDPRS, Goemec2.SILCAT, Goemec2.PHSPHT, 0, 0, 1, 10, 1, 2, 2); % convert from pH25atm to pHin
Goemec2.pHin = reshape(trex(:,20), size(Goemec2.CTDPRS));
Goemec2.pHin(Goemec2.pHin == -999) = NaN;

%

% Compare

% GOEMEC2
idxstart = 50; idxend = 53;
idx_cast = Goemec2.STNNBR >= idxstart & Goemec2.STNNBR <= idxend & Goemec2.PH_TOT ~= -999;
pHtot_Goemec2 = Goemec2.PH_TOT(idx_cast); 
pHin_Goemec2 = Goemec2.pHin(idx_cast); 
depth_cast = Goemec2.DEPTH(idx_cast);

% Spray
idx_spray = [10:12, 28:40];
SpraypHin = s.pHin(:,idx_spray);
SprayDepth = -s.depth(:,idx_spray);
SpraypHin = SpraypHin(:);
SprayDepth = SprayDepth(:);
% CanyonB
pHin_CanyonB = s.pHin_canb(:,idx_spray);
pHin_CanyonB = pHin_CanyonB(:);
% ESPER
pHin_ESPER = s.pHin_ESPER(:,idx_spray);
pHin_ESPER = pHin_ESPER(:);
figure
scatter(pHin_Goemec2,depth_cast,'cyan','o','filled');
hold on;
scatter(pHin_CanyonB,SprayDepth,'magenta','.');
scatter(pHin_ESPER,SprayDepth,'blue','.');
scatter(SpraypHin,SprayDepth,'black','.');

% colororder("reef")
set(gca,'Ydir','reverse');
grid on;
legend('GOEMEC2','CanyonB','ESPER','GliderFET003','Location','southwest')
title('Comparison all nearby data: pHin')

figure
geobasemap('satellite')
geoscatter(Goemec2.LATITUDE(idx_cast),Goemec2.LONGITUDE(idx_cast),'cyan','filled','o')
hold on;
geoscatter(s.lat(idx_spray),s.lon(idx_spray),'red','filled','o')
% legend('GOEMEC 2', 'Spray2',Location='northwest')

%CALCDISTANCESPRAYGOEMEc Calculate distance between Spray profile and BOG
% wgs84 = wgs84Ellipsoid("m");
% distance2spray = distance(Goemec2.LATITUDE(idx_cast),Goemec2.LONGITUDE(idx_cast),s.lat(idx_spray),s.lon(idx_spray),wgs84) / 1000; % km
% title(['Distance: ', num2str(distance2spray(1)), ' km'])

%

% Compare

% GOEMEC2
idxstart = 51; idxend = 51;
idx_cast = Goemec2.STNNBR >= idxstart & Goemec2.STNNBR <= idxend & Goemec2.PH_TOT ~= -999;
pHtot_Goemec2 = Goemec2.PH_TOT(idx_cast); 
pHin_Goemec2 = Goemec2.pHin(idx_cast); 
depth_cast = Goemec2.DEPTH(idx_cast);

% Spray
idx_spray = [12]; % Xover at 1000 m
SpraypHin = s.pHin(:,idx_spray);
SprayDepth = -s.depth(:,idx_spray);
SpraypHin = SpraypHin(:);
SprayDepth = SprayDepth(:);
CanyonB
pHin_CanyonB = s.pHin_canb(:,idx_spray);
pHin_CanyonB = pHin_CanyonB(:);
% ESPER
pHin_ESPER = s.pHin_ESPER(:,idx_spray);
pHin_ESPER = pHin_ESPER(:);
figure
scatter(pHin_Goemec2,depth_cast,'.');
hold on;
scatter(pHin_CanyonB,SprayDepth,'.');
scatter(pHin_ESPER,SprayDepth,'.');
scatter(SpraypHin,SprayDepth,'.');

% colororder("reef")
set(gca,'Ydir','reverse');
grid on;
legend('GOEMEC2','CanyonB','ESPER','Location','southwest')
title('Comparison: pHin')

figure
geobasemap('satellite')
geoscatter(Goemec2.LATITUDE(idx_cast),Goemec2.LONGITUDE(idx_cast),'cyan','filled','o')
hold on;
geoscatter(s.lat(idx_spray),s.lon(idx_spray),'red','filled','o')
legend('GOEMEC 2', 'Spray2',Location='northwest')

%CALCDISTANCESPRAYGOEMEc Calculate distance between Spray profile and BOG
wgs84 = wgs84Ellipsoid("m");
distance2spray = distance(Goemec2.LATITUDE(idx_cast),Goemec2.LONGITUDE(idx_cast),s.lat(idx_spray),s.lon(idx_spray),wgs84) / 1000; % km
title(['Distance: ', num2str(distance2spray(1)), ' km'])