startTime = tic;
EPSG = load('EPSG');
OPT.raw_path     = 'D:\Repositories\OpenEarthRawData\trunk\rijkswaterstaat\vaklodingen\raw\grid';
OPT.netcdf_path  = 'D:\Repositories\OpenEarthRawData\trunk\rijkswaterstaat\vaklodingen\nc';
mkpath(OPT.netcdf_path)
fns = dir(fullfile(OPT.netcdf_path,'*.nc'));
for ii = 1:length(fns)
    delete(fullfile(OPT.netcdf_path,fns(ii).name));
end

fns = dir(fullfile(OPT.raw_path,'*.asc'));

codes = char(fns.name);
[dummy,dummy,position] = unique(codes(:,1:11),'rows');

totalBytes = 0;
bytesRead  = 0;
for ii = 1:length(fns)
    totalBytes = totalBytes + fns(ii).bytes;
end

for ii = 1:max(position);
    subset = find(position==ii);
    fprintf('Reading data ... ');

    %% pre allocate
    time         =  nan(size(subset));
    ncols        =  nan(size(subset));
    nrows        =  nan(size(subset));
    xllcorner    =  nan(size(subset));
    yllcorner    =  nan(size(subset));
    cellsize     =  nan(size(subset));
    nodata_value =  nan(size(subset));
    C            = cell(size(subset));
    %% read contents
    for jj = 1:length(subset);
        time(jj) =                     datenum(fns(subset(jj)).name(12:19),'yyyymmDD'); 
        fid      = fopen(fullfile(OPT.raw_path,fns(subset(jj)).name));
        s        = fgetl(fid);

        if strcmp(s,'/* arc grid file created by Matlab'); 
            s = fgetl(fid); 
        end
                        ncols        (jj) = strread(s,       'ncols %d');
        s = fgetl(fid); nrows        (jj) = strread(s,       'nrows %d');
        s = fgetl(fid); xllcorner    (jj) = strread(s,   'xllcorner %f');
        s = fgetl(fid); yllcorner    (jj) = strread(s,   'yllcorner %f');
        s = fgetl(fid); cellsize     (jj) = strread(s,    'cellsize %f');
        s = fgetl(fid);
        try             nodata_value (jj) = strread(s,'nodata_value %f');
        catch;          nodata_value (jj) = strread(s,'NODATA_value %f'); %#ok<CTCH>
        end
        C{jj} = textscan(fid,'%f');

        fclose(fid);
        bytesRead = bytesRead + fns(subset(jj)).bytes;
        switch fns(subset(jj)).name
            case 'KB115_3938_19950501.asc'
                xllcorner(jj)    =  xllcorner(jj)-10;
                yllcorner(jj)    =  yllcorner(jj)+10; 
            case 'KB115_4140_19950601.asc'
                xllcorner(jj)    =  xllcorner(jj)-10;
                yllcorner(jj)    =  yllcorner(jj)+10; 
            case 'KB116_3938_19950401.asc'
                xllcorner(jj)    =  xllcorner(jj)+10;
                yllcorner(jj)    =  yllcorner(jj)+10;
            case 'KB116_4140_19950401.asc'
                xllcorner(jj)    =  xllcorner(jj)+10;
                yllcorner(jj)    =  yllcorner(jj)+10;
            case 'KB117_3938_19950401.asc'
                xllcorner(jj)    =  xllcorner(jj)+10;
                yllcorner(jj)    =  yllcorner(jj)+10;
            case 'KB117_4140_19950301.asc'
                xllcorner(jj)    =  xllcorner(jj)+10;
                yllcorner(jj)    =  yllcorner(jj)-10;
            case 'KB122_2322_19940401.asc'
                yllcorner(jj)    =  yllcorner(jj)-10;
            case 'KB129_1312_19930301.asc'
                yllcorner(jj)    =  yllcorner(jj)-10;
        end
    end

    if ~all(cellsize==20);
        error('irregular cell size')
    end
    %% reshape data
    x =         min(xllcorner):20:max(xllcorner + cellsize.*(ncols-1));
    y = flipud((min(yllcorner):20:max(yllcorner + cellsize.*(nrows-1)))');

    [xRD,yRD] = meshgrid(x,y);
    [lon,lat] = convertCoordinates(xRD,yRD,EPSG,'CS1.code',28992,'CS2.code',4326);
    Z = nan(length(subset),length(y),length(x));

    for jj = 1:length(subset); 
        x1 = find(x == xllcorner(jj)): 1:find(x == xllcorner(jj))+ncols(jj)-1;
        y1 = find(y == yllcorner(jj)):-1:find(y == yllcorner(jj))-nrows(jj)+1;
        z1 = reshape(C{jj}{1},ncols(jj),nrows(jj));
        z1(z1==nodata_value(jj)) = nan;
        Z(jj,y1,x1) = flipud(z1')/100;
    end

    allNans = all(all(isnan(Z),3),2);
    
    if all(allNans)
        % then skip
        fprintf(2,['creation of ' codes(subset(1),1:10) '.nc skipped due to lack of data. \n']);
    else 
    
        subset       (allNans)     = [];
        Z            (allNans,:,:) = [];
        time         (allNans)     = [];
        ncols        (allNans)     = [];
        nrows        (allNans)     = [];
        xllcorner    (allNans)     = [];
        yllcorner    (allNans)     = [];
        cellsize     (allNans)     = [];
        nodata_value (allNans)     = [];

        %% *** create empty outputfile ***
        % indicate NetCDF outputfile name and create empty structure
        fprintf(['Creating nc-file ',         codes(subset(1),1:10) '.nc ... ']);
        ncfile   = fullfile(OPT.netcdf_path,[codes(subset(1),1:10) '.nc']);
        ncid     = netcdf.create(ncfile,0);
        globalID = netcdf.getConstant('NC_GLOBAL');

        %% *** add global attributes ***
        % netcdf.putAtt(ncid,globalID, 'conventions',     OPT.Conventions);
        % netcdf.putAtt(ncid,globalID, 'title',           OPT.title);
        % netcdf.putAtt(ncid,globalID, 'institution',     OPT.institution);
        % netcdf.putAtt(ncid,globalID, 'source',          OPT.source);
        % netcdf.putAtt(ncid,globalID, 'history',         OPT.history);
        % netcdf.putAtt(ncid,globalID, 'references',      OPT.references);
        % netcdf.putAtt(ncid,globalID, 'comment',         OPT.comment);
        % netcdf.putAtt(ncid,globalID, 'email',           OPT.email);
        % netcdf.putAtt(ncid,globalID, 'version',         OPT.version);
        % netcdf.putAtt(ncid,globalID, 'terms_for_use',   OPT.terms_for_use);
        % netcdf.putAtt(ncid,globalID, 'disclaimer',      OPT.disclaimer);

        %% *** add dimensions ***
        netcdf.defDim(ncid,     'time',            length(subset));
        netcdf.defDim(ncid,     'x',               length(x)     );
        netcdf.defDim(ncid,     'y',               length(y)     );

        %% *** add variables ***

        nc_oe_standard_names('ncid', ncid, 'nc_library', 'matlab', 'varname', {'time'}, 'oe_standard_name', {'time'},                    'dimension', {'time'},      'timezone', '+01:00');
        nc_cf_standard_names('ncid', ncid, 'nc_library', 'matlab', 'varname', {'lon'},  'cf_standard_name', {'longitude'},               'dimension', {'y','x'});
        nc_cf_standard_names('ncid', ncid, 'nc_library', 'matlab', 'varname', {'lat'},  'cf_standard_name', {'latitude'},                'dimension', {'y','x'});
        nc_cf_standard_names('ncid', ncid, 'nc_library', 'matlab', 'varname', {'x'},    'cf_standard_name', {'projection_x_coordinate'}, 'dimension', {'x'});
        nc_cf_standard_names('ncid', ncid, 'nc_library', 'matlab', 'varname', {'y'},    'cf_standard_name', {'projection_y_coordinate'}, 'dimension', {'y'});
        nc_cf_standard_names('ncid', ncid, 'nc_library', 'matlab', 'varname', {'z'},    'cf_standard_name', {'altitude'},                'dimension', {'time','y','x'});

        netcdf.endDef(ncid)

        netcdf.putVar(ncid, netcdf.inqVarID(ncid,'x'   ),                       x);
        netcdf.putVar(ncid, netcdf.inqVarID(ncid,'y'   ),                       y);
        netcdf.putVar(ncid, netcdf.inqVarID(ncid,'time'),time-datenum('1-1-1970'));
        netcdf.putVar(ncid, netcdf.inqVarID(ncid,'lat' ),                     lat);
        netcdf.putVar(ncid, netcdf.inqVarID(ncid,'lon' ),                     lon);
        netcdf.putVar(ncid, netcdf.inqVarID(ncid,'z'   ),                       Z);

        %% *** close nc-file ***
        netcdf.close(ncid);
        fprintf('% 6.2f%% done, ETA is %s\n',...
            (bytesRead/totalBytes)*100,datestr((toc(startTime)*(1/(bytesRead/totalBytes)-1)/3600/24)+now,'HH:MM:SS'))
    end
    %% *** clear some memory ***
    clear Z C lat lon xRD yRD z1
end


