clear
clc
%% read all the .nc files in the origin folder

% dirName = 'd:\openearthrawdata\trunk\rijkswaterstaat\vaklodingen\interpolation\products\nc\rijkswaterstaat\vaklodingen';  % folder path
dirName = 'c:\Users\santinel\products\nc\rijkswaterstaat\vaklodingen';  % folder path
newdir = 'c:\Users\santinel\products_interp\nc\rijkswaterstaat\vaklodingen_interp'; % new path

beg_year = 1926;
end_year = 2016;

fl = dir( fullfile(dirName,'*.nc') );   % list all files
files = {fl.name}';                      % file names
files = cell2mat(files(1:end));
% copy them in a different folder
mkdir(newdir);
for ii=1:length(files)
    copyfile([dirName filesep files(ii,:)],[newdir filesep files(ii,:)]);
end
for i =1:length(fl);
    disp([num2str(i) '/' num2str(length(fl))])
    Zi = nc_varget([newdir filesep files(i,:)],'z');
    lon = nc_varget([newdir filesep files(i,:)],'lon');
    lat = nc_varget([newdir filesep files(i,:)],'lat');
    date = nc_cf_time([newdir filesep files(i,:)],'time');
    oldfile = copyfile([newdir filesep files(i,:)],newdir);
    % create new variables to be used
    timesr = NaN(1,length(beg_year:2016));
    Zf = NaN(length(timesr),length(lon(:,1)),length(lon(1,:)));
    
    for j=1:length(beg_year:2016);
        timesr(j)=datenum(['01-01-',num2str(beg_year-1+j)]);
    end
    
    % FIND INDEX OF ~NaN VALUES INSIDE ALTITUDE DATA
    [r,c,v] = ind2sub(size(Zi),find(~isnan(Zi)));
    t = find(ismember(timesr,date));
    t = t(r)';
    
    % ASSIGN DATA AVAILABLE TO THE NEW ALTITUDE VARIABLE
    for x=1:length(r);
        Zf(t(x),c(x),v(x)) = Zi(r(x),c(x),v(x));
    end
    
    % INTERPOLATE TO OBTAIN NEW DATA
    Zfmod = NaN(length(timesr)+1, size(Zf,2), size(Zf,3));
    ZfN = logical(zeros(length(timesr)+1, size(Zf,2), size(Zf,3)));
    
    ZfN(2:end,:,:) = ~isnan(Zf);
    ZfN(1,:,:)=logical(ones(1,size(Zf,2), size(Zf,3)));
    Zfmod(2:end,:,:) = Zf;
    Zfmod(isnan(Zfmod)) = interp1(find(ZfN), Zfmod(ZfN), find(isnan(Zfmod)));
    Zf1= Zfmod(2:end,:,:);
 
    % ASSIGN THE LAST VALUE TO THE LAST NaNs
    Zf2 = Zf1;
    for a =1:length(c);
        ix = find(~isnan(Zf2(:,c(a),v(a))));
        Zf2(ix(end):length(timesr),c(a),v(a))= Zf2(ix(end),c(a),v(a));
    end
    %% change content of nc inside copy folder [keeping the same name of nc files, different folder]
    % global attributes
    if ispc
        userdir = getenv('USERPROFILE');
    else
        userdir = getenv('HOME');
    end
    publisher_email    = 'Giorgio.Santinelli@deltares.nl';
    publisher_name     = getenv('USERNAME');
    publisher_domain   = getenv('USERDOMAIN');
    publisher_computer = getenv('COMPUTERNAME');
    publish_script_id  = '$Header: https://svn.oss.deltares.nl/repos/openearthrawdata/trunk/rijkswaterstaat/vaklodingen/scripts/ncgen_interpolated/ncgen_vaklodingen.m 8899 2016-04-16 07:58:24Z heijer $';
    fmt = 'yyyy-mm-ddTHH:MMZ';
    ncwriteatt([newdir filesep files(i,:)],'/','id',['vaklodingen_interpolated_release', datestr(nowutc,'yyyymmdd')]);
    ncwriteatt([newdir filesep files(i,:)],'/','metadata_link','https://publicwiki.deltares.nl/display/OET/Vaklodingen+Interpolated+Dataset')
    ncwriteatt([newdir filesep files(i,:)],'/','title','Vaklodingen_interpolated')
    ncwriteatt([newdir filesep files(i,:)],'/','summary','Linear interpolation in time of bathymetry and topography measurements along the Dutch coast')
    ncwriteatt([newdir filesep files(i,:)],'/','history',sprintf('Created on %s by %s on computer %s\\%s with script %s',datestr(nowutc, fmt), publisher_name, publisher_domain, publisher_computer, publish_script_id))
    ncwriteatt([newdir filesep files(i,:)],'/','date_issued', datestr(nowutc, fmt))
	ncwriteatt([newdir filesep files(i,:)],'/','publisher_name', publisher_name)
    ncwriteatt([newdir filesep files(i,:)],'/','publisher_email', publisher_email)
    ncwriteatt([newdir filesep files(i,:)],'/','date_created', datestr(nowutc, fmt))
    ncwriteatt([newdir filesep files(i,:)],'/','date_modified', datestr(nowutc, fmt))
    ncwriteatt([newdir filesep files(i,:)],'/','timeCoverage', [num2str(beg_year), '-01-01T00:00:00 - ', num2str(end_year), '-01-01T00:00:00'])
    ncwriteatt([newdir filesep files(i,:)],'/','flag_years', num2str(find(ismember(timesr,date)))) % add years
    % variables
    nc_varput([newdir filesep files(i,:)],'time',timesr - datenum('01-01-1970'));
    nc_attput([newdir filesep files(i,:)],'time','actual_range',[num2str(beg_year), '-01-01T00:00:00 - ', num2str(end_year), '-01-01T00:00:00']);
    nc_varput([newdir filesep files(i,:)],'z',Zf2);
    nc_attput([newdir filesep files(i,:)],'z','definition','Time interpolated altitude based on measured altitude data');
end
