function grainsize2netcdf(varargin)
%GRAINSIZE2NETCDF  create netcdf file of VTV grainsize data
%
%   Read relevant grain size data from Excel file and store it in a netCDF
%   file. This are the grain sizes as used for the VTV (guidelines for the
%   safety assessment of the Dutch dune coast)
%
%   Syntax:
%   grainsize2netcdf(varargin)
%
%   Input:
%   varargin  =
%
%   Example
%   grainsize2netcdf
%
%   See also

%   --------------------------------------------------------------------
%   Copyright (C) 2009 Delft University of Technology
%       C.(Kees) den Heijer
%
%       C.denHeijer@TUDelft.nl	
%
%       Faculty of Civil Engineering and Geosciences
%       P.O. Box 5048
%       2600 GA 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
%   or http://www.gnu.org/licenses/licenses.html, http://www.gnu.org/, http://www.fsf.org/
%   --------------------------------------------------------------------

% Created: 16 Feb 2009
% Created with Matlab version: 7.4.0.287 (R2007a)

% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords:

%% initialise default arguments
% define default arguments and apply input
OPT = struct(...
...    'filename', fullfile(fileparts(which(mfilename)), 'grainsize2.nc'),...
    'sourcefile', 'GrainSizes.xls',...
    'EPSGcode', 28992);

OPT = setproperty(OPT, varargin{:});

%% meta data
data_title = 'Grain size';
data_summary = 'Grain size measurements along the Dutch coast, to be used for the dune safety assessment (VTV).';
data_keywords = 'grain size, coast, VTV';
path_main = '.';
path_sub_netcdf = '';
path_sub_kml = 'rijkswaterstaat';

publisher_email = 'Kees.denHeijer@deltares.nl';
institution_longname = 'Rijkswaterstaat';
institution_shortname = 'Rijkswaterstaat';

creator_name = institution_longname;
creator_url = 'www.rijkswaterstaat.nl';
creator_email = 'info@rijkswaterstaat.nl';

processing_level = 'final';

publisher_name = getenv('USERNAME');
publisher_domain = getenv('USERDOMAIN');
publisher_computer = getenv('COMPUTERNAME');
publish_script_id = '$Id$';
path_data_title = lower(strrep(data_title, ' ', '_'));
path_netcdf = fullfile(path_main, 'nc', path_sub_netcdf);
path_kml = fullfile(path_main, 'kml', path_data_title, path_sub_kml);

%% read source data from Excel file
Grainsizes = xlsread(OPT.sourcefile);
Grainsizes(:,1) = [];
% filter for empty lines
id = ~isnan(Grainsizes(:,1)) & ~isnan(Grainsizes(:,3));
Grainsizes = Grainsizes(id,:);
% express grain sizes in meters
Grainsizes(:,end-2:end) = Grainsizes(:,end-2:end) * 1e-6;

% find area names, corresponding with the area codes, and derive the
% maximum stringsize
areacodes = Grainsizes(:,1);
areaname = cellfun(@rws_kustvak, num2cell(areacodes));

stringsize = max(cellfun(@length, areaname));

id = Grainsizes(:,1)*1e6 + Grainsizes(:,2)*100;

%% read jarkus transect info, to get coordinates of RSP points
tr = jarkus_transects(...
    'output', {'x' 'y' 'lat' 'lon' 'id'},...
    'cross_shore', 0);

[d ix] = cellfun(@(x) min(abs(x - tr.id)), num2cell(id));
if any(d~=0)
    fprintf('%i locations correspond to jarkus transects\n', sum(d==0));
    fprintf('The remaining %i locations deviate from the closest transects:\n', sum(d~=0));
    fprintf('\t%.1f m\n', d(d~=0));
end

loc = cellfun(@(x) x(ix), struct2cell(rmfield(tr, 'url')), 'uniformoutput', false);
[x y lat lon id] = deal(loc{:});

%% create schema
dimstruct        = nccreateDimstruct(...
    'Name', 'alongshore',...
    'Length', length(id));
dimstruct(end+1) = nccreateDimstruct(...
    'Name','stringsize',...
    'Length', stringsize);

varstruct = nccreateVarstruct(...
    'Name','id',...
    'Dimensions',{'alongshore'},...
    'Datatype', 'int32',...
    'Attributes',{'long_name', 'identifier',...
    'comment', 'sum of area code (*1e6) and alongshore coordinate'});
varstruct(end+1) = nccreateVarstruct(...
    'Name','areacode',...
    'Dimensions',{'alongshore'},...
    'Datatype', 'int16',...
    'Attributes',{'long_name', 'area code',...
    'comment', 'codes for the 17 areas as defined by Rijkswaterstaat'});
varstruct(end+1) = nccreateVarstruct(...
    'Name','areaname',...
    'Dimensions',{'alongshore' 'stringsize'},...
    'Datatype', 'char',...
    'Attributes',{'long_name', 'area name',...
    'comment', 'names for the 17 areas as defined by Rijkswaterstaat'});

