-
Notifications
You must be signed in to change notification settings - Fork 0
/
extractSalinity.m
199 lines (155 loc) · 6.21 KB
/
extractSalinity.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
% Here, the following is done:
% - observational data from Irminger Sea is loaded
% - practical salinity (PS) is extracted for EACH MOORING
% - absolute salinity (SA) is calculated and output into a file. If these
% files already exist this method can be skipped.
% - a plot is made of SP vs SA. This is just for illustrative purposes.
% Add subfolders to path
addpath(genpath(pwd));
clear; clc; close all;
set(groot,'defaultAxesXGrid','on');
set(groot,'defaultAxesYGrid','on');
set(groot, 'defaultFigureUnits', 'centimeters', 'defaultFigurePosition', [5 5 40 20]);
%% Load data
data = load('Data/merged_hourly_unfiltered_data_20142020.mat');
time = datetime(data.time_grid,'ConvertFrom','datenum'); % time [h], datetime format
latitude = data.lat;
lon = data.lon';
dgrid = data.DGRID;
SP = data.S;
% For each of the five moorings IC1, IC2, IC3, IC4, and M1, we extract the
% levels at which data is recorded, merge data that was recorded by a
% single sensor to create a full six-year time series, and find the
% (weighted) average depth that corresponds to each sensor time series.
%% Extract IC3 Salinity
% IC3: four sensors at the following depth levels
% 1: 5,9,19
% 2: 47,48,49
% 3: 99,100
% 4: 147,155,156
SP_IC3 = squeeze(SP(4,:,:));
SP_IC3_1 = nan(1,52873); SP_IC3_2 = nan(1,52873); SP_IC3_3 = nan(1,52873);
SP_IC3_4 = nan(1,52873);
for i=1:length(time)
SP_IC3_1(i) = sum([SP_IC3(i,5),SP_IC3(i,9),SP_IC3(i,19)],'omitnan');
SP_IC3_2(i) = sum([SP_IC3(i,47),SP_IC3(i,48),SP_IC3(i,49)],'omitnan');
SP_IC3_3(i) = sum([SP_IC3(i,99),SP_IC3(i,100)],'omitnan');
SP_IC3_4(i) = sum([SP_IC3(i,147),SP_IC3(i,155),SP_IC3(i,156)],'omitnan');
end
SP_IC3_avg = cat(1,SP_IC3_1,SP_IC3_2,SP_IC3_3,SP_IC3_4);
SP_IC3_avg(SP_IC3_avg==0) = NaN;
% Weighted average depth
z_Sp_IC3 = [(dgrid(5,4) + 4*dgrid(9,4) + dgrid(19,4))/6 ...
(dgrid(47,4) + dgrid(48,4) + 4*dgrid(49,4))/6 ...
(2*dgrid(99,4) + 4*dgrid(100,4))/6 ...
(2*dgrid(147,4) + dgrid(155,4) + 3*dgrid(156,4))/6];
clear SP_IC3_1 SP_IC3_2 SP_IC3_3 SP_IC3_4;
%% Extract M1 Salinity
% M1: ten sensors at the following depth levels
% 1: 6
% 2: 11,12
% 3: 21,22
% 4: 36,37
% 5: 51,52
% 6: 72
% 7: 92
% 8: 122,123
% 9: 141
% 10: 156,157
SP_M1_1 = nan(1,52873); SP_M1_2 = nan(1,52873); SP_M1_3 = nan(1,52873);
SP_M1_4 = nan(1,52873); SP_M1_5 = nan(1,52873); SP_M1_6 = nan(1,52873);
SP_M1_7 = nan(1,52873); SP_M1_8 = nan(1,52873); SP_M1_9 = nan(1,52873);
SP_M1_10 = nan(1,52873);
SP_M1 = squeeze(SP(6,:,:));
for i=1:length(time)
SP_M1_1(i) = sum([SP_M1(i,6)],'omitnan');
SP_M1_2(i) = sum([SP_M1(i,11),SP_M1(i,12)],'omitnan');
SP_M1_3(i) = sum([SP_M1(i,21),SP_M1(i,22)],'omitnan');
SP_M1_4(i) = sum([SP_M1(i,36),SP_M1(i,37)],'omitnan');
SP_M1_5(i) = sum([SP_M1(i,51),SP_M1(i,52)],'omitnan');
SP_M1_6(i) = sum([SP_M1(i,72)],'omitnan');
SP_M1_7(i) = sum([SP_M1(i,92)],'omitnan');
SP_M1_8(i) = sum([SP_M1(i,122),SP_M1(i,123)],'omitnan');
SP_M1_9(i) = sum([SP_M1(i,141)],'omitnan');
SP_M1_10(i) = sum([SP_M1(i,156),SP_M1(i,157)],'omitnan');
end
SP_M1_avg = cat(1,SP_M1_1,SP_M1_2,SP_M1_3,SP_M1_4,SP_M1_5,SP_M1_6,SP_M1_7,SP_M1_8,SP_M1_9,SP_M1_10);
SP_M1_avg(SP_M1_avg==0) = NaN;
% Weighted average depth
z_Sp_M1 = [dgrid(6,6) ...
(3*dgrid(11,6) + 3*dgrid(12,6))/6 ...
(2*dgrid(21,6) + 4*dgrid(22,6))/6 ...
(2*dgrid(36,6) + 4*dgrid(37,6))/6 ...
(2*dgrid(51,6) + 4*dgrid(52,6))/6 ...
dgrid(72,6) ...
dgrid(92,6) ...
(2*dgrid(122,6) + 4*dgrid(123,6))/6 ...
dgrid(141,6) ...
(2*dgrid(156,6) + 4*dgrid(157,6))/6];
clear SP_M1_1 SP_M1_2 SP_M1_3 SP_M1_4 SP_M1_5 SP_M1_6 SP_M1_7 SP_M1_8 SP_M1_9 SP_M1_10;
%% Save Practical Salinity from IC3 and M1
SP_IC3 = SP_IC3_avg;
SP_M1 = SP_M1_avg;
clear SP_IC3_avg SP_M1_avg;
save Matfiles/Sp.mat SP_IC3 SP_M1;
%% Practical Salinity: fill in missing values
SP_IC3f = fillMissingSalinity(4,SP_IC3',time);
SP_M1f = fillMissingSalinity(10,SP_M1',time);
%% Depths: save for later
save Matfiles/depthsS z_Sp_IC3 z_Sp_M1;
%% Pressures: needed for Sp->SA conversion
p_Sp_IC3 = sw_pres(z_Sp_IC3,latitude(4))';
p_Sp_M1 = sw_pres(z_Sp_M1,latitude(6))';
save Matfiles/p_S.mat p_Sp_IC3 p_Sp_M1;
%% Convert Practical Salinity (SP) into Absolute Salinity (SA)
% Note that the convertIntoAbsoluteSalinity function will save a mat file
% according to the filename entered.
% SA must be found before CT because CT is dependent on SA.
%% SP3->SA3
% Have now removed SP_IC3f, hopefully it still works with the unfilled
% data
if isfile('SA_IC3.mat')
disp('SA already calculated for IC3');
load SA_IC3.mat;
else
disp('Calculating SA for IC3');
SA_IC3 = convertIntoAbsoluteSalinity(SP_IC3,time,p_Sp_IC3,lon(4),latitude(4),"Matfiles/SA_IC3");
end
%% SP-M1 -> SA-M1
if isfile('SA_M1.mat')
disp('SA already calculated for M1');
load SA_M1.mat;
else
disp('Calculating SA for M1');
SA_M1 = convertIntoAbsoluteSalinity(SP_M1,time,p_Sp_M1,lon(6),latitude(6),"Matfiles/SA_M1");
end
%% SP -> SA: Filled Data
disp('Start conversion: filled SP-M1...');
SA_M1f = convertIntoAbsoluteSalinity(SP_M1f,time,p_Sp_M1,lon(6),latitude(6),"Matfiles/SA_M1f");
disp('Start conversion: filled SP-IC3...');
SA_IC3f = convertIntoAbsoluteSalinity(SP_IC3f,time,p_Sp_IC3,lon(4),latitude(4),'Matfiles/SA_IC3f');
%% Visualise
% SA_XXs = smoothed (two week mean),
% SA_XXv = vertical average (of all salinity measurement locations).
SA_M1s = nan(10,52873);
SA_IC3s = nan(4,52873);
for i = 1:10
SA_M1s(i,:) = movmean(SA_M1(i,:),24*14);
end
for i = 1:4
SA_IC3s(i,:) = movmean(SA_IC3(i,:),24*14);
end
SA_M1v = mean(SA_M1s,1,'omitnan');
SA_IC3v = mean(SA_IC3s,1,'omitnan');
ax1 = figure;
plot(time,SA_M1s,'Color','#a6cee3','Marker','.','MarkerSize',1.5,'HandleVisibility','off');
hold on
plot(time,SA_IC3s,'Color','#b2df8a','Marker','.','MarkerSize',2,'HandleVisibility','off');
plot(time,SA_M1v,'Color','#1f78b4','LineWidth',2,'DisplayName','M1');
plot(time,SA_IC3v,'Color','#33a02c','LineWidth',2,'DisplayName','IC3');
hold off
legend(); title('Time Evolution of $S_A$ at Moorings IC3 and M1','Interpreter','latex','FontSize',16);
ylabel('Absolute Salinity $S_A$','Interpreter','latex'); xlabel('Time','Interpreter','latex');
exportgraphics(ax1,'figures/main/extract/0_SA.png');
%% Clear unnecessary variables
clear latitude lon i dgrid ax1 SP SP_IC3 SP_M1;