Interpolate nc file from grid1 to grid2

ncfile_input is the file being interpolated, so must contain all the data.
ncfile_output is the target output file, only need to contain several time steps data.
You can use ncks -F -d time,1,10 input.nc example_fvcom_outfile_for_leem.nc to generate

%cd /Users/chenfuh/Documents/cleveland/mesh_files/truncated_river/closed_westbreakwall/generate_restart
%cd ~/Documents/LEEM_project/program_interp_fvcom2fvcom/

close all;clear all;clc;

%%
tic
%%%% start
clear ncfile_input ncfile_output
ncfile_input=['erie_avg_0001.nc']
ncfile_output=['example_fvcom_outfile_for_leem.nc']

%interp parameter
clear radius minum npow nquad
radius=0.10;
minum=0.0001;
npow=2;
nquad=2;
%
clear lon_* lat_* lonc_* latc_* siglay_* siglev_*
lon_old=double(ncread([ncfile_input],'lon'));
lat_old=double(ncread([ncfile_input],'lat'));
lonc_old=double(ncread([ncfile_input],'lonc'));
latc_old=double(ncread([ncfile_input],'latc'));
siglay_old=double(ncread([ncfile_input],'siglay'));
siglev_old=double(ncread([ncfile_input],'siglev'));

lon_new=double(ncread([ncfile_output],'lon'));
lat_new=double(ncread([ncfile_output],'lat'));
lonc_new=double(ncread([ncfile_output],'lonc'));
latc_new=double(ncread([ncfile_output],'latc'));
siglay_new=double(ncread([ncfile_output],'siglay'));
siglev_new=double(ncread([ncfile_output],'siglev'));

%check lon/lat lonc/latc to be consistent, either 0to360 or -180to180
if any(lon_old>180);lon_old=lon_old-360;end;
if any(lonc_old>180);lonc_old=lonc_old-360;end;
if any(lon_new>180);lon_new=lon_new-360;end;
if any(lonc_new>180);lonc_new=lonc_new-360;end;


clear nz_*
nz_lay_old=size(siglay_old,2);
nz_lev_old=size(siglev_old,2);

nz_lay_new=size(siglay_new,2);
nz_lev_new=size(siglev_new,2);
%
if (nz_lay_old~=nz_lay_new || nz_lev_old~=nz_lev_new)
    disp('the vertical dimension of new restart file DOES NOT consistent')
    disp('with the old restart file, this program is not able to solve such')
    disp('case, stopped')
end
    
%
clear aa bb nnvar nnvar2
aa=ncinfo([ncfile_input])
bb=ncinfo([ncfile_output])
nnvar=length(aa.Variables)
nnvar2=length(bb.Variables)
for i=1:nnvar
    cname{i}=aa.Variables(i).Name;
end
    cname


