%LDB2NC  rewrites Deltares landboundary file with segment names file (line segments) into NetCDF file
%
%See also: SNCTOOLS, LANDBOUNDARY, POLYSPLIT, POLYJOIN, POLY_FUN

% $Id: ldb2nc.m 
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords: $

%% Initialize

   OPT.fillvalue         = 999.999; % NaNs do work in netcdf API
   OPT.dump              = 1;
   
   OPT.pause             = 0;
   
   OPT.refdatenum        = datenum(0000,0,0); % matlab datenumber convention: A serial date number of 1 corresponds to Jan-1-0000. Gives wring date sin ncbrowse due to different calenders. Must use doubles here.
   OPT.refdatenum        = datenum(1970,1,1); % lunix  datenumber convention
   
   OPT.directory_raw     = pwd;
   OPT.directory_nc      = pwd; 
   OPT.ext               = ''; % to add to output file name before *.nc

%% Set (1)

   OPT.file.landboundary = 'ref20OR2.ldb';
   OPT.file.names        = '';
   OPT.file.nc           = 'refline_OR';
   OPT.title             = 'USGS reference coastline of the state of Oregon';
   

%% 0 Read raw data
   
  [D.xpar,D.ypar]        = landboundary('read',OPT.file.landboundary);% ,999.999
  [D.lon ,D.lat ]        = convertCoordinates(D.xpar,D.ypar,'CS1.code',32610,'CS2.code',4326);
   D.version             = '1.0';
   
%% 0 Create file
   
   outputfile    = fullfile(OPT.directory_nc,[OPT.file.nc,'.nc']);
   
   nc_create_empty (outputfile)
   
%% 1 Add global meta-info to file
% Add overall meta info:
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#description-of-file-contents
   
   nc_attput(outputfile, nc_global, 'title'           , OPT.title);
   nc_attput(outputfile, nc_global, 'institution'     , 'Deltares');
   nc_attput(outputfile, nc_global, 'source'          , 'USGS Woodshole');
   nc_attput(outputfile, nc_global, 'history'         , ['Original filename: ',OPT.file.landboundary,...
                                  ', version:' ,D.version,...
                                  ', tranformation to NetCDF: $HeadURL$']);
   nc_attput(outputfile, nc_global, 'references'      , '<http://openearth.deltares.nl>');
   nc_attput(outputfile, nc_global, 'email'           , 'afarris@usgs.gov');
   
   nc_attput(outputfile, nc_global, 'comment'         , 'Converted to netCDF by Ben de Sonneville');
   nc_attput(outputfile, nc_global, 'version'         , D.version);
   
   nc_attput(outputfile, nc_global, 'Conventions'     , 'CF-1.4');
   nc_attput(outputfile, nc_global, 'CF:featureType'  , '');  % https://cf-pcmdi.llnl.gov/trac/wiki/PointObservationConventions
   
   nc_attput(outputfile, nc_global, 'terms_for_use'   , 'These data can be used freely for research purposes provided that the following source is acknowledged: Deltares.');
   nc_attput(outputfile, nc_global, '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.');
    
   nc_attput(outputfile, nc_global, 'spatial_ref', 'PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],AUTHORITY["EPSG","32610"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]');% UTM zone 10N; see spatialreference.org using wkt. approach used by gdal.

%% 2 Create dimensions
   
   nc_add_dimension(outputfile, 'points'     , length(D.lon))  ; % incl nan-separators
   
%% 3 Create variables
   
   clear nc
   ifld = 0;
   
   % Longitude:
   % http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#longitude-coordinate
   
   ifld = ifld + 1;
   nc(ifld).Name         = 'lon';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'points'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'longitude');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', 'degrees_east');
   nc(ifld).Attribute(3) = struct('Name', 'standard_name'  ,'Value', 'longitude');
   nc(ifld).Attribute(4) = struct('Name', 'comment'        ,'Value', 'segments are separated by NaN');
   
   % Latitude:
   % http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#latitude-coordinate
   
   ifld = ifld + 1;
   nc(ifld).Name         = 'lat';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'points'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'latitude');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', 'degrees_north');
   nc(ifld).Attribute(3) = struct('Name', 'standard_name'  ,'Value', 'latitude');
   nc(ifld).Attribute(4) = struct('Name', 'comment'        ,'Value', 'segments are separated by NaN');
   
%% Coordinate system
%  http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#appendix-grid-mappings
   
