From ca7e2f1fa11bc6ef76bc5a5653294a49dbf5c445 Mon Sep 17 00:00:00 2001 From: Pierre Cazenave Date: Fri, 2 Nov 2012 14:39:33 +0000 Subject: [PATCH] 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 --- fvcom_prepro/grid2fvcom.m | 104 +++++++++++++++++ fvcom_prepro/grid2fvcom_U10V10.m | 186 ------------------------------- 2 files changed, 104 insertions(+), 186 deletions(-) create mode 100644 fvcom_prepro/grid2fvcom.m delete mode 100644 fvcom_prepro/grid2fvcom_U10V10.m diff --git a/fvcom_prepro/grid2fvcom.m b/fvcom_prepro/grid2fvcom.m new file mode 100644 index 0000000..5533e26 --- /dev/null +++ b/fvcom_prepro/grid2fvcom.m @@ -0,0 +1,104 @@ +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 + diff --git a/fvcom_prepro/grid2fvcom_U10V10.m b/fvcom_prepro/grid2fvcom_U10V10.m deleted file mode 100644 index d11b839..0000000 --- a/fvcom_prepro/grid2fvcom_U10V10.m +++ /dev/null @@ -1,186 +0,0 @@ -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); - -- GitLab