grid2fvcom.m 4.37 KB
Newer Older
1
function fvcom = grid2fvcom(Mobj,vars,data)
2 3 4 5 6 7 8 9 10 11 12 13
% Interpolate regularly gridded wind speed data onto a given FVCOM grid
%
% grid2fvcom(Mobj,vars,wind,fvcom_forcing_file,infos)
% 
% DESCRIPTION:
%   Takes a given NCEP reanalysis grid file and interpolates the U10 and
%   V10 values onto the specified FVCOM grid file. 
%   
% 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.).
14
%   data - a struct which contains the following arrays:
15 16 17 18 19 20 21
%       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)
% 
% OUTPUT:
%   fvcom - struct of the interpolated data values at the model nodes and
22 23
%   element centres. Also includes any variables which were in the input
%   struct but which have not been interpolated (e.g. time).
24 25 26 27 28 29 30
%
% NOTE: 
%   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.
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 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 80
% 
% Author(s):
%   Pierre Cazenave (Plymouth Marine Laboratory)
% 
% 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:
%   it's easier to subsample and provide the subsampled data here). 
%   2012-10-17 Add outputs to the function for use in visualisation. 
%   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
%   write_FVCOM_forcing.m and made this purely an interpolation script. 
% 
%==========================================================================

warning off

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
81
    ntimes = numel(data.time);
82
catch
83
    ntimes = numel(data.(vars{1}).time);
84 85 86 87 88 89
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)
90 91 92
    if strcmpi(vars{vv}, 'time')
        fprintf('transferring variable %s as is\n', vars{vv})
        fvcom.(vars{vv}) = data.(vars{1});
93
        continue
94 95 96
    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})
        fvcom.(vars{vv}) = Mobj.(vars{1});
97 98 99 100 101 102 103
    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);
104
            currvar = data.(vars{vv}).data(:, :, i);
105 106 107 108
            % 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)
109
            ftsin = TriScatteredInterp(data.x(:), data.y(:), currvar(:), 'natural');
110 111 112 113 114 115 116
            fvcom.(vars{vv}).node(:,i) = ftsin(x,y);
            fvcom.(vars{vv}).data(:,i) = ftsin(xc,yc);
        end
        fprintf('interpolation of %s complete\n', vars{vv});
    end
end