grid2fvcom.m 4.81 KB
Newer Older
1
function fvcom = grid2fvcom(Mobj,vars,data)
2 3
% Interpolate regularly gridded surface forcing data onto a given FVCOM
% grid.
4
%
5
% grid2fvcom(Mobj,vars,data)
6
%
7 8
% DESCRIPTION:
%   Takes a given NCEP reanalysis grid file and interpolates the U10 and
9 10
%   V10 values onto the specified FVCOM grid file.
%
11 12 13 14
% INPUT:
%   Mobj - MATLAB mesh object
%   vars - a cell array of the variables to be interpolated on the FVCOM
%   grid in Mobj (e.g. uwnd, U10, vwnd, V10 etc.).
15
%   data - a struct which contains the following arrays:
16 17 18 19
%       x - x data (probably best in cartesian for the interpolation)
%       y - y data (probably best in cartesian for the interpolation)
%       The struct must also contain all the variables defined in vars.
%       time - time vector (in Modified Julian Days)
20
%
21 22
% OUTPUT:
%   fvcom - struct of the interpolated data values at the model nodes and
23 24
%   element centres. Also includes any variables which were in the input
%   struct but which have not been interpolated (e.g. time).
25
%
26
% NOTE:
27 28 29 30 31
%   The shape of the returned arrays for rhum and slp (via
%   get_NCEP_forcing.m) have sometimes differed from the other vairables
%   (they appear to be projected onto a different grid). Unless you
%   desperately need them, I would suggest omitting them from the
%   interpolation here as this assumes the arrays are all the same size.
32
%
33 34
% Author(s):
%   Pierre Cazenave (Plymouth Marine Laboratory)
35
%
36 37 38 39 40 41
% Revision history:
%   2012-10-15 First version based on ncep2fvcom_U10V10.m in the
%   fvcom-toolbox.
%   2012-10-16 Removed the code to read the NCEP file. Instead, farmed that
%   out to a new function (read_NCEP_wind) so that the relevant section can
%   be more readily extracted (rather than using the entire globe's data:
42 43
%   it's easier to subsample and provide the subsampled data here).
%   2012-10-17 Add outputs to the function for use in visualisation.
44 45 46
%   2012-10-19 Add wind struct as input rather than separate u, v, time and
%   lat/long arrays. Makes invocation a bit cleaner.
%   2012-11-01 Farmed out the creation of the NetCDF file to
47 48
%   write_FVCOM_forcing.m and made this purely an interpolation script.
%
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
%==========================================================================

if nargin ~= 3
    error('Incorrect number of arguments')
end

subname = 'grid2fvcom';

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

%--------------------------------------------------------------------------
% Get the relevant bits from the FVCOM mesh object
%--------------------------------------------------------------------------
x   = Mobj.x;
y   = Mobj.y;
nVerts = Mobj.nVerts;
nElems = Mobj.nElems;
if(ftbverbose);
    fprintf('info for FVCOM domain\n');
    fprintf('number of nodes: %d\n',nVerts);
    fprintf('number of elems: %d\n',nElems);
end

xc = nodes2elems(x, Mobj);
yc = nodes2elems(y, Mobj);

try
80
    ntimes = numel(data.time);
81
catch
82
    ntimes = numel(data.(vars{1}).time);
83 84 85 86 87 88
end

% Interpolate supplied regularly gridded data to FVCOM mesh. Use
% TriScatteredInterp to do the interpolation instead of griddata (should be
% faster).
for vv=1:length(vars)
89 90
    if strcmpi(vars{vv}, 'time')
        fprintf('transferring variable %s as is\n', vars{vv})
91
        fvcom.(vars{vv}) = data.(vars{vv});
92
        continue
93 94
    elseif strcmpi(vars{vv}, 'lat') || strcmpi(vars{vv}, 'lon') || strcmpi(vars{vv}, 'x') || strcmpi(vars{vv}, 'y')
        fprintf('reassigning variable %s from unstructured grid\n', vars{vv})
95
        fvcom.(vars{vv}) = Mobj.(vars{vv});
96 97 98 99 100 101 102
    else
        % Preallocate the output arrays
        fvcom.(vars{vv}).data = zeros(nElems,ntimes);
        fvcom.(vars{vv}).node = zeros(nVerts,ntimes);

        for i=1:ntimes
            fprintf('interpolating %s, frame %d of %d\n', vars{vv}, i, ntimes);
103
            currvar = data.(vars{vv}).data(:, :, i);
104 105 106 107
            % griddata way (cubic interpolation)
            %fvcom.(vars{vv}).node(:,i) = griddata(wind.x,wind.y,currvar,x,y,'cubic');
            %fvcom.(vars{vv}).data(:,i) = griddata(wind.x,wind.y,currvar,xc,yc,'cubic');
            % TriScatteredInterp way (with natural neighbour interpolation)
108
            ftsin = TriScatteredInterp(data.x(:), data.y(:), currvar(:), 'natural');
109 110
            fvcom.(vars{vv}).node(:,i) = ftsin(x,y);
            fvcom.(vars{vv}).data(:,i) = ftsin(xc,yc);
111 112 113 114 115 116 117
            nnans(1) = sum(isnan(fvcom.(vars{vv}).node(:,i)));
            nnans(2) = sum(isnan(fvcom.(vars{vv}).data(:,i)));
            if  nnans(1) > 0
                warning('%i NaNs in the interpolated node data. This won''t work with FVCOM.', nnans(1))
            end
            if nnans(2) > 0
                warning('%i NaNs in the interpolated element data. This won''t work with FVCOM.', nnans(2))
118
            end
119 120 121 122 123
        end
        fprintf('interpolation of %s complete\n', vars{vv});
    end
end

124 125 126
if ftbverbose;
    fprintf(['end   : ' subname '\n'])
end