Page History

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    = nc_varget(url_line,'lon');
L.lat    = nc_varget(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

  Code Block
  nc_dump(ncfile) G.lon = nc_varget(ncfile,'lon' ); % 1D G.lat = nc_varget(ncfile,'lat' ); % 1D

• Find the subset-indices within the bounding box

  Code Block
  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]

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

• Request data subset

  Code Block
  G.lat = nc_varget(ncfile,'lat' ,start(1),count(1),stride(1)); % 1D G.lon = nc_varget(ncfile,'lon' ,start(2),count(2),stride(2)); % 1D G.topo = nc_varget(ncfile,'topo',start(:),count(:),stride(:)); % 2D G.title = nc_attget(ncfile,nc_global,'title')

• Plot data subset

  Code Block
  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).