%GEBCO2CF  rewrite 1 arc-minute/30 arc-second GEBCO worldwide DEM vector to netCDF-CF grid
%
% The bare netCDF files should be downloaded manually from www.gebco.net
%
% Example on acces of this data: OpenEarthTools: OPeNDAP_subsetting_with_Matlab_tutorial.m
%
%See also: SNCTOOLS, NETCDF, NC_CF_GRID, OPENDAP_SUBSETTING_WITH_MATLAB_TUTORIAL
%          http://geoport.whoi.edu/thredds/bathy_catalog.html

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

OPT.netcdfversion = 3;  % 4 does not make it any smaller for 'topo' variable, but for 'SID' variable it would
OPT.gebco         = 30; % 1 = one arc minute, 30 = 30 arc second

if OPT.gebco==1
   ncfileraw      = '..\raw\gridone.nc';
   ncfile         = '..\raw\gebco_1min.nc';
   OPT.references = 'The GEBCO One Minute Grid, version 2.0, http://www.gebco.net';
   OPT.version    = '2.0';
   OPT.chunksize  = 1000; % for 1GB memory on win32 machine
elseif OPT.gebco==30
   ncfileraw      = '..\raw\gebco_08.nc';
   ncfileSID      = '..\raw\gebco_SID.nc';
   ncfile         = '..\raw\gebco_30sec.nc';
   OPT.references = 'The GEBCO_08 Grid version 20091120, http://www.gebco.net';
   OPT.version    = '20091120';
   OPT.chunksize  =  500; % for 1GB memory on win32 machine
end

%% The GEBCO netCDF file will be structured after the 6 
%  global DEMs available through web services via the 
%  the USGS/WHOI OPeNDAP server at http://geoport.whoi.edu/thredds/bathy_catalog.html, e.g.:
% smith_sandwell_v11 {
% 
% dimensions:
%	lon = 43200 ;
%	lat = 21600 ;
%//	lon = 21600 ; // gridone.nc
%//	lat = 17280 ; // gridone.nc
% 
% variables:
% 	// Preference 'PRESERVE_FVD':  false,
% 	// dimensions consistent with ncBrowse, not with native MATLAB netcdf package.
% 	single lon(lon), shape = [21600]
% 		lon:long_name = "Uniformly spaced longitudes (-179.9917E - 179.9917E.01667E)" 
% 		lon:units = "degrees_east" 
% 	single lat(lat), shape = [17280]
% 		lat:long_name = "Mercator spaced latitudes (80.738N - 80.738S)" 
% 		lat:units = "degrees_north" 
% 	int16 topo(lat,lon), shape = [17280 21600]
% 		topo:_CoordinateAxes = "lat lon " 
% 		topo:units = "meters" 
% 		topo:long_name = "topography" 
% 
% //global Attributes:
% 		:Conventions = "COARDS" 
% 		:title = "Smith and Sandwell v11.1 Global DEM (1 min)" 
% 		:institution = "University of California San Diego" 
% 		:history = " 28-Apr-2008: Converted to NetCDF using gdal_translate from FWTOOLS, then rearranged from 0:360 to -180:180 using Matlab (Rich Signell: rsignell@usgs.gov) " 
% 
% }

%% load: the GEBCO grid as bare netCDF
%   NetCDF-3 Classic gebco_08.nc { 
%// NetCDF-3 Classic gridone.nc  {

% 
% dimensions:
%	side = 2 ;
%	xysize = 933120000 ; 
%//	xysize = 233312401 ; // gridone.nc
% 
% variables:
%	// Preference 'PRESERVE_FVD':  false,
%	// dimensions consistent with ncBrowse, not with native MATLAB netcdf package.
%	double x_range(side), shape = [2]
%		x_range:units = "user_x_unit" 
%	double y_range(side), shape = [2]
%		y_range:units = "user_y_unit" 
%	int16 z_range(side), shape = [2]
%		z_range:units = "user_z_unit" 
%	double spacing(side), shape = [2]
%	int32 dimension(side), shape = [2]
%//	int16 dimension(side), shape = [2] // gridone.nc
%	int16 z(xysize), shape = [933120000]
%//	int16 z(xysize), shape = [233312401] // gridone.nc
%		z:scale_factor = 1 
%		z:add_offset = 0 
%		z:node_offset = 1 d // 1 = grid line registration
%//		z:node_offset = 0 d // 0 = pixel registration // gridone.nc
%
% //global Attributes:
%		:title = "GEBCO_08 Grid" 
%//		:title = "GEBCO One Minute Grid" 
%		:source = "20100927" 
%//		:source = "2.00" 
%
% }

   NCidraw     = netcdf.open(ncfileraw,netcdf.getConstant('NOWRITE'));

   varid = netcdf.inqVarID(NCidraw,'x_range'  );G.lon_range   = netcdf.getVar(NCidraw,varid);
   varid = netcdf.inqVarID(NCidraw,'y_range'  );G.lat_range   = netcdf.getVar(NCidraw,varid);
   varid = netcdf.inqVarID(NCidraw,'spacing'  );G.spacing     = netcdf.getVar(NCidraw,varid);
   varid = netcdf.inqVarID(NCidraw,'dimension');G.dimension   = netcdf.getVar(NCidraw,varid,'double');
   varid = netcdf.inqVarID(NCidraw,'z'        );G.node_offset = netcdf.getAtt(NCidraw,varid,'node_offset');
   G.title                                                    = netcdf.getAtt(NCidraw,netcdf.getConstant('GLOBAL'),'title' );   
   G.source                                                   = netcdf.getAtt(NCidraw,netcdf.getConstant('GLOBAL'),'source');   

   % copied from GMT manual: node_offset: 
   % 0 for grid line registration (default), 
   % 1 for pixel registration
   
   if G.node_offset==0
   G.lon       = linspace(G.lon_range(1),G.lon_range(2),G.dimension(1));
   G.lat       = linspace(G.lat_range(2),G.lat_range(1),G.dimension(2)); % upside down !
   else
   G.lon       = corner2center1(linspace(G.lon_range(1),G.lon_range(2),G.dimension(1)+1));
   G.lat       = corner2center1(linspace(G.lat_range(2),G.lat_range(1),G.dimension(2)+1)); % upside down !
   end
   
   % add 30 arc-second grid: copied from manual gebco_08.pdf:
   % The complete data sets give global coverage and each file consists 
   % of 21,600 rows x 43,200 columns, resulting in 933,120,000 data points. 
   % The data start at the Northwest corner of the files, i.e. for the global 
   % files, position 89° 59’ 45’’N, 179° 59’ 45’’W and
   % are arranged in latitudinal bands of 360 degrees x 120 points/degree = 43,200
   % values. The data range eastward from 179° 59’ 45’’W to 179° 59’ 45’’E. Thus, the 
   % first band contains 43,200 values for 89° 59’ 45’’N, then followed by a band of 
   % 43,200 values at 89°59’ 15’’N and so on at 30 arc-second latitude intervals down 
   % to 89° 59’ 45’’S. The data values are pixel centre registered i.e. they refer to 
   % data values at the centre of grid cells.   
   
   if OPT.gebco==1
     if rms(degrees2dms(G.lon(1))-[-180 0 0]) > 1e-9;error('error in longitude');end
     if rms(degrees2dms(G.lat(1))-[  90 0 0]) > 1e-9;error('error in latitude' );end
   elseif OPT.gebco==30
     if rms(degrees2dms(G.lon(1))-[-179 59 45.0000]) > 1e-9;error('error in longitude');end
     if rms(degrees2dms(G.lat(1))-[  89 59 45.0000]) > 1e-9;error('error in latitude' );end
   end
   netcdf.close(NCidraw);
   
%% add global attributes to the dataset
   
   OPT.num_bytes = 20000; % pad header, to speed up
   
   mode = netcdf.getConstant('NOCLOBBER'); % do not overwrite existing files
   switch  OPT.netcdfversion
       case 3
           OPT.deflate = false;
           % we stay within 32 bit file size, so 64 bit offset mode not required
       case 4
           OPT.deflate = true;
           mode        = bitor(mode,netcdf.getConstant('NETCDF4'));
   end
   
   NCid     = netcdf.create(ncfile,mode);
   globalID = netcdf.getConstant('NC_GLOBAL');

   netcdf.putAtt(NCid,globalID, 'Conventions',     'CF-1.5 COARDS');
   netcdf.putAtt(NCid,globalID, 'title',           G.title);
   netcdf.putAtt(NCid,globalID, 'institution',     'GEBCO');
   netcdf.putAtt(NCid,globalID, 'source',          G.source);
   netcdf.putAtt(NCid,globalID, 'history',         '$HeadURL$ $Date$ $Author$ $Revision$');
   netcdf.putAtt(NCid,globalID, 'references',      OPT.references);
   netcdf.putAtt(NCid,globalID, 'comment',         '');
   netcdf.putAtt(NCid,globalID, 'email',           'enquiries@bodc.ac.uk');
   netcdf.putAtt(NCid,globalID, 'version',         [OPT.version,', modified to CF by $HeadURL$ $Date$ $Author$ $Revision$']);
   % manually copied from pdf accompanying GEBCO files
   netcdf.putAtt(NCid,globalID, 'terms_of_use',    'Data within this GEBCO grid are subject to copyright and database right restrictions.	Reproduction of the gridded bathymetry sets in derivative form for scientific research, environmental conservation, or education or other non-commercial purposes is authorised without prior permission, providing the source material is properly credited.	The production of these gridded data sets is the result of an international collaboration of numerous scientists and hydrographers who have devoted mich of their time and effort, often on a voluntary basis. The work was stimulated by a wish to create an authoritative, high-quality bathymetry of the world''s oceans for the benefit of us all. Therefore we ask you to contact us first, should you wish to pass on the data to third parties of use the data for commercial purposes.	In the first instance, please contact the British Oceanographic Data Centre (BODC) <enquiries@bodc.ac.uk> and include a clear statement of the purpose for which the material will be used and the manner in which it will be reproduced. 	In the case of commercial activities, a contribution may be requested for the further improvement of GEBCO''s datasets.');
   netcdf.putAtt(NCid,globalID, 'attribution',     ['Any material reproduced from GEBCO gridded data sets should be accompanied by appropriate attribution to the source of the material. The version number of the grid is given in the header information within the grid data file.	For the GEBCO bathymetry data, this should be of the form (including the appropriate version number): Reference for this grid: ''',OPT.version,'''	For non-GEBCO material, such as land data from the SRTM30 dataset, the ackowledgement should include reference to the original source material.']);
   netcdf.putAtt(NCid,globalID, 'disclaimer',      'THIS GEBCO GRID IS NOT TO BE USED FOR NAVIGATION OR FOR ANY OTHER PURPOSE INVOLVING SAFETY AT SEA.	Data use in the creation of the GEBCO 08 Grid and GEBCO One Minute Grid grid have been obtained from sources believed to be relaible but their accuracy and completeness cannot be guaranteed. Whilst every effort has been made to ensure their reliability within the limits of present knowlegde, no responsibility can be acccepted by those involved in their compilation or publication for  any consequential loss or damage arasing from their use.	This GEBCO grid is essentially a deep water product and does not include detailed bathymetry for shallow shelf waters. Even to the present day, most areas of the world''s oceans have not been fully surveyed and, for the most part, bathymetric mapping is an interpolation based on random tracklines of data from many different sources. The quality and coverage of data from these sources is highly variable. Although the GEBCO 08 Grid and GEBCO One Minute Grid grid are presented at 30 arc-second and one arc-minute intervals of latitude and longitude respectively, this does not imply that knowlegde is available on sea floor depth at this resolution - the depth in most 30 arc-second squares of the world has yet to be measured.');

%% specify dimensions

   netcdf.defDim(NCid,'lon',G.dimension(1));
   netcdf.defDim(NCid,'lat',G.dimension(2));

%% define coordinate variables
  
   nc.Name             = 'lon';
   nc.Nctype           = nc_double; % vector so double does not occupy much space, and accuracy is required
   nc.Dimension        = {'lon'};
   nc.Attribute(    1) = struct('Name', 'standard_name'  ,'Value', 'longitude');
   nc.Attribute(end+1) = struct('Name', 'long_name'      ,'Value', 'Uniformly spaced longitudes (-180E - +180E)');
   nc.Attribute(end+1) = struct('Name', 'units'          ,'Value', 'degrees_east');
   nc.Attribute(end+1) = struct('Name', 'axis'           ,'Value', 'X');
   nc.Attribute(end+1) = struct('Name', 'spacing'        ,'Value', G.spacing(1));
   nc.Attribute(end+1) = struct('Name', 'actual_range'   ,'Value', [G.lon(1) G.lon(end)]);
   netcdf_addvar(NCid, nc);
   
   nc.Name             = 'lat';
   nc.Nctype           = nc_double; % vector so double does not occupy much space, and accuracy is required
   nc.Dimension        = {'lat'};
   nc.Attribute(    1) = struct('Name', 'standard_name'  ,'Value', 'latitude');
   nc.Attribute(end+1) = struct('Name', 'long_name'      ,'Value', 'Uniformly spaced latitudes (-90E - +90E)');
   nc.Attribute(end+1) = struct('Name', 'units'          ,'Value', 'degrees_north');
   nc.Attribute(end+1) = struct('Name', 'axis'           ,'Value', 'Y');
   nc.Attribute(end+1) = struct('Name', 'spacing'        ,'Value', G.spacing(2));
   nc.Attribute(end+1) = struct('Name', 'actual_range'   ,'Value', [G.lat(1) G.lat(end)]);
   netcdf_addvar(NCid, nc);
   
   nc.Name             = 'topo';
   nc.Nctype           = 'short'; % int16
   nc.Dimension        = {'lon','lat'}; % flipped compared to snctools and ncBrowse
   nc.Attribute(    1) = struct('Name', 'standard_name'  ,'Value', 'surface_altitude');
   nc.Attribute(end+1) = struct('Name', 'long_name'      ,'Value', 'topography');
   nc.Attribute(end+1) = struct('Name', 'units'          ,'Value', 'meters');
   nc.Attribute(end+1) = struct('Name', '_CoordinateAxes','Value', 'lat lon'); % to make it match WHOI DEMs, not in any known standard
   nc.Attribute(end+1) = struct('Name', 'positive'       ,'Value', 'up');
   nc.Attribute(end+1) = struct('Name', 'axis'           ,'Value', 'Z');
   nc.Attribute(end+1) = struct('Name', 'actual_range'   ,'Value', [nan nan]);
   netcdf_addvar(NCid, nc);
   varid  = netcdf.inqVarID(NCid,'topo');
   if OPT.deflate
      netcdf.defVarDeflate(NCid,varid,false,true,2);
   end
   
   if OPT.gebco==30
  % map int16 to int8 to save almost 1GB of disk space
  % the GEBCO codes are in too wide an integer range
  % for their purpose
   SID.mapped_values   = [1    2     3   4   5   6    7    8    9   10]; % choose 0/1-based?
   SID.original_values = [0 9999 -8888 110 120 130 1000 1100 1200 1300];
   SID.values          = '0 9999 -8888 110 120 130 1000 1100 1200 1300';
   SID.short_labels    = ['gravimetry + sounding guidance;',...
                          'gravimetry + sounding constraint; ',...
                          'land (SRTM30);',...
                          'Caspian Sea, sounding;',...
                          'Black Sea, sounding;',...
                          'Weddell Sea, sounding;',...
                          'IBCAO;',...
                          'Caspian Sea, interpolated;',...
                          'Black Sea, interpolated;',...
                          'Weddell Sea, interpolated'];
   SID.long_labels     = ['interpolated ship-track soundings guided by gravimetry data (0);',...
                          'interpolated ship-track soundings guided by gravimetry data constrained by bathymetric soundings (9999);',...
                          'land from SRTM30 (-8888);',...
                          'the Caspian Sea (110);',...
                          'the Black Sea (120);',...
                          'the Weddell Sea (130);',...
                          'IBCAO (1000);',...
                          'the Caspian Sea without sounding (1100);',...
                          'the Black Sea without sounding (1200);',...
                          'the Weddell Sea without sounding (1300);'];

   nc.Name             = 'SID';
   nc.Nctype           = 'char'; % int8
   nc.Dimension        = {'lon','lat'}; % flipped compared to snctools and ncBrowse
   nc.Attribute(    1) = struct('Name', 'standard_name'    ,'Value', '');
   nc.Attribute(end+1) = struct('Name', 'long_name'        ,'Value', '');
   nc.Attribute(end+1) = struct('Name', 'units'            ,'Value', '');
   nc.Attribute(end+1) = struct('Name', '_CoordinateAxes'  ,'Value', 'lat lon');  % to make it match WHOI DEMs, not in any known standard
   % note: mapped to int18
   nc.Attribute(end+1) = struct('Name', 'flag_values'      ,'Value', SID.mapped_values);  % CF-1.5 proposal
   nc.Attribute(end+1) = struct('Name', 'flag_meanings'    ,'Value', SID.values);         % CF-1.5 proposal
   nc.Attribute(end+1) = struct('Name', 'flag_short_labels','Value', SID.short_labels);   % not a CF standard, for plotting leged
   nc.Attribute(end+1) = struct('Name', 'flag_long_labels' ,'Value', SID.long_labels);    % not a CF standard, full text from GEBCO manual
   nc.Attribute(end+1) = struct('Name', 'flag_comments'    ,'Value', 'For detailed SID descriptiond see GEBCO documentation');
   netcdf_addvar(NCid, nc);
   end
   
%% Expand NC file

   netcdf.endDef(NCid,OPT.num_bytes,4,0,4); 

%% add coordinate data

   varid = netcdf.inqVarID(NCid,'lat');netcdf.putVar(NCid,varid,G.lat');
   varid = netcdf.inqVarID(NCid,'lon');netcdf.putVar(NCid,varid,G.lon');

%% close NC file

   netcdf.close(NCid)
   
   nc_dump(ncfile); % snctools: see http://mexcdf.sourceforge.net/

%%  add altitude data: now we are going to add the data
%   chunkwise to circumvent memory issues
%   This code cell can be run stand-alone from the
%   above steps that created the netcd file: this code cell 
%   cell just overwrites the contents of variable 'topo'
%                          N
%        +------------------------------------+  o-> cx = dimension(1)
%        |                                    |  |                  
% ic=1   |                                    |  v cy = chunksize
%        |                                    |  
%        +------------------------------------+  
%        |                                    |  o-> cx = dimension(1)
% ic=2   |W                                  E|  |                  
%        |                                    |  v cy = chunksize
%        +------------------------------------+
%        :                                    :
%        +------------------------------------+
% ic=nc  |cylast<=cy, chunksize of last chunk | 
%        +------------------------------------+
% |                         S
% v nc = number of chunks

NCid   = netcdf.open(ncfile,netcdf.getConstant('WRITE'));

if OPT.gebco==30
   nvar=2;
else
   nvar=1;
end

for ivar=1:nvar

   if ivar==1
   nc_from    = ncfileraw
   NCidraw    = netcdf.open(ncfileraw,netcdf.getConstant('NC_NOWRITE'));
   varnameraw = 'z';
   varname    = 'topo';
   varidraw   = netcdf.inqVarID(NCidraw,varnameraw);
   else
   nc_from    = ncfileSID;
   NCidraw    = netcdf.open(ncfileSID,netcdf.getConstant('NC_NOWRITE'));
   varnameraw = 'z';
   varname    = 'SID';
   varidraw   = netcdf.inqVarID(NCidraw,varnameraw);
   end

   varid   = netcdf.inqVarID(NCid,varname);

  [varname,xtype,dimids,natts] = netcdf.inqVar(NCid,varid);
  [~,nx]  = netcdf.inqDim(NCid,dimids(1));
  [~,ny]  = netcdf.inqDim(NCid,dimids(2));

   cx     = nx;            % chunk size in lon direction: entire zonal line
   cy     = OPT.chunksize; % chunk size in lat, depends on PC memory limitations
   nc     = ceil(ny/cy);   % number of chuncks (note, last one can have odd chunk size)
   cylast =  mod(ny,cy);   % chunk size in lat of last chunk

   zmin = +Inf; % startvalues for keeping ...
   zmax = -Inf; % ... track of actual_range
   
   for ic=1:nc
       
      start0 = (ic-1)*cx*cy;
      start  = [0 (ic-1)*cy];
      if ic==nc
         count0 =  cx*cylast;
         count  = [cx cylast];
      else
         count0 =  cx*cy;
         count  = [cx cy];
      end
     %G.z   = reshape(nc_varget    (nc_from,'z'     ,start0,count0),[cx count0/cx]); % snctools: see http://mexcdf.sourceforge.net/
      G.z   = reshape(netcdf.getVar(NCidraw,varidraw,start0,count0),[cx count0/cx]);
      
      % map int16 to int8
      if strcmpi(varname,'SID')
         G.z = int16(G.z);
         for imap=1:length(SID.original_values)
            mask = (G.z==SID.original_values(imap));
            G.z(mask) = SID.mapped_values(imap); % note: in-variable mapping works because SID.mapped_values are not elements of SID.values
         end
         G.z = char(G.z); % matlab netcdf library needs this to save it as 8-bit variable
      else
         G.z = int16(G.z);
      end
      
      netcdf.putVar(NCid,varid,start,count,G.z);
      
      zmin = min(zmin,min(G.z(:)));
      zmax = max(zmax,max(G.z(:)));
      
      disp(['done: ',num2str(ic/nc*100,'%.2f'),'%'])
   
   end
   
   netcdf.putAtt(NCid,varid,'actual_range',[zmin zmax]);

netcdf.close(NCidraw);

end   

netcdf.close(NCid);

%%
fid = fopen(strrep(ncfile,'nc','cdl'))
nc_dump(ncfile,fid); % snctools: see http://mexcdf.sourceforge.net/
fclose(fid)
