grid2fvcom.m 5.84 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
%       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.
19 20 21 22
%       time - time vector (in Modified Julian Days). If you're using some
%       of the NCEP surface products (e.g. relative humitidy, sea level
%       pressure), you need to supply x and y coordinates for their grids
%       as .xalt and .yalt).
23
%
24 25
% OUTPUT:
%   fvcom - struct of the interpolated data values at the model nodes and
26 27
%   element centres. Also includes any variables which were in the input
%   struct but which have not been interpolated (e.g. time).
28
%
29
% NOTE:
30 31 32 33 34
%   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.
35
%
36 37
% Author(s):
%   Pierre Cazenave (Plymouth Marine Laboratory)
38
%
39 40 41 42 43 44
% 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:
45 46
%   it's easier to subsample and provide the subsampled data here).
%   2012-10-17 Add outputs to the function for use in visualisation.
47 48 49
%   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
50
%   write_FVCOM_forcing.m and made this purely an interpolation script.
51 52 53 54
%   2013-02-14 Add support for interpolating data on two different grids
%   through the .xalt and .yalt fields in the input data structure. This is
%   handy if you've got data from NCEP from both the Surface and Gaussian
%   products, each of which are on different grids.
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
%==========================================================================

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
87
    ntimes = numel(data.time);
88
catch
89
    ntimes = numel(data.(vars{1}).time);
90 91 92 93 94 95
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)
96 97
    if strcmpi(vars{vv}, 'time')
        fprintf('transferring variable %s as is\n', vars{vv})
98
        fvcom.(vars{vv}) = data.(vars{vv});
99
        continue
100 101
    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})
102
        fvcom.(vars{vv}) = Mobj.(vars{vv});
103 104 105 106 107 108 109
    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);
110
            currvar = data.(vars{vv}).data(:, :, i);
111 112 113 114
            % 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)
115 116 117 118 119 120 121 122 123 124 125
            % Use a try/catch to account for the different grids over which
            % the humidity and sealevel pressure data are sampled.
            try
                ftsin = TriScatteredInterp(data.x(:), data.y(:), currvar(:), 'natural');
            catch err
                warning([err.identifier, ': Some NCEP data are projected' ...
                    ' onto a different grid. Check you have specified' ...
                    ' data.xalt and data.yalt arrays which are on the' ...
                    ' same grid as the data to be interpolated.'])
                ftsin = TriScatteredInterp(data.xalt(:), data.yalt(:), currvar(:), 'natural');
            end
126 127
            fvcom.(vars{vv}).node(:,i) = ftsin(x,y);
            fvcom.(vars{vv}).data(:,i) = ftsin(xc,yc);
128 129 130 131 132 133 134
            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))
135
            end
136 137 138 139 140
        end
        fprintf('interpolation of %s complete\n', vars{vv});
    end
end

141 142
if ftbverbose;
    fprintf(['end   : ' subname '\n'])
143
end