Extract tidal harmonic phases and amplitudes from POLPRED ASCII files.
get_POLPRED_spectide(Mobj, POLPRED)
DESCRIPTION:
Extract POLPRED harmonic amplitude and phases for the nearest point in
the POLPRED grid to the open boundary nodes in Mobj.
INPUT:
Mobj = MATLAB mesh object (see read_sms_mesh.m)
POLPRED = ASCII file of the POLPRED harmonics
OUTPUT:
Mobj = MATLAB mesh object with two new arrays:
phase - Harmonic phases at each open boundary point
amp - Harmonic amplitudes at each open boundary point
EXAMPLE USAGE
Mobj = get_POLPRED_spectide(Mobj, '/path/to/POLPRED.txt')
Author(s):
Pierre Cazenave (Plymouth Marine Laboratory)
Revision history
2012-11-15 First version. Based in part on tide_tools.py from the
fvcom-py Python toolbox (https://bitbucket.org/pwcazenave/fvcom-py)
and Ricardo Torres' searchtides.m.
==========================================================================0001 function [Mobj] = get_POLPRED_spectide(Mobj, POLPRED) 0002 % Extract tidal harmonic phases and amplitudes from POLPRED ASCII files. 0003 % 0004 % get_POLPRED_spectide(Mobj, POLPRED) 0005 % 0006 % DESCRIPTION: 0007 % Extract POLPRED harmonic amplitude and phases for the nearest point in 0008 % the POLPRED grid to the open boundary nodes in Mobj. 0009 % 0010 % INPUT: 0011 % Mobj = MATLAB mesh object (see read_sms_mesh.m) 0012 % POLPRED = ASCII file of the POLPRED harmonics 0013 % 0014 % OUTPUT: 0015 % Mobj = MATLAB mesh object with two new arrays: 0016 % phase - Harmonic phases at each open boundary point 0017 % amp - Harmonic amplitudes at each open boundary point 0018 % 0019 % EXAMPLE USAGE 0020 % Mobj = get_POLPRED_spectide(Mobj, '/path/to/POLPRED.txt') 0021 % 0022 % Author(s): 0023 % Pierre Cazenave (Plymouth Marine Laboratory) 0024 % 0025 % Revision history 0026 % 2012-11-15 First version. Based in part on tide_tools.py from the 0027 % fvcom-py Python toolbox (https://bitbucket.org/pwcazenave/fvcom-py) 0028 % and Ricardo Torres' searchtides.m. 0029 % 0030 %========================================================================== 0031 0032 subname = 'get_POLPRED_spectide'; 0033 0034 global ftbverbose; 0035 if ftbverbose 0036 fprintf('\n') 0037 fprintf(['begin : ' subname '\n']) 0038 end 0039 0040 % Check we have spherical coordinates in Mobj (we can't extract harmonics 0041 % without them (easily)). 0042 if ~Mobj.have_lonlat 0043 error('Spherical coordinates absent from Mobj. Cannot extract harmonics from cartesian coordinates.') 0044 end 0045 0046 % Read the POLPRED header into a struct of header names plus their values. 0047 fid = fopen(POLPRED,'rt'); 0048 if(fid < 0) 0049 error(['file: ' POLPRED ' does not exist']); 0050 end 0051 0052 if ftbverbose 0053 fprintf(['reading from: ' POLPRED '\n']) 0054 fprintf('extracting header\n') 0055 end 0056 0057 readingHeader = true; 0058 header = struct(); 0059 while readingHeader 0060 lin = fgetl(fid); 0061 if isempty(lin) 0062 % Got to the end of the header 0063 readingHeader = false; 0064 else 0065 % We have some header information. Load it into a struct. 0066 key = regexp(lin, ':', 'split'); 0067 header.(strtrim(regexprep(key{1}, ' ', '_'))) = strtrim(key{2}); 0068 end 0069 end 0070 0071 % Reformat the list of harmonics into a more sensible format 0072 header.Harmonics = regexp(header.Harmonics, ' ', 'split'); 0073 0074 % Get the positions in header.Harmonics for the constituents in which we're 0075 % interested. 0076 0077 pInd = 1:length(header.Harmonics); 0078 pIndUse = nan(length(Mobj.Components), 2); 0079 for i=1:length(Mobj.Components) 0080 pPos = pInd(strcmp(Mobj.Components{i}, header.Harmonics)); 0081 if isempty(pPos) 0082 warning('Supplied constituent (%s) is not present in the POLPRED data', Mobj.Components{i}) %#ok<WNTAG> 0083 else 0084 % Make index start at zero so the multiplication works, but 0085 % compensate for that once the offset has been applied. Also add 0086 % offset for the 2 columns (amplitude and phase). 0087 pIndUse(i, :) = (repmat((pPos - 1 ) * 6, 1, 2) + 1) + (0:1); 0088 end 0089 end 0090 % Add three to offset by the lat, lon and flag columns 0091 pIndUse = pIndUse + 3; 0092 0093 % Now we're at the data. Load it all into a massive array. 0094 if ftbverbose 0095 fprintf('extracting data\n') 0096 end 0097 0098 readingData = true; 0099 i = 0; 0100 % Preallocate data to something big and then cut back afterwards (if 0101 % necessary). Get the number of columns from the header and multiply by 6 0102 % (amplitude and phases for z, u and v). Add three for the lat, lon and 0103 % flag columns). The rows is the number of data lines in my files for the 0104 % northwest European shelf domain. 0105 nCols = 3 + (str2double(header.Number_of_harmonics) * 6); 0106 data = nan(358797, nCols); 0107 if ftbverbose 0108 tic 0109 end 0110 while readingData 0111 lin = fgetl(fid); 0112 if lin ~= -1 % EOF is -1 0113 i = i + 1; 0114 if ftbverbose 0115 if mod(i, 10000) == 0 0116 fprintf('line %i\n', i) 0117 end 0118 end 0119 % str2double doesn't work without a couple of calls to regexp, 0120 % which makes it ~20x slower than str2num on its own. The regexp 0121 % approach is still here if you don't believe me. 0122 data(i, :) = str2num(strtrim(lin)); %#ok<ST2NM> 0123 % data(i, :) = str2double(regexp(regexprep(strtrim(lin), '\s+', ' '), ' ', 'split')); 0124 else 0125 if ftbverbose 0126 fprintf('end of file at line %i\n', i) 0127 end 0128 readingData = false; 0129 end 0130 end 0131 if ftbverbose 0132 toc 0133 end 0134 0135 fclose(fid); 0136 0137 % Clear out any additional NaNs in data from preallocation. 0138 data = reshape(data(~isnan(data)), i, nCols); 0139 0140 % Now we have the data, identify the indices of data which correspond to 0141 % the nearest point to each open boundary point. This approach may not be 0142 % the best: it might instead be better to simply read in the positions and 0143 % create an index which we use to extract the harmonics of interest. 0144 % However, we've got this far so might as well carry on. 0145 0146 % Get the open boundary node IDs with which to extract their locations 0147 tmpObcNodes = Mobj.obc_nodes'; 0148 ObcNodes = tmpObcNodes(tmpObcNodes~=0)'; 0149 obc_lon = Mobj.lon(ObcNodes); 0150 obc_lat = Mobj.lat(ObcNodes); 0151 0152 % For each position, find the nearest POLPRED value. Use the 0153 % find_nearest_pt.m logic to get the nearest point (we can't use the 0154 % function here because the values for which we're searching aren't in 0155 % Mobj). 0156 distance = nan(size(obc_lon)); 0157 point = nan(size(distance)); 0158 % Omit the NaNs in the indices from POLPRED when calculating the output 0159 % array size. 0160 amp = nan(length(obc_lon), length(pIndUse(~isnan(pIndUse(:, 1)), 1))); 0161 phase = nan(size(amp)); 0162 for i=1:length(obc_lon) 0163 radvec = sqrt((obc_lon(i)-data(:,2)).^2 + (obc_lat(i)-data(:,1)).^2); 0164 [distance(i), point(i)] = min(radvec); 0165 % Get the amplitude and phase for each constituent (in order of 0166 % Mobj.Components). Check for and omit NaNs here (for missing tidal 0167 % constituents in the supplied list and what's given in POLPRED). 0168 amp(i, :) = data(point(i), pIndUse(~isnan(pIndUse(:, 1)), 1)); 0169 phase(i, :) = data(point(i), pIndUse(~isnan(pIndUse(:, 1)), 2)); 0170 end 0171 0172 % Check for and warn about NaNs (-999.9 in POLPRED data). 0173 if sum(amp(:)==-999.9) > 0 0174 warning('NaN values (-999.9 in POLPRED terms) in the amplitude data. Are your boundaries on land?') %#ok<WNTAG> 0175 end 0176 if sum(phase(:)==-999.9) > 0 0177 warning('NaN values (-999.9 in POLPRED terms) in the phase data. Are your boundaries on land?') %#ok<WNTAG> 0178 end 0179 0180 Mobj.amp = amp; 0181 Mobj.phase = phase; 0182 0183 % Plot the open boundary positions and the closest POLPRED point. 0184 % figure(1000) 0185 % plot(obc_lon, obc_lat, 'o') 0186 % hold on 0187 % plot(data(point,2), data(point,1), 'rx') 0188