In this example we will show how to benefit from OPeNDAP subsetting. We will request subsets from worldwide digital elevevation maps (DEMs) that are too big to request as a whole: resulting in approx. 500 MB (over 10,000 x 20,000 points) or 2 GB (over 20,000 x 40,000) for respectively a 1 minute or 30 second grids.

Add snctools http://mexcdf.sourceforge.net/, shipped with OpenEarthTools

run('...\openearthtools\matlab\oetsettings.m')

Worl digital elevation maps: Define data on an opendap server. These are 6 free datasets, which we will compare for the Southern North Sea. Please compare them to non-open datasets such as GEBCO. For details on the processing of the 6 DEM's, please refer to the meta-data in the THREDDS OPeNDPA server provided by USGS/WHOI.

url_grid{1} = 'http://geoport.whoi.edu/thredds/dodsC/bathy/etopo1_bed_g2';
url_grid{2} = 'http://geoport.whoi.edu/thredds/dodsC/bathy/etopo2_v2c.nc';
url_grid{3} = 'http://geoport.whoi.edu/thredds/dodsC/bathy/srtm30plus_v1.nc';
url_grid{4} = 'http://geoport.whoi.edu/thredds/dodsC/bathy/srtm30plus_v6';
url_grid{5} = 'http://geoport.whoi.edu/thredds/dodsC/bathy/smith_sandwell_v9.1.nc';
url_grid{6} = 'http://geoport.whoi.edu/thredds/dodsC/bathy/smith_sandwell_v11';

url_line    = 'http://opendap.deltares.nl/thredds/dodsC/opendap/noaa/gshhs/gshhs_i.nc';

World vector shoreline: Get line data: 1D vectors are small, so we can get all data

nc_dump(url_line)
L.lon    = ncread(url_line,'lon');
L.lat    = ncread(url_line,'lat');

Define bounding box

boundingbox.lon = [ 0 10];
boundingbox.lat = [50 55];

for i=1:length(url_grid)

   ncfile = url_grid{i}

  • Get full lat,lon vectors: 1D vectors are small, so we can get all data

    nc_dump(ncfile)
   G.lon    = ncread(ncfile,'lon' ); % 1D
   G.lat    = ncread(ncfile,'lat' ); % 1D

  • Find the subset-indices within the bounding box

    ilon     = find(G.lon > boundingbox.lon(1) & G.lon < boundingbox.lon(2));
ilat     = find(G.lat > boundingbox.lat(1) & G.lat < boundingbox.lat(2));

  • ranslate subset-indices to netCDF argument: [start,count,stride]

    stride   = [1 1]; % additionally specify a stride when the subset is still too big
start    = [min(ilon) min(ilat)]; % netCDF is 0-based, whereas matlab is 1-based
count    = ceil([length(ilon) length(ilat)]./stride); % use ceil to cover at least bounding box area

  • Request data subset

    G.lat    = ncread(ncfile,'lat' ,start(2),count(2),stride(2)); % 1D
G.lon    = ncread(ncfile,'lon' ,start(1),count(1),stride(1)); % 1D
G.topo   = ncread(ncfile,'topo',start(:),count(:),stride(:)); % 2D
G.title  = ncreadatt(ncfile,'/','title')

  • Plot data subset

    figure(i)
   pcolorcorcen(G.lon,G.lat,G.topo')
   hold on
   plot(L.lon,L.lat,'k')
   axis([boundingbox.lon boundingbox.lat])
   tickmap('ll','texttype','text','dellast',1)
   axislat % sets aspect ratio
   grid on
   clim ([-50 150])
   title(mktex(G.title))
   colorbarwithvtext('z [m]')
   print2screensize(mkvar(G.title))

end

Download the code of this Matlab example (repos,manual download). We also have a script that rewrites the GEBCO 1D vector matrix to a proper 2D lat,lon matrix with the same lay-out as the above OPeNDAP sets. (repos,manual download).

non-open GEBCO data source for comparison [The GEBCO_08 Grid version 20091120, http://www.gebco.net]
non-open GEBCO data for comparison [The GEBCO One Minute Grid, version 2.0, http://www.gebco.net]
non-open GEBCO data for comparison [The GEBCO_08 Grid version 20091120, http://www.gebco.net]

Disclaimer