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