Child pages
  • OPeNDAP subsetting with matlab

Versions Compared


  • This line was added.
  • This line was removed.
  • Formatting was changed.


Add snctools, shipped with OpenEarthTools

Code Block


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.

Code Block

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

Code Block

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

Define bounding box

Code Block

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

    Code Block
       G.lon    = nc_vargetncread(ncfile,'lon' ); % 1D    = nc_vargetncread(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( > & <;
  • 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(ilatilon)-1 min(ilonilat)-1]; % subtract one as netCDF is 0-based, whereas matlab is 1-basesbased
       count    = ceil([length(ilatilon) length(ilonilat)]./stride); % use ceil to cover at least bounding box area
  • Request data subset

    Code Block    = nc_vargetncread(ncfile,'lat' ,start(12),count(12),stride(12)); % 1D
       G.lon    = nc_vargetncread(ncfile,'lon' ,start(21),count(21),stride(21)); % 1D
       G.topo   = nc_vargetncread(ncfile,'topo',start(:),count(:),stride(:)); % 2D
       G.title  = nc_attgetncreadatt(ncfile,nc_global'/','title')
  • Plot data subset

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


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