grid_vert_interp.m 6.86 KB
Newer Older
1
function dataz = grid_vert_interp(Mobj, lon, lat, data, depth, mask, varargin)
2 3 4 5 6 7 8 9 10
% Child function to interpolate a given 3D array (x by y by sigma) values
% to the unstructured grid vertical layers.
%
% function grid_vert_interp(Mobj, lon, lat, data, depth, mask)
%
% DESCRIPTION:
%    Interpolate the regularly gridded data described by lon, lat, data
%    (whose size should be x by y for the lon and lat arrays, and x by y by
%    z for the data array) to the unstructured grid vertical layers defined
11 12 13
%    in Mobj.siglayz. The depths in depth are used to make sure the profile
%    is interpolated the right way up. The closest unstructured grid node
%    is used for the vertical interpolation of each regularly gridded node.
14 15 16
%
% INPUT:
%   Mobj        = MATLAB mesh structure which must contain:
17 18
%                   - Mobj.have_lonlat - boolean for spherical coordinate
%                   fields (Mobj.lon, Mobj.lat) presence (true = yes).
19 20 21 22 23 24 25
%                   - Mobj.siglayz - sigma layer depths for all model
%                   nodes.
%   lon, lat    = Rectangular arrays of longitude and latitude (see
%   meshgrid).
%   data, depth = x by y by z (where z is vertical layers) grids of the
%   data and water depths to be interpolated onto the vertical grid defined
%   by Mobj.siglayz.
26
%   mask        = logical array of positions outside the regularly gridded
27 28 29
%   domain (e.g. if the regular data contains NaNs or other undefined
%   values, create a logical array of those positions so they can be
%   omitted quickly).
30 31 32 33
%   'extrapolate' [optional] = keyword-argument pair in which the argument
%   specifies the coordinates to use for the extrapolation (e.g.
%   'extrapolate', [Mobj.lonc, Mobj.latc] to extrapolate onto the element
%   centres). Defaults to the element nodes (i.e. [Mobj.lon, Mobj.lat]).
34
%
35 36 37 38 39
% OUTPUT:
%   x by y by z array of vertically interpolated values at each regular
%   grid location.
%
% EXAMPLE USAGE
40 41 42 43 44 45
%   Basic usage:
%       grid_vert_interp(Mobj, lon, lat, data, depth, mask)
%
%   Extrapolate using the element centres.
%       grid_vert_interp(Mobj, lon, lat, data, depth, mask, ...
%           'extrapolate', [Mobj.lonc, Mobj.latc])
46
%
47 48 49 50 51
% Author(s):
%   Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
%   2013-02-08 First version.
52 53
%   2013-05-16 Add support for parallel for-loops (not mandatory, but
%   enabled if the Parallel Computing Toolbox is available).
54 55 56 57 58 59 60 61
%   2014-04-28 Add new argument to allow specifying different coordinates
%   for the extrapolation. This allows us to interpolate data which belongs
%   either to the element centres or element nodes (defaults to element
%   nodes). This is only really important when the model grid falls outside
%   the coverage of the supplied data and we're extrapolating data. Update
%   the parallel pool code to use the new parpool function instead of
%   matlabpool in anticipation of the latter's eventual removal from
%   MATLAB. Also update the help.
62 63
%   2014-06-12 Fix bug in interpolating in the vertical when HYCOM data has
%   only a single depth bin.
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
%
%==========================================================================

subname = 'grid_vert_interp';

global ftbverbose
if ftbverbose
    fprintf('\n')
    fprintf(['begin : ' subname '\n'])
end

if ~isfield(Mobj, 'siglayz')
    error('Error: missing required sigma layer depth values for the unstructured grid.')
end

if ~Mobj.have_lonlat
    error('Need spherical coordinates')
end

if sum(size(lon) ~= size(data(:, :, 1))) ~= 0 || sum(size(lat) ~= size(data(:, :, 1))) ~= 0
    error('Size of the longitude or latitude arrays do not match the supplied data array')