varstruct(end+1) = nccreateVarstruct(...
    'Name','alongshore',...
    'Dimensions',{'alongshore'},...
    'Attributes',{'long_name', 'alongshore coordinate',...
    'units', 'm'});
varstruct(end+1) = nccreateVarstruct(...
    'Name','meanD50',...
    'Dimensions',{'alongshore'},...
    'Attributes',{'long_name', 'mean of D50 grain size',...
    'units', 'm',...
    'actual_range',minmax(Grainsizes(:,3)')});
varstruct(end+1) = nccreateVarstruct(...
    'Name','sigmaD50',...
    'Dimensions',{'alongshore'},...
    'Attributes',{'long_name', 'standard deviation of D50 grain size',...
    'units', 'm',...
    'actual_range',minmax(Grainsizes(:,4)')});
varstruct(end+1) = nccreateVarstruct(...
    'Name','calcD50',...
    'Dimensions',{'alongshore'},...
    'Attributes',{'long_name', 'calculation value of D50 grain size',...
    'units', 'm',...
    'actual_range',minmax(Grainsizes(:,5)')});

varstruct(end+1) = nccreateVarstruct_standardnames_cf('projection_x_coordinate',...
    'Name', 'x',...
    'Dimensions', {'alongshore'},...
    'Attributes',{'grid_mapping','crs',...
    'actual_range',minmax(x')});
varstruct(end+1) = nccreateVarstruct_standardnames_cf('projection_y_coordinate',...
    'Name', 'y',...
    'Dimensions', {'alongshore'},...
    'Attributes',{'grid_mapping','crs',...
    'actual_range',minmax(y')});
varstruct(end+1) = nccreateVarstruct_crs(OPT.EPSGcode);
varstruct(end+1) = nccreateVarstruct_standardnames_cf('latitude',...
    'Name','lat',...
    'Dimensions',{'alongshore'},...
    'DeflateLevel',1,...
    'Attributes', {'actual_range', minmax(lat')});
varstruct(end+1) = nccreateVarstruct_standardnames_cf('longitude',...
    'Name','lon',...
    'Dimensions',{'alongshore'},...
    'DeflateLevel',1,...
    'Attributes', {'actual_range', minmax(lon')});
%     coordinates = 'lat lon';

ncattributes = {...
    ... % Identification
...    'id', 'none',...
    'naming_authority', 'deltares.nl',... % based on reverse DNS lookup (http://remote.12dt.com/)
    'Metadata_Conventions', 'CF-1.5',...
...    'metadata_link', 'none',...
    ... Text search
    'title', data_title,...
    'summary', data_summary,...
    'keywords', data_keywords,...
	'keywords_vocabulary', 'http://www.eionet.europa.eu/gemet',...
	'standard_name_vocabulary', 'http://cf-pcmdi.llnl.gov/documents/cf-standard-names/',...
    'history', sprintf('Created on %s by %s on computer %s\\%s with script %s',...
    datestr(now,31), publisher_name, publisher_domain, publisher_computer, publish_script_id),...
    'comment', sprintf(['TAW, May 1984. Guide on the assessment of dune safety as water defence (Leidraad voor de beoordeling van de veiligheid van duinen als waterkering), Technical Advisory Committee on Water Defences. http://repository.tudelft.nl/view/hydro/uuid%3A0348afd3-d846-4be9-934f-abb5f2d4e0b2/\n'...
    'Voorschrift toetsen op Veiligheid Primaire Waterkeringen (VTV2006) http://repository.tudelft.nl/view/hydro/uuid%3A15d29d22-862b-418e-903d-d90c01f04983/']),...
    'institution', institution_longname,...
	'cdm_data_type', 'points',...
...	'acknowledgment', 'none',...
    'creator_name', creator_name,...
    'creator_url', creator_url,...
    'creator_email', creator_email,...
    'publisher_name', publisher_name,...
    'publisher_url', 'www.deltares.nl',...
    'publisher_email', publisher_email,...
    ... % Other attributes
	'processing_level', processing_level,...
	'license',[sprintf('These data can be used freely for research purposes provided that the following source is acknowledged: %s. ', institution_shortname)...
                'disclaimer: This data is made available 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.'],...
    'projectionCoverage_x', minmax(x'),...
    'projectionCoverage_y', minmax(y'),...
    'geospatialCoverage_northsouth', minmax(lat'),...
    'geospatialCoverage_eastwest', minmax(lon'),...
    'geospatial_lon_units', 'degrees_east',...
    'geospatial_lon_min', min(lon),...
    'geospatial_lon_max', max(lon),...
    'geospatial_lat_units', 'degrees_north',...
    'geospatial_lat_min', min(lat),...
    'geospatial_lat_max', max(lat)};

schema = nccreateSchema(dimstruct,varstruct,...
    'Filename','a',...
    'Attributes', ncattributes,...
    'Format', 'classic');

%% write data
ncfile = fullfile(path_netcdf, 'grainsize_vtv.nc');

mkpath(path_netcdf)

if exist(ncfile, 'file')
    delete(ncfile)
end
ncwriteschema(ncfile, schema);

ncwrite(ncfile,'id', id(:));
ncwrite(ncfile,'areacode', Grainsizes(:,1));
ncwrite(ncfile,'areaname', char(areaname));
ncwrite(ncfile,'alongshore', Grainsizes(:,2)*100);
ncwrite(ncfile,'meanD50', Grainsizes(:,3));
ncwrite(ncfile,'sigmaD50', Grainsizes(:,4));
ncwrite(ncfile,'calcD50', Grainsizes(:,5));
ncwrite(ncfile,'x', x);
ncwrite(ncfile,'y', y);
ncwrite(ncfile,'lat',lat);
ncwrite(ncfile,'lon',lon);

%% check
meanD50 = ncread(ncfile, 'meanD50');
sigmaD50 = ncread(ncfile, 'sigmaD50');
calcD50 = ncread(ncfile, 'calcD50');

ok1 = all(Grainsizes(:,3) == meanD50)
ok2 = all(Grainsizes(:,4) == sigmaD50)
ok3 = all(Grainsizes(:,5) == calcD50)

%% kml
locname = cellfun(@num2str, num2cell(id), 'uniformoutput', false);
% html = cellfun(@(x,y,z) sprintf('mean D50: %.0f &mu;m</br>sigma D50: %.0f &mu;m</br>calculation D50: %.0f &mu;m', x, y, z),...
%     num2cell(meanD50*1e6), num2cell(sigmaD50*1e6), num2cell(calcD50*1e6), 'uniformoutput', false);

html = cellfun(@(varargin) sprintf(['<table border="0">'...
    '<tr>'...
    '<td>Transect id</td>'...
    '<td>%.0f</td>'...
    '</tr>'...
    '<tr>'...
    '<td>Area (kustvak)</td>'...
    '</tr>'...
    '<tr>'...
    '<td>Name</td>'...
    '<td>%s</td>'...
    '</tr>'...
    '<tr>'...
    '<td>Code</td>'...
    '<td>%i</td>'...
    '</tr>'...
    '<tr>'...
    '<td>Latitude</td>'...
    '<td>%.2f ?N</td>'...
    '</tr>'...
    '<tr>'...
    '<td>Longitude</td>'...
    '<td>%.2f ?E</td>'...
    '</tr>'...
    '<tr>'...
    '<td>Mean D<sub>50</sub></td>'...
    '<td>%.0f &mu;m</td>'...
    '</tr>'...
    '<tr>'...
    '<td>Sigma D<sub>50</sub></td>'...
    '<td>%.0f &mu;m</td>'...
    '</tr>'...
    '<tr>'...
    '<td>Calculation D<sub>50</sub></td>'...
    '<td>%.0f &mu;m</td>'...
    '</tr>'...
    '</table>'...
    ...
    '<h3>Meta information</h3>'...
    '<table border="0" padding="0" width="200">'...
    '<tr><td>THREDDS server</td>'...
    '<td><a href="http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/grainsize_vtv/grainsize_vtv.nc.html">meta-info     </a></td></tr>'...
    '<tr><td>THREDDS server</td>'...
    '<td><a href="http://opendap.deltares.nl/thredds/fileServer/opendap/rijkswaterstaat/grainsize_vtv/grainsize_vtv.nc">netCDF via ftp</a></td></tr>'...
    '</table>'...
    ], varargin{:}),...
    num2cell(id'),...
    areaname, num2cell(Grainsizes(:,1)),...
    num2cell(lat), num2cell(lon),...
    num2cell(meanD50*1e6), num2cell(sigmaD50*1e6), num2cell(calcD50*1e6), 'uniformoutput', false);

kmlfile = fullfile(path_kml, 'grainsize_vtv.kml');
mkpath(path_kml)
KMLscatter(lat, lon, meanD50*1e6,...
    'fileName', kmlfile,...
    'name',locname,...
...    'snippets', snippets,...
    'html', html,...
    'scalenormalState',.5,...
    'scalehighlightState',1.5)

strrep_in_files(kmlfile,...
    '<name>grainsize_vtv</name>', '<name>VTV grain size</name>')
strrep_in_files(kmlfile,...
    '<description></description>', ['<description><![CDATA[<a href="http://www.rijkswaterstaat.nl">Rijkswaterstaat</a> grain sizes provided by <a href="http://www.openearth.nl"> OpenEarthTools</a><br>'...
    'see <a href="http://repository.tudelft.nl/view/hydro/uuid%3A15d29d22-862b-418e-903d-d90c01f04983/">VTV2006</a> and '...
    '<a href="http://repository.tudelft.nl/view/hydro/uuid%3A0348afd3-d846-4be9-934f-abb5f2d4e0b2/">TAW1984</a> for more info'...
    ']]></description>'])

KML2kmz(kmlfile, fullfile(path_kml, 'grainsize_vtv_ver_lft.png'))
delete(kmlfile, fullfile(path_kml, 'grainsize_vtv_ver_lft.png'))
