Skip to content

GitLab

fvcom-toolbox
fvcom-toolbox
Commit ca7e2f1f authored by Pierre Cazenave's avatar Pierre Cazenave
Browse files
Options

Move the code to write the NetCDF file into a new function... 

Move the code to write the NetCDF file into a new function (write_FVCOM_forcing.m). Also improve the speed of the interpolation in two ways: 1. Use TriScatteredInterp to do the interpolation (with natural neighbour interpolation), and 2. Do the interpolation onto the element centres using their positions directly, rather than averaging the nodal values for each element. This last step significantly improved the speed of the interpolation
parent a6f70532
Hide whitespace changes
Inline Side-by-side
Showing with 104 additions and 186 deletions
fvcom_prepro/grid2fvcom.m 0 → 100644
View file @ ca7e2f1f
function fvcom = grid2fvcom(Mobj,vars,wind)
% 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.).
% wind - a struct which contains the following arrays:
% 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
% element centres.
%
% 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
ntimes = numel(wind.time);
catch
ntimes = numel(wind.(vars{1}).time);
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)
if strcmpi(vars{vv}, 'lat') || strcmpi(vars{vv}, 'lon') || strcmpi(vars{vv}, 'x') || strcmpi(vars{vv}, 'y') || strcmpi(vars{vv}, 'time')
fprintf('skipping variable %s\n', vars{vv})
continue
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);
currvar = wind.(vars{1}).data(:, :, 1);
% 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)
ftsin = TriScatteredInterp(wind.x(:), wind.y(:), currvar(:), 'natural');
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
fvcom_prepro/grid2fvcom_U10V10.m deleted 100644 → 0
View file @ a6f70532
function [fvcom_u10_node, fvcom_v10_node] = grid2fvcom_U10V10(Mobj,wind,fvcom_forcing_file,infos)
% Interpolate regularly gridded wind speed data onto a given FVCOM grid
%
% grid2fvcom_U10V10(Mobj,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
% wind - a struct which contains the following arrays:
% x - x data (probably best in cartesian for the interpolation)
% y - y data (probably best in cartesian for the interpolation)
% u10 - u-component wind data
% v10 - v-component wind data
% time - time vector (in Modified Julian Days)
% fvcom_forcing_file - FVCOM forcing file name
% infos - Additional remarks to be written to the "infos" NetCDF variable
%
% OUTPUT:
% fvcom_u10_node - interpolated u vector values at the nodes
% fvcom_v10_node - interpolated v vector values at the nodes
% FVCOM wind speed forcing NetCDF file
%
% 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.
%
%==========================================================================
warning off
if nargin ~= 4
error('Incorrect number of arguments')
end
subname = 'ncep2fvcom_U10V10';
global ftbverbose;
if(ftbverbose);
fprintf('\n')
fprintf(['begin : ' subname '\n'])
end;
%--------------------------------------------------------------------------
% Get the relevant bits from the FVCOM mesh object
%--------------------------------------------------------------------------
tri = Mobj.tri;
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);
ntimes = numel(wind.time);
% Interpolate NCEP data to FVCOM mesh
fvcom_u10 = zeros(nElems,ntimes);
fvcom_v10 = zeros(nElems,ntimes);
fvcom_u10_node = zeros(nVerts,ntimes);
fvcom_v10_node = zeros(nVerts,ntimes);
for i=1:ntimes
fprintf('interpolating frame %d of %d\n', i, ntimes);
fvcom_u10_node(:,i) = griddata(wind.x,wind.y,wind.u10(:,:,i),x,y);
fvcom_v10_node(:,i) = griddata(wind.x,wind.y,wind.v10(:,:,i),x,y);
for j=1:nElems
fvcom_u10(j,i) = mean(fvcom_u10_node(tri(j,1:3)));
fvcom_v10(j,i) = mean(fvcom_v10_node(tri(j,1:3)));
end
end
fprintf('interpolation complete\n');
%--------------------------------------------------------------------------
% Dump header and data for netcdf FVCOM forcing file
%--------------------------------------------------------------------------
nc = netcdf.create(fvcom_forcing_file, 'clobber');
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'type','FVCOM U10/V10 Forcing File')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'source','fvcom grid (unstructured) surface forcing')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'references','http://fvcom.smast.umassd.edu, http://codfish.smast.umassd.edu')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'institution','Plymouth Marine Laboratory')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'institution','ncep2fvcom_U10V10.m')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'infos',infos)
% Dimensions
three_dimid=netcdf.defDim(nc,'three',3);
nele_dimid=netcdf.defDim(nc,'nele',nElems);
node_dimid=netcdf.defDim(nc,'node',nVerts);
time_dimid=netcdf.defDim(nc,'time',netcdf.getConstant('NC_UNLIMITED'));
% Time vars
time_varid=netcdf.defVar(nc,'time','NC_FLOAT',time_dimid);
netcdf.putAtt(nc,time_varid,'long_name','time');
netcdf.putAtt(nc,time_varid,'units','days since 1858-11-17 00:00:00');
netcdf.putAtt(nc,time_varid,'format','modified julian day (MJD)');
netcdf.putAtt(nc,time_varid,'time_zone','UTC');
itime_varid=netcdf.defVar(nc,'Itime','NC_INT',time_dimid);
netcdf.putAtt(nc,itime_varid,'units','days since 1858-11-17 00:00:00');
netcdf.putAtt(nc,itime_varid,'format','modified julian day (MJD)');
netcdf.putAtt(nc,itime_varid,'time_zone','UTC');
itime2_varid=netcdf.defVar(nc,'Itime2','NC_INT',time_dimid);
netcdf.putAtt(nc,itime2_varid,'units','msec since 00:00:00');
netcdf.putAtt(nc,itime2_varid,'time_zone','UTC');
% Space vars
x_varid=netcdf.defVar(nc,'x','NC_FLOAT',node_dimid);
netcdf.putAtt(nc,x_varid,'long_name','nodal x-coordinate');
netcdf.putAtt(nc,x_varid,'units','m');
y_varid=netcdf.defVar(nc,'y','NC_FLOAT',node_dimid);
netcdf.putAtt(nc,y_varid,'long_name','nodal y-coordinate');
netcdf.putAtt(nc,y_varid,'units','m');
nv_varid=netcdf.defVar(nc,'nv','NC_FLOAT',[nele_dimid, three_dimid]);
netcdf.putAtt(nc,nv_varid,'long_name','nodes surrounding element');
% Wind vars
u10_varid=netcdf.defVar(nc,'U10','NC_FLOAT',[nele_dimid, time_dimid]);
netcdf.putAtt(nc,u10_varid,'long_name','Eastward 10-m Velocity');
netcdf.putAtt(nc,u10_varid,'standard_name','Eastward Wind Speed');
netcdf.putAtt(nc,u10_varid,'units','m/s');
netcdf.putAtt(nc,u10_varid,'grid','fvcom_grid');
netcdf.putAtt(nc,u10_varid,'type','data');
v10_varid=netcdf.defVar(nc,'V10','NC_FLOAT',[nele_dimid, time_dimid]);
netcdf.putAtt(nc,v10_varid,'long_name','Northward 10-m Velocity');
netcdf.putAtt(nc,v10_varid,'standard_name','Northward Wind Speed');
netcdf.putAtt(nc,v10_varid,'units','m/s');
netcdf.putAtt(nc,v10_varid,'grid','fvcom_grid');
netcdf.putAtt(nc,v10_varid,'type','data');
u10_node_varid=netcdf.defVar(nc,'U10_node','NC_FLOAT',[node_dimid, time_dimid]);
netcdf.putAtt(nc,u10_node_varid,'long_name','Eastward 10-m Velocity');
netcdf.putAtt(nc,u10_node_varid,'standard_name','Eastward Wind Speed');
netcdf.putAtt(nc,u10_node_varid,'units','m/s');
netcdf.putAtt(nc,u10_node_varid,'grid','fvcom_grid');
netcdf.putAtt(nc,u10_node_varid,'type','data');
v10_node_varid=netcdf.defVar(nc,'V10_node','NC_FLOAT',[node_dimid, time_dimid]);
netcdf.putAtt(nc,v10_node_varid,'long_name','Northward 10-m Velocity');
netcdf.putAtt(nc,v10_node_varid,'standard_name','Northward Wind Speed');
netcdf.putAtt(nc,v10_node_varid,'units','m/s');
netcdf.putAtt(nc,v10_node_varid,'grid','fvcom_grid');
netcdf.putAtt(nc,v10_node_varid,'type','data');
% End definitions
netcdf.endDef(nc);
% Write data
netcdf.putVar(nc,nv_varid, tri');
netcdf.putVar(nc,time_varid,0,ntimes,wind.time);
netcdf.putVar(nc,itime_varid,0,ntimes,floor(wind.time));
netcdf.putVar(nc,itime2_varid,0,ntimes,mod(wind.time,1)*24*3600*1000);
netcdf.putVar(nc,x_varid,x);
netcdf.putVar(nc,y_varid,y);
netcdf.putVar(nc,u10_varid,[0,0],[nElems,ntimes],fvcom_u10);
netcdf.putVar(nc,v10_varid,[0,0],[nElems,ntimes],fvcom_v10);
netcdf.putVar(nc,u10_node_varid,[0,0],[nVerts,ntimes],fvcom_u10_node);
netcdf.putVar(nc,v10_node_varid,[0,0],[nVerts,ntimes],fvcom_v10_node);
% Close the NetCDF files
netcdf.close(nc);
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment