%function gshhs2nc
%GSHHS2NC  rewrites GSHHS binary files into NetCDF file
%
%See also: M_MAP, SNCTOOLS, LANDBOUNDARY, POLYSPLIT, POLYJOIN, POLY_FUN

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

%% Initialize

   OPT.fillvalue         = nan; % NaNs do work in netcdf API
   OPT.dump              = 0;
   
   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     = '../raw/gshhs_1.10/';
   OPT.directory_nc      = '../processed/';
   OPT.ext               = ''; % to add to output file name before *.nc
   
  %OPT.fname             = 'gshhs_f.b'; OPT.resolution = 'fine (finest out of 5 available resolutions)'; % does not work with java: memory issues
  %OPT.fname             = 'gshhs_h.b'; OPT.resolution = 'high (2nd finest out of 5 available resolutions)';
  %OPT.fname             = 'gshhs_i.b'; OPT.resolution = 'intermediate  (middle out of 5 available resolutions)';
  %OPT.fname             = 'gshhs_l.b'; OPT.resolution = 'low (2nd crudest out of 5 available resolutions)';
   OPT.fname             = 'gshhs_c.b'; OPT.resolution = 'crude (crudest out of 5 available resolutions)';

%% Load

   OPT.title             = 'GSHHS: A Global Self-consistent, Hierarchical, High-resolution Shoreline Database';

   D = gshhs_read(OPT.fname(7),[OPT.directory_raw,OPT.fname]);

%% 0 Create file
   
   outputfile    = fullfile(OPT.directory_nc,[filename(OPT.fname),'.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,' ',OPT.resolution]);%,' resolution'
   nc_attput(outputfile, nc_global, 'institution'       , 'Paul Wessel, SOEST, University of Hawaii, Honolulu, HI & Walter H. F. Smith, NOAA Geosciences Lab, National Ocean Service, Silver Spring, MD');
   nc_attput(outputfile, nc_global, 'source'            , 'Paul Wessel, SOEST, University of Hawaii, Honolulu, HI & Walter H. F. Smith, NOAA Geosciences Lab, National Ocean Service, Silver Spring, MD');
   nc_attput(outputfile, nc_global, 'history'           , ['Original filename: ',OPT.fname,...
                                                           ', tranformation to NetCDF: $HeadURL$']);
   nc_attput(outputfile, nc_global, 'references'        , '<http://www.soest.hawaii.edu/wessel/gshhs/gshhs.html>; Wessel, P., and W. H. F. Smith, A Global Self-consistent, Hierarchical, High-resolution Shoreline Database, J. Geophys. Res., 101, 8741-8743, 1996.');
   nc_attput(outputfile, nc_global, 'email'             , '');
   
   nc_attput(outputfile, nc_global, 'comment'           , 'Not original GSHHS data: function mu_coast.m from the m_map toolbox (http://www.eos.ubc.ca/~rich/map.html) of Rich Pawlowicz was used to solve issues with data at 180/-180 transition and Antarctica. This resules in multiple occurences of some objects, hence npoints > n_unique.');
   nc_attput(outputfile, nc_global, 'version'           , ['GSHHS_DATA_VERSION:',num2str(D.GSHHS_DATA_VERSION),',GSHHS_PROG_VERSION:',D.GSHHS_PROG_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: GSHHS.');
   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, 'GSHHS_DATA_VERSION', D.GSHHS_DATA_VERSION);
   nc_attput(outputfile, nc_global, 'GSHHS_PROG_VERSION', D.GSHHS_PROG_VERSION);
   nc_attput(outputfile, nc_global, 'n_unique'          , D.n_unique          );

%% Add discovery information (test):

   %  http://www.unidata.ucar.edu/projects/THREDDS/tech/catalog/InvCatalogSpec.html
  
   nc_attput(outputfile, nc_global, 'geospatial_lat_min'         , min(D.lat));
   nc_attput(outputfile, nc_global, 'geospatial_lat_max'         , max(D.lat));
   nc_attput(outputfile, nc_global, 'geospatial_lon_min'         , min(D.lon));
   nc_attput(outputfile, nc_global, 'geospatial_lon_max'         , max(D.lon));
  %nc_attput(outputfile, nc_global, 'time_coverage_start'        , '');
  %nc_attput(outputfile, nc_global, 'time_coverage_end'          , '');
   nc_attput(outputfile, nc_global, 'geospatial_lat_units'       , 'degrees_north');
   nc_attput(outputfile, nc_global, 'geospatial_lon_units'       , 'degrees_east' );

%% 2 Create dimensions
   
   nc_add_dimension(outputfile, 'npoints'           , length(D.lon)) ; % incl nan-separators
   nc_add_dimension(outputfile, 'closed_segments'   , length(D.area)); % incl nan-separators
   nc_add_dimension(outputfile, 'segment_separators', length(D.k))   ; % 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    = {'npoints'};
   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', 'all segments are closed contours, 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    = {'npoints'};
   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', 'all segments are closed contours, separated by NaN');
   
   % Parameters with standard names:
   % * http://cf-pcmdi.llnl.gov/documents/cf-standard-names/standard-name-table/current/standard-name-table/
   
   % Other parameters

   ifld = ifld + 1;
   nc(ifld).Name         = 'area';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'closed_segments'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'segment patch areas (>0 for land, <0 for lakes and inland seas)');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', 'degrees_north*degrees_east');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', 'Area=zeros(length(k)-1,1);for i=1:length(k)-1,x=lon([k(i)+1:(k(i+1)-1) k(i)+1]);y=lat([k(i)+1:(k(i+1)-1) k(i)+1]);nl=length(x);Area(i)=sum(diff(x).*(y(1:nl-1)+y(2:nl))/2 );end;');

   ifld = ifld + 1;
   nc(ifld).Name         = 'k';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'segment_separators'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'segment separator indices');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', '#');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', 'k=find(isnan(lon(:,1)))');

   ifld = ifld + 1;
   nc(ifld).Name         = 'id';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'closed_segments'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'segment id');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', '');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', '');
   nc(ifld).Attribute(4) = struct('Name', 'coordinates'    ,'Value', 'lat lon');

   ifld = ifld + 1;
   nc(ifld).Name         = 'n';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'closed_segments'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'number of points per segment');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', '');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', '');
   nc(ifld).Attribute(4) = struct('Name', 'coordinates'    ,'Value', 'lat lon');

   ifld = ifld + 1;
   nc(ifld).Name         = 'level';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'closed_segments'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'level of segment');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', '');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', '1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake');
   nc(ifld).Attribute(4) = struct('Name', 'coordinates'    ,'Value', 'lat lon');

   ifld = ifld + 1;
   nc(ifld).Name         = 'version';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'closed_segments'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'version of segment');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', '');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', '4 for GSHHS version 1.4');
   nc(ifld).Attribute(4) = struct('Name', 'coordinates'    ,'Value', 'lat lon');

   ifld = ifld + 1;
   nc(ifld).Name         = 'greenwich';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'closed_segments'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'whether segment crosses greenwich');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', '');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', 'Greenwich is 1 if Greenwich is crossed');
   nc(ifld).Attribute(4) = struct('Name', 'coordinates'    ,'Value', 'lat lon');

   ifld = ifld + 1;
   nc(ifld).Name         = 'source';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'closed_segments'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'segment source');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', '');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', '0 = CIA WDBII, 1 = WVS');
   nc(ifld).Attribute(4) = struct('Name', 'coordinates'    ,'Value', 'lat lon');

   ifld = ifld + 1;
   nc(ifld).Name         = 'west';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'closed_segments'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'segment source');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', '');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', '0 = CIA WDBII, 1 = WVS');
   nc(ifld).Attribute(4) = struct('Name', 'coordinates'    ,'Value', 'lat lon');

   ifld = ifld + 1;
   nc(ifld).Name         = 'east';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'closed_segments'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'segment source');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', '');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', '0 = CIA WDBII, 1 = WVS');
   nc(ifld).Attribute(4) = struct('Name', 'coordinates'    ,'Value', 'lat lon');

   ifld = ifld + 1;
   nc(ifld).Name         = 'south';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'closed_segments'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'segment source');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', '');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', '0 = CIA WDBII, 1 = WVS');
   nc(ifld).Attribute(4) = struct('Name', 'coordinates'    ,'Value', 'lat lon');

   ifld = ifld + 1;
   nc(ifld).Name         = 'north';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {'closed_segments'};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'segment source');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', '');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', '0 = CIA WDBII, 1 = WVS');
   nc(ifld).Attribute(4) = struct('Name', 'coordinates'    ,'Value', 'lat lon');

   ifld = ifld + 1;
   nc(ifld).Name         = 'n_unique';
   nc(ifld).Nctype       = 'float'; % no double needed
   nc(ifld).Dimension    = {};
   nc(ifld).Attribute(1) = struct('Name', 'long_name'      ,'Value', 'number of unique segments');
   nc(ifld).Attribute(2) = struct('Name', 'units'          ,'Value', '');
   nc(ifld).Attribute(3) = struct('Name', 'comment'        ,'Value', 'Dimension npoints > n_unique because function mu_coast.m from the m_map toolbox (http://www.eos.ubc.ca/~rich/map.html) of Rich Pawlowicz was used to solve issues with data at 180/-180 transition and Antarctica. This resules in multiple occurences of some objects.');
   