for i=1:nnvar
    [num2str(i) ' -- ' cname{i}]
    lonlat_coord=[strcmp(cname{i},'lon') strcmp(cname{i},'lonc') ...
                  strcmp(cname{i},'lat') strcmp(cname{i},'latc')];
    xy_coord=[strcmp(cname{i},'x') strcmp(cname{i},'xc') ...
              strcmp(cname{i},'y') strcmp(cname{i},'yc')];
          
    sigma_coord=[strcmp(cname{i},'siglay') strcmp(cname{i},'siglev')];
          
    depth=strcmp(cname{i},'h');
    
    control_area=[strcmp(cname{i},'art1') strcmp(cname{i},'art2')];
    
    mesh_grid_info=[strcmp(cname{i},'nv') strcmp(cname{i},'nbsn') ...
                       strcmp(cname{i},'nbve') strcmp(cname{i},'ntsn') ...
                       strcmp(cname{i},'ntve') strcmp(cname{i},'ntsn') ...
                       strcmp(cname{i},'a1u') strcmp(cname{i},'a2u') ...
                       strcmp(cname{i},'aw0') strcmp(cname{i},'awx') ...
                       strcmp(cname{i},'awy') strcmp(cname{i},'nbe')];
                   
    no_needed_var=[strcmp(cname{i},'latent_heat_flux') strcmp(cname{i},'sensible_heat_flux') ...
                   strcmp(cname{i},'long_wave') ...
                   strcmp(cname{i},'uwind_speed') strcmp(cname{i},'vwind_speed') ...
                   strcmp(cname{i},'uwind_stress') strcmp(cname{i},'vwind_stress')];             
    
    var=ncread([ncfile_input],cname{i});
    nndim=length(aa.Variables(i).Dimensions);
    %
    if (nndim==0)
        % doing nothing, keep unchanged in ncfile_output
        disp(['variable: ' cname{i} ' -- doing nothing, keep unchanged in ncfile_output']);
    elseif (nndim>0)
        dimname=aa.Variables(i).Dimensions(1).Name;
        stat1=strcmp(dimname,'time');
        stat2=strcmp(dimname,'node');
        stat3=strcmp(dimname,'nele');
        stat4=strcmp(dimname,'DateStrLen');
    end
    %
    [num2str(i) ' -- ' cname{i}]
    if (~any(lonlat_coord) & ~any(xy_coord) & ~depth & nndim>0 & ~any(no_needed_var) & ...
        ~any(mesh_grid_info) & ~any(control_area) & ~any(sigma_coord))
        if (stat1) % time variables
            clear var_new
            if strcmp(cname{i},'time') || strcmp(cname{i},'Itime')
                clear tdiff
                tdiff=[datenum(1970,1,1)-datenum(1858,11,17)]; %diff in reference time
                var_new=var+tdiff;
            else
                var_new=var;
            end
            ncwrite([ncfile_output],cname{i},var_new)
            
        elseif (stat2)
            clear var_new
            if (nndim==1) % grid info variables
                disp(['grid variable: ' cname{i} ' -- doing nothing, keep unchanged in ncfile_output']);
            elseif (nndim==2)
                dim2nd=aa.Variables(i).Dimensions(2).Name;
                if strcmp(dim2nd,'siglay') 
                    % should only be siglay, but it has already been excluded
                else
                    % dim2nd must be "time"
                    clear nt
                    nt=size(var,2);
                    for t=1:nt
                        clear tmp
                        tmp=interp_scatter_quick(lon_old,lat_old,double(var(:,t)),lon_new,lat_new,radius,minum,npow,nquad);
                        var_new(:,t)=tmp;
                    end
                end
                    ncwrite([ncfile_output],cname{i},var_new);
                    
            elseif (nndim==3)
              dim2nd=aa.Variables(i).Dimensions(2).Name;
              if strcmp(dim2nd,'siglay') % 3d var(node,siglay,time): silinity,temp
                  clear nt
                  nt=size(var,3);
                for t=1:nt
                for k=1:nz_lay_old
                    clear tmp
                    tmp=interp_scatter_quick(lon_old,lat_old,double(squeeze(var(:,k,t))),lon_new,lat_new,radius,minum,npow,nquad);
                    var_new(:,k,t)=tmp;
                end
                end
                    ncwrite([ncfile_output],cname{i},var_new)
                    
              elseif strcmp(dim2nd,'siglev') %3d var(node,siglev,time): omega
                  clear nt
                  nt=size(var,3);
                for t=1:nt  
                for k=1:nz_lev_old
                    tmp=interp_scatter_quick(lon_old,lat_old,double(squeeze(var(:,k,t))),lon_new,lat_new,radius,minum,npow,nquad);
                    var_new(:,k,t)=tmp;
                end  
                end
                    ncwrite([ncfile_output],cname{i},var_new)
              end
            end
                
        elseif (stat3)
            clear var_new
            if (nndim==1) %1d var: partition
                disp(['variable: ' cname{i} ' -- doing nothing, keep unchanged in ncfile_output']);

            elseif (nndim==2) % 2d var(nele,time): tauc,ua,va, (and not used) uwind_speed/stress,vwind_speed/stress
                clear nt
                nt=size(var,2);
                for t=1:nt
                    clear tmp
                    tmp=interp_scatter_quick(lonc_old,latc_old,double(var(:,t)),lonc_new,latc_new,radius,minum,npow,nquad);
                    var_new(:,t)=tmp;
                end
                ncwrite([ncfile_output],cname{i},var_new)
                
            elseif (nndim==3)
              dim2nd=aa.Variables(i).Dimensions(2).Name
              if strcmp(dim2nd,'siglay') %3d(nele,siglay,time): u,v,ww
                  clear nt
                  nt=size(var,3)
                for t=1:nt  
                for k=1:nz_lay_old
                    clear tmp
                    tmp=interp_scatter_quick(lonc_old,latc_old,double(squeeze(var(:,k,t))),lonc_new,latc_new,radius,minum,npow,nquad);
                    var_new(:,k,t)=tmp;  
                end
                end
                    ncwrite([ncfile_output],cname{i},var_new)
                    
              elseif strcmp(dim2nd,'siglev') %3d(nele,siglev,time),should be none
                  clear nt
                  nt=size(var,3)
                for t=1:nt  
                for k=1:nz_lev_old
                    clear tmp
                    tmp=interp_scatter_quick(lonc_old,latc_old,double(squeeze(var(:,k,t))),lonc_new,latc_new,radius,minum,npow,nquad);
                    var_new(:,k,t)=tmp;
                end  
                end
                    ncwrite([ncfile_output],cname{i},var_new)
              end
            end
        elseif (stat4) % var(DateStrLen,time): Times,file_date
                clear var_new
                var_new=var;
                ncwrite([ncfile_output],cname{i},var_new)
        end
    end
    disp(['finished interpolation: ' cname{i}]);