%      ifld = ifld + 1;
%    nc(ifld).Name         = 'Oblique_Stereographic';
%    nc(ifld).Nctype       = 'int';
%   %nc(ifld).Dimension    = {'x_cen'}; % no dimension, dummy variable
%    nc(ifld).Attribute(1) = struct('Name', 'grid_mapping_name'            ,'Value', 'Oblique_Stereographic');
%    nc(ifld).Attribute(2) = struct('Name', 'latitude_of_projection_origin','Value', 52.15616055555555);
%    nc(ifld).Attribute(3) = struct('Name', 'false_easting'                ,'Value', 155000.0);
%    nc(ifld).Attribute(4) = struct('Name', 'false_northing'               ,'Value', 463000.0);
%    nc(ifld).Attribute(5) = struct('Name', 'scale_factor'                 ,'Value', 0.9999079);
%    nc(ifld).Attribute(6) = struct('Name', 'central_meridian'             ,'Value', 5.38763888888889);
%    nc(ifld).Attribute(7) = struct('Name', 'comment'                      ,'Value', 'value set to the EPSG code');
%    nc(ifld).Attribute(8) = struct('Name', 'spatial_ref'                  ,'Value', ...
%    ['"PROJCS[\"Amersfoort / RD New\"',...
%         'GEOGCS[\"Amersfoort\",',...
%         '    DATUM[\"Amersfoort\",',...
%         '        SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,',...
%         '            AUTHORITY[\"EPSG\",\"7004\"]],',...
%         '        AUTHORITY[\"EPSG\",\"6289\"]],',...
%         '    PRIMEM[\"Greenwich\",0,',...
%         '        AUTHORITY[\"EPSG\",\"8901\"]],',...
%         '    UNIT[\"degree\",0.01745329251994328,',...
%         '        AUTHORITY[\"EPSG\",\"9122\"]],',...
%         '    AUTHORITY[\"EPSG\",\"4289\"]],',...
%         'UNIT[\"metre\",1,',...
%         '    AUTHORITY[\"EPSG\",\"9001\"]],',...
%         'PROJECTION[\"Oblique_Stereographic\"],',...
%         'PARAMETER[\"latitude_of_origin\",52.15616055555555],',...
%         'PARAMETER[\"central_meridian\",5.38763888888889],',...
%         'PARAMETER[\"scale_factor\",0.9999079],',...
%         'PARAMETER[\"false_easting\",155000],',...
%         'PARAMETER[\"false_northing\",463000],',...
%         'AUTHORITY[\"EPSG\",\"28992\"],',...
%         'AXIS[\"X\",EAST],',...
%         'AXIS[\"Y\",NORTH]]"']);
   
   % x:
   % 
   
   ifld = ifld + 1;
   nc(ifld).Name         = 'x';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'points'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'x coordinate');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', 'm');
   nc(ifld).Attribute(3) = struct('Name', 'standard_name'  ,'Value', 'projection_x_coordinate');
   nc(ifld).Attribute(4) = struct('Name', 'comment'        ,'Value', 'segments are separated by NaN');
   
   % y:
   % 
   
   ifld = ifld + 1;
   nc(ifld).Name         = 'y';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'points'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'y coordinate');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', 'm');
   nc(ifld).Attribute(3) = struct('Name', 'standard_name'  ,'Value', 'projection_y_coordinate');
   nc(ifld).Attribute(4) = struct('Name', 'comment'        ,'Value', 'segments are separated by NaN');

   % Parameters with standard names:
   % * http://cf-pcmdi.llnl.gov/documents/cf-standard-names/standard-name-table/current/standard-name-table/
   
%% 4 Create variables with attibutes
% When variable definitons are created before actually writing the
% data in the next cell, netCDF can nicely fit all data into the
% file without the need to relocate any info.
   
   for ifld=1:length(nc)
       nc_addvar(outputfile, nc(ifld));
   end
   
%% 5 Fill variables
    
   nc_varput(outputfile, 'lon'                  , D.lon);
   nc_varput(outputfile, 'lat'                  , D.lat);
   nc_varput(outputfile, 'x'                    , D.xpar);
   nc_varput(outputfile, 'y'                    , D.ypar);
%    nc_varput(outputfile, 'Oblique_Stereographic', 28992);
   
%% 6 Check
   
   if OPT.dump
       nc_dump(outputfile);
   end
   
%    x        = nc_varget(outputfile,'x');
%    y        = nc_varget(outputfile,'y');
%    fillPolygon([x y],'k','y',1000,1000);
%    lon        = nc_varget(outputfile,'lon');
%    lat        = nc_varget(outputfile,'lat');
%    kmlline(lat ,lon ,'fileName','KMLline_test.kml','lineColor',jet(3) ,'kmlName','horizontal')

   
%% EOF
