source: altifloat/matlab_toolbox/interp_wind.m @ 199

function [] = interp_wind (windfile,meshfile,outfile)
%[] = interp_wind (windfile,meshfile,outfile)
%outfile is a netcdf file


%Test the function
%windfile = '../../data/wind/wind_20130825-20130930.nc';
%meshfile = '../exp_dan_lr/meshgrid_aviso.dat';
%outfile  = '../exp_dan_lr/wind_interp.nc';

ufield = 'u10';
vfield = 'v10';

fields = {ufield,vfield};

copyfile (windfile,outfile);

%mehsgrid
mg=load(meshfile);
Nlat=length(unique(mg(:,2))); %Phi
Nlong=length(unique(mg(:,1))); %Lambda
GR_long=reshape(mg(:,1),Nlong,Nlat); 
GR_lat=reshape(mg(:,2),Nlong,Nlat); 

%Grid point size
R_earth=6371229; %in meters
delta{1}=zeros(Nlong,Nlat); %deltax
delta{2}=zeros(Nlong,Nlat); %deltay (should correspond to fields order)

for ilon=1:Nlong-1;
    for jlat=1:Nlat-1;
        
    delta{1}(ilon,jlat)=R_earth*(2*pi/360)*(GR_long(ilon+1,jlat)-GR_long(ilon,jlat))...
        .*cos(2*pi*GR_lat(ilon,jlat)/360); 
    delta{2}(ilon,jlat)=R_earth*(2*pi/360)*(GR_lat(ilon,jlat+1)-GR_lat(ilon,jlat));  
    
    
end 
end



delta{1}(Nlong,1:Nlat)=delta{1}(Nlong-1,1:Nlat);
delta{2}(Nlong,1:Nlat)=delta{2}(Nlong-1,1:Nlat);


delta{1}(1:Nlong,Nlat)=delta{1}(1:Nlong,Nlat-1);
delta{2}(1:Nlong,Nlat)=delta{2}(1:Nlong,Nlat-1);


%Open netcdf
ncid = netcdf.open(outfile,'WRITE');

%lon/lat

varid = netcdf.inqVarID(ncid,'latitude');
lat = netcdf.getVar(ncid,varid);
varid = netcdf.inqVarID(ncid,'longitude');
lon = netcdf.getVar(ncid,varid);


[LAT,LON]=meshgrid(lat,lon); %C'est à l'envers à cause de matlab

timeid = netcdf.inqDimID(ncid,'time');

%Make new dimensions 
netcdf.reDef(ncid);
latdimid = netcdf.defDim(ncid,'lat_interp',Nlat);
londimid = netcdf.defDim(ncid,'lon_interp',Nlong);

latid = netcdf.defVar(ncid,'lat_interp','float',latdimid);
lonid = netcdf.defVar(ncid,'lon_interp','float',londimid);

netcdf.endDef(ncid);

%Make new lat/lon variables in netcdf
netcdf.putVar(ncid,latid,unique(mg(:,2)));
netcdf.putVar(ncid,lonid,unique(mg(:,1)));

for k=1:length(fields)
%for k=1:1
  
  %Read var and attributes
  varid = netcdf.inqVarID(ncid,fields{k});
  data = netcdf.getVar(ncid,varid); %time is the last dimension
  scale_factor = netcdf.getAtt(ncid,varid,'scale_factor');
  add_offset = netcdf.getAtt(ncid,varid,'add_offset');
  FillValue = netcdf.getAtt(ncid,varid,'_FillValue');
  missing_value = netcdf.getAtt(ncid,varid,'missing_value');
  %units = netcdf.getAtt(ncid,varid,'units');
  %long_name = netcdf.getAtt(ncid,varid,'long_name');
  
  %Convert to floating point
  data_f = double(data)*scale_factor + add_offset ;
  data_f(data == FillValue | data == missing_value) = nan ;
  
  %interp on meshgrid
  data_fi = nan*ones(Nlong,Nlat,size(data,3));
  for it=1:size(data,3)
  data_fi(:,:,it) = interp2(LAT,LON,data_f(:,:,it),GR_lat,GR_long,'nearest')./delta{k};
  end
  
  % Save in the ncfile
  netcdf.reDef(ncid);
  dimids = [londimid latdimid timeid];
  newid = netcdf.defVar(ncid,[fields{k} '_f'],'float',dimids);
  %netcdf.copyAtt(ncid,varid,'units',ncid,newid);
  netcdf.putAtt(ncid,newid,'units','gridpoint per second');
  netcdf.copyAtt(ncid,varid,'long_name',ncid,newid);
  netcdf.putAtt(ncid,newid,'Description',['Convert to floating point and interpolated using' ...
                      meshfile ' file - ' datestr(now)]);
  netcdf.endDef(ncid);
  netcdf.putVar(ncid,newid,data_fi);
  
  
end %for k



netcdf.close(ncid);