-
Notifications
You must be signed in to change notification settings - Fork 0
/
runall_tropicalpacificonlygridded.m
244 lines (212 loc) · 9.58 KB
/
runall_tropicalpacificonlygridded.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
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
% This runall_tropicalpacificonlygridded.m MATLAB script
% loads and computes variables for re-creating
% figures in Leung et al. (2019).
% It runs the following figure creation scripts...
% mainfig1andsuppfig1de.m
% mainfig2andsuppfig4.m
% mainfig3.m
% mainfig4.m
% suppfig1abc.m
% suppfig2.m
% suppfig3.m
% suppfig6.m
% suppfig8.m
% suppfig9.m
% ...Thus, creating main figures 1-4 and
% suppplementary figures 1-3,6,8-9.
% Figure creation scripts with multiple figure
% names in their titles generate all figures in the
% title, but only one at a time, which must be specified
% (e.g., to generate Supp. Fig. 4, you will have to
% go into mainfig2andsuppfig4.m and specify that you
% want Supp. Fig. 4 rather than Main Fig. 2 plotted.
% Same for mainfig1andsuppfig1de.m).
clear all;
addpath(genpath('.'));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1.) Define and load spatial grid, time variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
begdate = datetime(1955,1,1); % beginning date of WOD dataset
enddate = datetime(2017,5,1); % end date of WOD dataset
numdegres = 5; % horiz grid resolution in degrees
initmonum = 1; % WOD dataset's initial month (jan=1,feb=2,etc.)
%--------------------------
% Load the following space/time grid variables:
% (note that "tponly" mean Tropical Pacific only)
% lon (0-360), lat, depth (m), timemoWOD, timeyrWOD
%--------------------------
load('data/WODtropicalpacificonlygridded/WOD_gridandtime_195501-201705_tponly_5deg_700m.mat','lon','lat','depth','timemoWOD','timeyrWOD');
%--------------------------
% Load/calculate the following ENSO index (ONI) variables:
% oni (degC), onitime,
% onien (0 = non-El Nino month, 1 = El Nino month),
% oniln (0 = non-La Nina month, 1 = La Nina month)
%--------------------------
loadoni;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2.) Load and/or compute gridded variables from WOD data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------------
% Load the following variables:
% temp (degC), sal (psu), po2 (kPa)
%--------------------------
load('data/WODtropicalpacificonlygridded/WOD_tempsalpo2_195501-201705_tponly_5deg_700m.mat','temp','sal','po2');
%--------------------------
% Compute monthly anomalies.
% (note that "anwod" = monthly anomalies from WOD climatology)
% Get the following variables:
% tempanwod (degC), salanwod (psu), po2anwod (kPa)
% po250 (po2 at 50 m depth, kPa),
% po2100 (po2 at 100 m depth, kPa)
%--------------------------
tempclimwod = calcWODclimatology(temp,initmonum);
tempanwod = calcmonthlyanoms(temp,tempclimwod,initmonum);
salclimwod = calcWODclimatology(sal,initmonum);
salanwod = calcmonthlyanoms(sal,salclimwod,initmonum);
po2climwod = calcWODclimatology(po2,initmonum);
po2anwod = calcmonthlyanoms(po2,po2climwod,initmonum);
po250 = squeeze(po2(:,:,find(depth==50),:));
po2100 = squeeze(po2(:,:,find(depth==100),:));
%--------------------------
% Compute tuna hypoxic depth.
% Get the following variables:
% tunahypoxicpressure (kPa),
% thd (tuna hypoxic depth, m), thdanwod (m),
%--------------------------
load('data/WODtropicalpacificonlygridded/WOD_tunahypoxicdepth_195501-201705_tponly_5deg_700m.mat','thd','thdinfo','tunahypoxicpressure');
% --> Code to calculate the variables saved in WOD_tunahypoxicdepth_195501-201705_tponly_5deg_700m.mat:
%tunahypoxicpressure = 15; % kPa
%[thd,thdinfo] = calcthd(tunahypoxicpressure,po2,depth);
%save('data/WODtropicalpacificonlygridded/WOD_tunahypoxicdepth_195501-201705_tponly_5deg_700m.mat','thd','thdinfo','tunahypoxicpressure','-v7.3');
% Use only pO2 profiles where computing thd makes most sense
% (i.e., profiles shouldn't start inside an OMZ, there
% shouldn't be highly fluctuating pO2 values that would make
% deciding where thd really is difficult, etc.)
indxs = sub2ind(size(thd),...
thdinfo.mtabovebelowprofidxs(:,1),...
thdinfo.mtabovebelowprofidxs(:,2),...
thdinfo.mtabovebelowprofidxs(:,3));
newthd = nan(size(thd));
newthd(indxs) = thd(indxs);
thd = newthd; clear newthd;
[thdclimwod,~] = calcthd(tunahypoxicpressure,po2climwod,depth);
thdanwod = calcmonthlyanoms(thd,thdclimwod,initmonum);
%--------------------------
% Compute thermocline depth.
% Get the following variables:
% tcdtype (thermocline depth calculation method),
% tcd (thermocline depth, m), tcdanwod (m)
%--------------------------
load('data/WODtropicalpacificonlygridded/WOD_thermoclinedepth_195501-201705_tponly_5deg_700m.mat','tcd','tcdinfo','tcdtype');
% --> Code to calculate the variables saved in WOD_thermoclinedepth_195501-201705_tponly_5deg_700m.mat:
%tcdtype='varit';
%[tcd,tcdinfo] = calctcd(tcdtype,temp,depth);
%save('data/WODtropicalpacificonlygridded/WOD_thermoclinedepth_195501-201705_tponly_5deg_700m.mat','tcd','tcdinfo','tcdtype');
% Use only temp profiles where computing tcd makes most sense
indxs = sub2ind(size(tcd),...
tcdinfo.mtdecprofidxs(:,1),...
tcdinfo.mtdecprofidxs(:,2),...
tcdinfo.mtdecprofidxs(:,3));
newtcd = nan(size(tcd));
newtcd(indxs) = tcd(indxs);
tcd = newtcd; clear newtcd;
[tcdclimwod,tcdinfo] = calctcd(tcdtype,tempclimwod,depth);
tcdanwod = calcmonthlyanoms(tcd,tcdclimwod,initmonum);
%--------------------------
% Compute correlations.
% Get the following variable:
% thdanwodvstcdanwod (correlation coefficients
% and p-values between thdanwod and tcdanwod)
%--------------------------
[thdanwodvstcdanwod.tempccmap,~,~,~,~,~,thdanwodvstcdanwod.pvalmap]=tempcorrmapnanwithintandslopeuncertainty(tcdanwod,thdanwod);
%--------------------------
% Compute means and ENSO composites.
% (note that "enoni" = monthly El Nino composite anomalies,
% "lnoni" = monthly La Nina composite anomalies)
% Get the following variables:
% tcdmean (m), tcdanwod_enoni (m), tcdanwod_lnoni (m)
%--------------------------
tcdmean = nanmean(tcd,3);
[tcdanwodmean,tcdanwod_enoni,tcdanwod_lnoni] = calcmeanandENLNcompmaps(tcdanwod,onien,oniln);
%--------------------------
% Regrid from depth space to isopycnal space.
% (note that "rho" = variable is in isopycnal space)
% Get the following variables:
% rho (kg/m^3), po2rho (kPa),
% po2rhoanwod (kPa), tcdrho (m)
%--------------------------
load('data/WODtropicalpacificonlygridded/WOD_varsinrhospace_195501-201705_tponly_5deg_700m.mat','rhoedges','rho','po2rho','po2rhoanwod','tcdrho');
% --> Code to calculate the variables saved in WOD_varsinrhospace_195501-201705_tponly_5deg_700m.mat:
%[pottemp_refsurf,potdens_refsurf]=calcpottempanddens(temp,sal,depth,0);
%potdens_refsurfclimwod = calcWODclimatology(potdens_refsurf,initmonum);
%rhoedges = [1020:0.2:1022.8 1023:1:1026 1026.5 1027 1027.2:0.2:1028];
%rho = (rhoedges(1:end-1)+rhoedges(2:end))/2;
%po2rho = calcrhospace4d(po2,rhoedges,potdens_refsurf);
%po2rhoclimwod = calcrhospace4d(po2climwod,rhoedges,potdens_refsurfclimwod);
%po2rhoanwod = calcmonthlyanoms(po2rho,po2rhoclimwod,initmonum);
%tcdrho = calcrhofromdepth(tcd,depth,potdens_refsurf);
%save('data/WODtropicalpacificonlygridded/WOD_varsinrhospace_195501-201705_tponly_5deg_700m.mat','rhoedges','rho','po2rho','po2rhoanwod','tcdrho','-v7.3');
%--------------------------
% Count how much data there is.
% (note that "np" = number of data points,
% "tot" = total)
% Get the following variables (all w/ units of # months):
% po2nptot, po2np_enoni, po2np_lnoni,
% thdnptot, thdnp_enoni, thdnp_lnoni,
% tcdnptot, tcdnp_enoni, tcdnp_lnoni
%--------------------------
[po2nptot,po2np_enoni,po2np_lnoni] = calcmeanandENLNcompnumpts(po2,onien,oniln);
[thdnptot,thdnp_enoni,thdnp_lnoni] = calcmeanandENLNcompnumpts(thd,onien,oniln);
[tcdnptot,tcdnp_enoni,tcdnp_lnoni] = calcmeanandENLNcompnumpts(tcd,onien,oniln);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3.) Prepare helpful things for plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------------
% Set up for plotting EEZ outlines.
% Get the following variables:
% allpid, nupid, eezmap0pt25, loneez0pt25, lateez0pt25
%--------------------------
% Load in EEZ shapefile attributes to get Country-PolygonID association
wlon = 100; elon = 300; slat = -30; nlat = 30;
[s,a] = shaperead('data/EEZs/EEZshapefile/World_EEZ_v9_2016_HR_0_360','UseGeoCoords',true,'BoundingBox',[wlon,slat;elon,nlat]);
c = struct2cell(a);
allctry = c(10,:); allpid = c(1,:);
nuctries = {'Micronesia','Kiribati','Marshall Islands','Nauru','Palau','Papua New Guinea','Solomon Islands','Tuvalu'};
nupid = cell(length(nuctries),1);
for ictry = 1:length(nuctries)
cidx = find(strcmp(allctry,nuctries{ictry}));
nupid{ictry} = allpid(cidx);
end
% Load nc files of gridded EEZs w/ PolygonID as identifier
fnamenc = 'data/EEZs/World_EEZ_v9_2016_HR_0_360_30Nto30S_polyid_0pt25deg.nc';
ncid=netcdf.open(fnamenc,'nc_nowrite');
varid = netcdf.inqVarID(ncid,'z'); % dims: lon, lat
eezmap0pt25 = netcdf.getVar(ncid,varid,'double');,
eezmap0pt25 = eezmap0pt25'; % x = lon, y = lat
netcdf.close(ncid);
numdegreseez = 0.25;
loneez0pt25 = 0:numdegreseez:360;
lateez0pt25 = (-30:numdegreseez:30)';
%--------------------------
% Shift axes for proper plotting w/ pcolor
% (essentially defines bottom/left border of grid cell).
% Get the following variables:
% plon, plat, pdepth (m), prho (kg/m^3)
%--------------------------
plat = lat-numdegres/2;
plon = lon-numdegres/2;
pdepth = [(depth(1)-depth(2))/2; depth(2:end)-diff(depth/2)];
prho = rhoedges(1:end-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4.) Create and save figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mainfig1andsuppfig1de; % go to figurescripts/mainfig1andsuppfig1de.m to specify the figure you want to create
mainfig2andsuppfig4; % go to figurescripts/mainfig2andsuppfig4.m to specify the figure you want to create
mainfig3;
mainfig4;
suppfig1abc;
suppfig2;
suppfig3;
suppfig6;
suppfig8;
suppfig9;