/
tmd_data.m
263 lines (231 loc) · 10.7 KB
/
tmd_data.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
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
function [Z,x_or_lon,y_or_lat,conList] = tmd_data(filename,variable,varargin)
% tmd_data loads tide model data into the Matlab workspace.
%
%% Syntax
%
% [Z,x_or_lon,y_or_lat] = tmd_data(filename,variable)
% [...] = tmd_data(...,'constituents',conList)
% [...] = tmd_data(...,'bounds',[xi yi])
% [Z,lon,lat] = tmd_data(...,'geo')
%
%% Description
%
% [Z,x_or_lon,y_or_lat] = tmd_data(filename,variable) loads any of the following
% variables from a given tide model file along with projected x,y model
% coordinates:
%
% * 'h' complex tidal height (m)
% * 'hRe' real part of tidal height
% * 'hIm' imaginary part of tidal height
% * 'hAm' amplitude of tidal height
% * 'hPh' phase of tidal height
% * 'u' complex zonal velocity (m/s)
% * 'uRe' real part of zonal velocity
% * 'uIm' imaginary part of zonal velocity
% * 'uAm' amplitude of zonal velocity
% * 'uPh' phase of zonal velocity
% * 'U' complex zonal transport (m^2/s)
% * 'URe real part of zonal transport
% * 'UIm' imaginary part of zonal transport
% * 'UAm' amplitude of zonal transport
% * 'UPh' phase of zonal transport
% * 'v' complex meridional velocity (m/s)
% * 'vRe' real part of meridional velocity
% * 'vIm' imaginary part of meridional velocity
% * 'vAm' amplitude of meridional velocity
% * 'vPh' phase of meridional velocity
% * 'V' complex meridional transport (m^2/s)
% * 'VRe' real part of meridional transport
% * 'VIm' imaginary part of meridional transport
% * 'VAm' amplitude of meridional transport
% * 'VPh' phase of meridional transport
% * 'wct' water column thickness (m)
% * 'mask' binary land/ocean mask
% * 'flexure' ice shelf flexure coefficient from a linear elastic model applied to BedMachine ice thickness (can slightly exceed 1). Only for CATS model.
%
% [...] = tmd_data(...,'constituents',conList) specifies tidal constituents as a
% cell array (e.g, {'m2','s2'}. If constituents are not specified, all constituents
% from the model are returned.
%
% [...] = tmd_data(...,'bounds',[xi yi]) enter an Mx2 matrix of coordinates
% to only load data in a rectangle around the specified [xi(:) yi(:)] where
% xi and yi are projected coordinates for tide models that are in projected
% coordinates. For tide models in geo coordinates (global tide models,
% generally), xi refers to longitudes of interest and yi is latitudes of interest.
% The advantage of specifying bounds is to minimize the amount of data that
% are loaded, when only a small area within a larger tide model is of
% interest.
%
% [Z,lon,lat] = tmd_data(...,'geo') returns grid coordinates as
% geographic coordinates. This is the default behavior for global models,
% whereas regional models return projected coordinates by default.
%
% [...,cons] = tmd_data(...) returns a list of constituents in the model.
%
%% Examples
% For examples type
%
% tmd tmd_data
%
%% Author Info
% This function was written by Chad A. Greene in 2022.
%
% See also tmd_interp and tmd_predict.
%% Error checks
assert(contains(filename,'.nc'),'Input filename must end in .nc.')
assert(exist(filename,'file'),['Cannot find ',filename,'. Check the path and try again.'])
% Ensure the model file is TMD3.0 compatible:
try
tmd_version = ncreadatt(filename,'/','tmd_version');
assert(tmd_version>=3.0,[filename,' is not compatible with TMD3.0+.'])
catch
error([filename,' is not compatible with TMD3.0+.'])
end
%% Input parsing
narginchk(2,Inf)
% Set defaults
spatialSubset = false;
conSubset = false;
geo = false; % output geographic coordinates
stride = Inf;
NCons = 1;
bounds = [[-Inf;Inf] [-Inf;Inf]];
conList = strsplit(ncreadatt(filename,'constituents','constituent_order'));
proj4 = ncreadatt(filename,'mapping','spatial_proj4');
if strcmpi(proj4,'+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
GlobalModel = true;
else
GlobalModel = false;
end
if nargin>2
tmp = strncmpi(varargin,'bounds',3);
if any(tmp)
spatialSubset = true;
bounds = varargin{find(tmp)+1};
assert(size(bounds,2)==2,'Spatial bounds must be Mx2 in the form [xi(:) yi(:)] or [loni(:) lati(:)].')
% Wrap longitudes to 360 for global model:
if GlobalModel
tmp=bounds(:,1);
tmp(tmp<0) = tmp(tmp<0)+360;
bounds(:,1) = tmp;
end
end
tmp = strncmpi(varargin,'constituents',3);
if any(tmp)
conSubset = true;
assert(all(ismember(varargin{find(tmp)+1},conList)),'Error: Some of the requested constituents do not exist in the model file.')
conList = varargin{find(tmp)+1};
if ischar(conList)
conList = cellstr(conList);
end
stride = 1;
NCons = length(conList);
end
if any(strcmpi(varargin,'geo'))
geo = true;
end
end
%% Load data
if GlobalModel
x_or_lon = double(ncread(filename,'lon'));
y_or_lat = double(ncread(filename,'lat'));
else
x_or_lon = double(ncread(filename,'x'));
y_or_lat = double(ncread(filename,'y'));
end
cons = strsplit(ncreadatt(filename,'constituents','constituent_order'));
dx = abs(diff(x_or_lon(1:2)));
dy = abs(diff(y_or_lat(1:2)));
if spatialSubset
% Get xlimits (xl) and ylimits (yl) of input coordinates + buffer:
xl = [min(bounds(:,1))-2*dx max(bounds(:,1))+2*dx];
yl = [min(bounds(:,2))-2*dy max(bounds(:,2))+2*dy];
% Region of rows and columns of pixels to read:
ri=find((y_or_lat>=yl(1))&(y_or_lat<=yl(2)));
ci=find((x_or_lon>=xl(1))&(x_or_lon<=xl(2)));
x_or_lon = x_or_lon(ci);
y_or_lat = y_or_lat(ri);
else
ri = 1:length(y_or_lat);
ci = 1:length(x_or_lon);
end
if conSubset
Z = NaN(numel(ri),numel(ci),NCons);
end
% In case there's no data:
assert(numel(ci)>0,'No tide data available in the region of interest.')
for k = 1:NCons
if conSubset
conind = find(ismember(cons,conList{k}));
placInd = k;
else
conind = 1;
placInd = 1:length(cons);
end
switch variable
case 'mask'
Z = logical(permute(ncread(filename,variable,[ci(1) ri(1)],[numel(ci) numel(ri)]),[2 1]));
case 'flexure'
Z = double(permute(ncread(filename,variable,[ci(1) ri(1)],[numel(ci) numel(ri)]),[2 1]))/100;
case 'wct'
Z = double(permute(ncread(filename,variable,[ci(1) ri(1)],[numel(ci) numel(ri)]),[2 1]));
case {'hRe','hIm','URe','UIm','VRe','VIm'}
Z(:,:,placInd) = double(permute(ncread(filename,variable,[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3]));
case 'uRe'
Z(:,:,placInd) = double(permute(ncread(filename,'URe',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3]));
case 'uIm'
Z(:,:,placInd) = double(permute(ncread(filename,'UIm',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3]));
case 'vRe'
Z(:,:,placInd) = double(permute(ncread(filename,'VRe',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3]));
case 'vIm'
Z(:,:,placInd) = double(permute(ncread(filename,'VIm',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3]));
case 'h'
Z(:,:,placInd) = complex(double(permute(ncread(filename,'hRe',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3])),...
double(permute(ncread(filename,'hIm',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3])));
case 'hAm'
Z(:,:,placInd) = abs(complex(double(permute(ncread(filename,'hRe',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3])),...
double(permute(ncread(filename,'hIm',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3]))));
case 'hPh'
Z(:,:,placInd) = angle(complex(double(permute(ncread(filename,'hRe',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3])),...
double(permute(ncread(filename,'hIm',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3]))));
case {'u','U'}
Z(:,:,placInd) = complex(double(permute(ncread(filename,'URe',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3])),...
double(permute(ncread(filename,'UIm',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3])));
case {'uAm','UAm'}
Z(:,:,placInd) = abs(complex(double(permute(ncread(filename,'URe',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3])),...
double(permute(ncread(filename,'UIm',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3]))));
case {'uPh','UPh'}
Z(:,:,placInd) = angle(complex(double(permute(ncread(filename,'URe',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3])),...
double(permute(ncread(filename,'UIm',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3]))));
case {'v','V'}
Z(:,:,placInd) = complex(double(permute(ncread(filename,'VRe',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3])),...
double(permute(ncread(filename,'VIm',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3])));
case {'vAm','VAm'}
Z(:,:,placInd) = abs(complex(double(permute(ncread(filename,'VRe',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3])),...
double(permute(ncread(filename,'VIm',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3]))));
case {'vPh','VPh'}
Z(:,:,placInd) = angle(complex(double(permute(ncread(filename,'VRe',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3])),...
double(permute(ncread(filename,'VIm',[ci(1) ri(1) conind],[numel(ci) numel(ri) stride]),[2 1 3]))));
otherwise
error(['The requested variable ',variable,' is not a variable in file ',filename,'.'])
end
end
% Convert transports to velocities if requested:
if ismember(variable,{'uRe','uIm','u','uAm','vRe','vIm','v','vAm'})
wct = double(permute(ncread(filename,'wct',[ci(1) ri(1)],[numel(ci) numel(ri)]),[2 1]));
Z = Z./max(wct,10); % divide transports by at least 10 m wct to prevent insanely high velocities.
% % ALTERNATE APPROACH: Ensure velocities don't get insanely high by dividing transports by at least 10 m, and here's where we *could* NaN out land for currents.
% mask = double(permute(ncread(filename,'mask',[ci(1) ri(1)],[numel(ci) numel(ri)]),[2 1]));
% wct(wct<10) = 10;
% wct(~mask) = nan;
% Z = Z./wct;
end
if geo
if GlobalModel
% do nothing, x_or_lon is already lon and y_or_lat is already lat.
else
x_or_lon = double(permute(ncread(filename,'lon',[ci(1) ri(1)],[numel(ci) numel(ri)]),[2 1]));
y_or_lat = double(permute(ncread(filename,'lat',[ci(1) ri(1)],[numel(ci) numel(ri)]),[2 1]));
end
end
end