%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: GEBCO2CF.m 5150 2011-11-01 20:48:31Z boer_g $ % $Date: 2011-11-01 21:48:31 +0100 (di, 01 nov 2011) $ % $Author: boer_g $ % $Revision: 5150 $ % $HeadURL: https://svn.oss.deltares.nl/repos/openearthrawdata/trunk/gebco/scripts/GEBCO2CF.m $ % $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: https://svn.oss.deltares.nl/repos/openearthrawdata/trunk/gebco/scripts/GEBCO2CF.m $ $Date: 2011-11-01 21:48:31 +0100 (di, 01 nov 2011) $ $Author: boer_g $ $Revision: 5150 $'); 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: https://svn.oss.deltares.nl/repos/openearthrawdata/trunk/gebco/scripts/GEBCO2CF.m $ $Date: 2011-11-01 21:48:31 +0100 (di, 01 nov 2011) $ $Author: boer_g $ $Revision: 5150 $']); % 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) 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)