Child pages
  • OPeNDAP subsetting with matlab
Skip to end of metadata
Go to start of metadata

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, shipped with OpenEarthTools


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} = '';
url_grid{2} = '';
url_grid{3} = '';
url_grid{4} = '';
url_grid{5} = '';
url_grid{6} = '';

url_line    = '';

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

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

Define bounding box

boundingbox.lon = [ 0 10]; = [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

       G.lon    = ncread(ncfile,'lon' ); % 1D    = 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( > & <;
  • 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    = 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

       hold on
       axislat % sets aspect ratio
       grid on
       clim ([-50 150])
       colorbarwithvtext('z [m]')

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