grid_vert_interp.m 4.48 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
function dataz = grid_vert_interp(Mobj, lon, lat, data, depth, mask)
% 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 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
%
% INPUT:
%   Mobj        = MATLAB mesh structure which must contain:
%                   - Mobj.siglayz - sigma layer depths for all model
%                   nodes.
%                   - Mobj.lon, Mobj.lat - node coordinates (long/lat).
%   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.
%   mesh        = logical array of positions outside the regularly gridded
%   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).
%   
% OUTPUT:
%   x by y by z array of vertically interpolated values at each regular
%   grid location.
%
% EXAMPLE USAGE
%   grid_vert_interp(Mobj, lon, lat, data, depth, mask)
% 
% Author(s):
%   Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
%   2013-02-08 First version.
42 43
%   2013-05-16 Add support for parallel for-loops (not mandatory, but
%   enabled if the Parallel Computing Toolbox is available).
44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
%
%==========================================================================

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

67 68 69 70 71 72 73 74 75 76
wasOpened = false;
if license('test', 'Distrib_Computing_Toolbox')
    % We have the Parallel Computing Toolbox, so launch a bunch of workers.
    if matlabpool('size') == 0
        % Force pool to be local in case we have remote pools available.
        matlabpool open local
        wasOpened = true;
    end
end

77 78 79 80 81 82
[~, fz] = size(Mobj.siglayz);
[nx, ny, ~, ~] = size(data);

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

83 84 85 86 87
if ftbverbose
    tic
end

parfor xi = 1:nx
88
    % Get all the y and z dimension data for the current x position
89
    % (temperature, salinity and depth).
90 91 92 93
    xdata = squeeze(data(xi, :, :));
    xdepth = squeeze(depth(xi, :, :));
    xmask = mask(xi, :);
    
94
    % Preallocate the arrays for the inner loop.
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
    ydata = nan(ny, fz);
    for yi = 1:ny
        if xmask(yi)
            continue
        end

        % Find the nearest sigma layer z values from the unstructured grid.
        [md, mi] = min(sqrt((Mobj.lon - lon(xi, yi)).^2 + (Mobj.lat - lat(xi, yi)).^2));

        % 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
            % Use the FVCOM node's sigma depths to interpolate this POLCOMS
            % position's temperature and salinity data.
            
            % Get the FVCOM depths closest to this POLCOMS grid position.
            tfz = Mobj.siglayz(mi, :);
117 118 119
            % Now get the POLCOMS depths at this node for all the vertical
            % layers and linearly interpolate through the water column onto
            % the FVCOM vertical profile.
120 121 122 123 124 125
            tpz = xdepth(yi, :);
            ydata(yi, :) = interp1(tpz, xdata(yi, :), tfz, 'linear', 'extrap');
        end
    end
    dataz(xi, :, :) = ydata;
end
126 127 128 129 130 131 132 133 134

% Close the MATLAB pool if we opened it.
if wasOpened
    matlabpool close
end

if ftbverbose
    toc
end
135 136 137 138

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