%WATERWEG2NETCDF  Script to transform the raw data from Delflandse Kust to NetCDF format
%

% --------------------------------------------------------------------
% Copyright (C) 2004-2009 Delft University of Technology
% Version:      Version 1.0, February 2004
%     Mark van Koningsveld
%
%     m.vankoningsveld@tudelft.nl	
%
%     Hydraulic Engineering Section
%     Faculty of Civil Engineering and Geosciences
%     Stevinweg 1
%     2628CN Delft
%     The Netherlands
%
% This library is free software; you can redistribute it and/or
% modify it under the terms of the GNU Lesser General Public
% License as published by the Free Software Foundation; either
% version 2.1 of the License, or (at your option) any later version.
%
% This library is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% Lesser General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public
% License along with this library; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA
% --------------------------------------------------------------------

% $Id$
% $Date$
% $Author$
% $Revision$

clear all; close all; clc

cd(fileparts(mfilename('fullpath')))

%% *** get data ***
filespath = [fileparts(fileparts(mfilename('fullpath'))) filesep 'raw' filesep];
fns       = dir([filespath '*.xyz']);
files     = {fns.name}';

if dbstate
    files = {
        'GM_Update_NL IN RD WATERWEG_20090203.xyz';
        'GM_Update_NL IN RD WATERWEG_20090204.xyz'};
end

%% determine the outline of the dataset
[minx, maxx, miny, maxy] = deal([]);
for i = 1:length(files)
    disp(num2str(i))
    data_temp = load([filespath files{i}]);
    minx = min([min(data_temp(:,1)) minx]);
    maxx = max([max(data_temp(:,1)) maxx]);
    miny = min([min(data_temp(:,2)) miny]);
    maxy = max([max(data_temp(:,2)) maxy]);
end

%% construct x and y vectors
data.x = minx:2.5:maxx;
data.y = [miny:2.5:maxy]'; %#ok<NBRAK>

%% *** create empty outputfile

% indicate NetCDF outputfile name and create empty structure
outputfile = 'Waterweg.nc';
nc_create_empty(outputfile);

% add attributes global to the dataset
nc_attput(outputfile, nc_global, 'convention',  'CF-1.4');
nc_attput(outputfile, nc_global, 'title',       'Delflandse Kust');
nc_attput(outputfile, nc_global, 'institution', 'Van Oord & Boskalis');
nc_attput(outputfile, nc_global, 'comment',     'Data surveyed by Van Oord & Boskalis for the project Delflandse Kust');
nc_attput(outputfile, nc_global, 'source',      'Bathymetry measured with survey vessel.');
nc_attput(outputfile, nc_global, 'history',     'Created with: $Id$ $HeadURL$');
nc_attput(outputfile, nc_global, 'references',  'No reference material available');

%% *** add dimensions ***
nc_add_dimension(outputfile,     'x',           size(data.x,2))
nc_add_dimension(outputfile,     'y',           size(data.y,1))
nc_add_dimension(outputfile,     't',           0);             % 0 = 'unlimited'

%% *** add variables ***
% add variable: coordinate system reference (this variable contains all projection information, e.g. for use in ArcGIS)
crsVariable = struct(...
    'Name', 'crs', ...
    'Nctype', 'int', ...
    'Dimension', {{}}, ...
    'Attribute', struct( ...
        'Name', ...
        	{'spatial_ref'}, ...
        'Value', ...
        	{'COMPD_CS["Amersfoort / RD New + NAP",PROJCS["Amersfoort / RD New",GEOGCS["Amersfoort",DATUM["Amersfoort",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[565.04,49.91,465.84,-0.40939438743923684,-0.35970519561431136,1.868491000350572,0.8409828680306614],AUTHORITY["EPSG","6289"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.017453292519943295],AXIS["Geodetic latitude",NORTH],AXIS["Geodetic longitude",EAST],AUTHORITY["EPSG","4289"]],PROJECTION["Oblique Stereographic",AUTHORITY["EPSG","9809"]],PARAMETER["central_meridian",5.387638888888891],PARAMETER["latitude_of_origin",52.15616055555556],PARAMETER["scale_factor",0.9999079],PARAMETER["false_easting",155000.0],PARAMETER["false_northing",463000.0],UNIT["m",1.0],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","28992"]],VERT_CS["Normaal Amsterdams Peil",VERT_DATUM["Normaal Amsterdams Peil",2005,AUTHORITY["EPSG","5109"]],UNIT["m",1.0],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","5709"]],AUTHORITY["EPSG","7415"]]'} ...
            ) ...
    );
nc_addvar(outputfile, crsVariable);

% add variable: time
tVariable = struct(...
    'Name', 't', ...
    'Nctype', 'float', ...
    'Dimension', {{'t'}}, ...
    'Attribute', struct( ...
        'Name', ...
            {'standard_name', 'long_name',  'units'}, ...
        'Value', ...
            {'time', 'time', 'days since 0000-01-01 00:00:00 +1:00'} ...
         ) ...
    );
nc_addvar(outputfile, tVariable);

% add variable: x-coordinate information
xVariable = struct(...
    'Name', 'x', ...
    'Nctype', 'float', ...
    'Dimension', {{'x'}}, ...
    'Attribute', struct( ...
        'Name', ... 
            {'long_name', 'units'}, ...
        'Value', ... 
            {'x-coordinate', 'meter'} ...
        ) ...
    );
nc_addvar(outputfile, xVariable);

% add variable: y-coordinate information
yVariable = struct(...
    'Name', 'y', ...
    'Nctype', 'float', ...
    'Dimension', {{'y'}}, ...
    'Attribute', struct( ...
        'Name', ... 
            {'long_name', 'units'}, ...
        'Value', ... 
            {'y-coordinate', 'meter'} ...
        ) ...
    );
nc_addvar(outputfile, yVariable);

% add variable: z-coordinate information
zVariable = struct(...
    'Name', 'z', ...
    'Nctype', 'float', ...
    'Dimension', {{'t','y','x'}}, ...
    'Attribute', struct( ...
        'Name', ...
            {'standard_name', 'long_name', 'units', 'positive', 'grid_mapping'}, ...
        'Value', ...
            {'surface_altitude', 'elevation', 'meter', 'up', 'crs'} ...
        ) ...
    );
nc_addvar(outputfile, zVariable);

%% *** put data into nc file ***
nc_varput(outputfile, 'x', data.x)
nc_varput(outputfile, 'y', data.y)

for j = 1:length(files)
    disp(['Processing file: ' num2str(j)])
    data_temp = load([filespath files{j}]);        
    data.z = zeros(size(data.y,1), size(data.x,2))*nan;
    for i = 1:size(data_temp,1)
        ii = find(data_temp(i,1) == data.x);
        jj = find(data_temp(i,2) == data.y);
        data.z(jj,ii) = data_temp(i,3);
    end
    [start_idx, end_idx, extents, matches, tokens, names, splits]  = regexp(files{j}, '\d+'); % 20090311.xyz
    data.t = datenum(matches, 'yyyymmdd');
    
    nc_varput(outputfile, 't', data.t, j-1, 1)
    nc_varput(outputfile, 'z', data.z, [j-1 0 0], [1 size(data.z)])
end

%% dump outputfile
nc_dump(outputfile)