end 
disp(['finished all interpolation from ' ncfile_input ' to ' ncfile_output]);
%
toc

%% check

ncload('example_fvcom_outfile_for_gleDA.nc')
x1=lon;
y1=lat;
xc1=lonc;
yc1=latc;
cc1=u;

ncload('example_fvcom_outfile_for_leem.nc')
x2=lon;
y2=lat;
xc2=lonc;
yc2=latc;
cc2=u;

figure(1);clf;hold on;
scattercontourf(xc1,yc1,cc1(1,:)',[-1e-0:1e-1:1e-0])

figure(2);clf;hold on;
scattercontourf(xc2,yc2,cc2(1,:)',[-1e-0:1e-1:1e-0])


k=12
t=6

figure(1);clf;hold on;
scattercontourf(lon_old,lat_old,double(squeeze(var(:,k,t))),[0:1:30]); caxis([0 30]);

figure(2);clf;hold on;
scattercontourf(lon_new,lat_new,double(squeeze(var_new(:,k,1))),[0:1:30]); caxis([0 30]);

%
t=1
cmin=-0.1;cinc=0.01;cmax=0.1;
figure(1);clf;hold on;
scattercontourf(lonc_old,latc_old,double(squeeze(var(:,t))),[cmin:cinc:cmax]); caxis([cmin cmax]);

figure(2);clf;hold on;
scattercontourf(lonc_new,latc_new,double(squeeze(var_new(:,t))),[cmin:cinc:cmax]); caxis([cmin cmax]);

%
k=1
t=1
cmin=-0.1;cinc=0.01;cmax=0.1;
figure(1);clf;hold on;
scattercontourf(lonc_old,latc_old,double(squeeze(var(:,k,t))),[cmin:cinc:cmax]); caxis([cmin cmax]);

figure(2);clf;hold on;
scattercontourf(lonc_new,latc_new,double(squeeze(var_new(:,k,t))),[cmin:cinc:cmax]); caxis([cmin cmax]);

Notice: The “time” and “Itime” in the interpolated file is wrong, and need to be modified. Please use the following code to modify.

clc;clear all;close all;

clear ncfile_input ncfile_output
ncfile_input=['erie_avg_0001.nc']
ncfile_output=['example_fvcom_outfile_for_leem.nc']


time_old=double(ncread([ncfile_input],'time'));
time_new=double(ncread([ncfile_output],'time'));
ncwrite(ncfile_output,'time',time_old);

Itime_old=double(ncread([ncfile_input],'Itime'));
Itime_new=double(ncread([ncfile_output],'Itime'));
ncwrite(ncfile_output,'Itime',Itime_old);

Run matlab script by linux command

matlab -batch "script_name"