%% 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 (snctools now requires 1D variables to have correct shape, so we need to do (:))
   
   nc_varput(outputfile, 'lon' , D.lon(:));
   nc_varput(outputfile, 'lat' , D.lat(:));
   nc_varput(outputfile, 'k'   , D.k(:));
   nc_varput(outputfile, 'area', D.area(:));
   
   nc_varput(outputfile,       'id',        D.id(:));
   nc_varput(outputfile,        'n',         D.n(:));
   nc_varput(outputfile,    'level',     D.level(:));
   nc_varput(outputfile,  'version',   D.version(:));
   nc_varput(outputfile,'greenwich', D.greenwich(:));
   nc_varput(outputfile,   'source',    D.source(:));
   nc_varput(outputfile,     'west',      D.west(:));
   nc_varput(outputfile,     'east',      D.east(:));
   nc_varput(outputfile,    'south',     D.south(:));
   nc_varput(outputfile,    'north',     D.north(:));
   
%% 6 Check
   
   if OPT.dump
       nc_dump(outputfile);
   end

%% Check by plotting

   %% prevent memory issues

   USE_JAVA = getpref('SNCTOOLS','USE_JAVA');
   setpref('SNCTOOLS','USE_JAVA',0);

   N.lon = nc_varget(outputfile,'lon');
   N.lat = nc_varget(outputfile,'lat');

   plot(N.lon,N.lat)
   pausedisp
   close all
   clear N

   %% prevent memory issues

   setpref('SNCTOOLS','USE_JAVA',USE_JAVA);
   
%% EOF