end

87 88 89 90 91 92 93 94 95 96 97 98
% Extract the extrapolation coordinates. Default to nodes but use
% whatever's given if we have the 'extrapolate' argument.
ulon = Mobj.lon;
ulat = Mobj.lat;
for aa = 1:2:length(varargin)
    switch varargin{aa}
        case 'extrapolate'
            ulon = varargin{aa + 1}(:, 1);
            ulat = varargin{aa + 1}(:, 2);
    end
end

99 100 101
wasOpened = false;
if license('test', 'Distrib_Computing_Toolbox')
    % We have the Parallel Computing Toolbox, so launch a bunch of workers.
102 103
    try
        % New version for MATLAB 2014a (I think) onwards.
104
        if isempty(gcp('nocreate'))
105 106 107 108 109 110 111 112 113 114
            pool = parpool('local');
            wasOpened = true;
        end
    catch
        % Version for pre-2014a MATLAB.
        if matlabpool('size') == 0
            % Force pool to be local in case we have remote pools available.
            matlabpool open local
            wasOpened = true;
        end
115 116 117
    end
end

118 119 120 121 122 123
[~, fz] = size(Mobj.siglayz);
[nx, ny, ~, ~] = size(data);

% Preallocate the output arrays
dataz = nan(nx, ny, fz);

124 125 126 127 128
if ftbverbose
    tic
end

parfor xi = 1:nx
129
    % Get all the y and z dimension data for the current x position
130
    % (temperature, salinity and depth).
131 132 133
    xdata = squeeze(data(xi, :, :));
    xdepth = squeeze(depth(xi, :, :));
    xmask = mask(xi, :);
134

135
    % Preallocate the arrays for the inner loop.
136 137 138 139 140 141 142
    ydata = nan(ny, fz);
    for yi = 1:ny
        if xmask(yi)
            continue
        end

        % Find the nearest sigma layer z values from the unstructured grid.
143
        [md, mi] = min(sqrt((ulon - lon(xi, yi)).^2 + (ulat - lat(xi, yi)).^2));
144 145 146 147 148 149 150 151 152

        % Skip data point if the closest FVCOM node is more than 10 minutes
        % away.
        if md > 10 / 60
%             if ftbverbose
%                 fprintf('%s : skipping %f, %f (more than 10 minutes from the nearest unstructured grid node),\n', subname, lon(yi, xi), lat(yi, xi))
%             end
            continue
        else
153 154
            % Use the FVCOM node's sigma depths to interpolate this regular
            % grid position's temperature and salinity data.
155

156 157 158 159 160 161 162 163 164
            % Get the FVCOM depths closest to this regular grid position.
            try
                tfz = Mobj.siglayz(mi, :);
            catch
                tfz = Mobj.siglayzc(mi, :);
            end
            % Now get the regular grid depths at this node for all the
            % vertical layers and linearly interpolate through the water
            % column onto the FVCOM vertical profile.
165
            tpz = xdepth(yi, :);
166 167 168 169 170

            % Remove any NaN values in the vertical depths (as is the case
            % with the HYCOM data, and interpolate only the data we have
            % that are finite).
            nidx = isnan(tpz);
171 172 173 174 175 176 177 178
            % If we have only a single value, repeat it down the entire
            % water column, otherwise, interpolate linearly (with
            % extrapolation).
            if length(tpz(~nidx)) == 1
                ydata(yi, :) = tpz(~nidx);
            else
                ydata(yi, :) = interp1(tpz(~nidx), xdata(yi, ~nidx), tfz, 'linear', 'extrap');
            end
179 180 181 182
        end
    end
    dataz(xi, :, :) = ydata;
end
183 184 185

% Close the MATLAB pool if we opened it.
if wasOpened
186 187 188 189 190
    try
        pool.delete
    catch
        matlabpool close
    end
191 192 193 194 195
end

if ftbverbose
    toc
end
196 197 198 199

if ftbverbose
    fprintf(['end   : ' subname '\n'])
end