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.

#!/usr/bin/env python

#OPENDAP_SUBSETTING_WITH_PYTHON_TUTORIAL how to benefit from OPeNDAP subsetting in python
# This tutorial is also available for Matlab
#See also:

# $Id:  $
# $Date:  $
# $Author:  $
# $Revision:  $
# $HeadURL: $
# $Keywords: $

# This document is also posted on a wiki:

# Read data from an opendap server
import urllib
import numpy as np
import netCDF4
import pydap
import matplotlib
import matplotlib.pyplot as plt
import pylab

# from opendap import opendap # OpenEarthTools module, see above that makes pypdap quack like netCDF4

# Define data on an opendap server
# for converting non-open access GECBO grids to same netCDF structure: see
gridurls =  ['',

lineurl  = '';

# Get line data: 1D vectors are small, so we can get all data
linedata = netCDF4.Dataset(lineurl) # opendap(url_line) # when netCDF4 was not compiled with OPeNDAP

lines = dict(

# Define bounding box
BB = dict(
   lon=[ 0, 10],
   lat=[50, 55]

for i,ncfile in enumerate(gridurls):

   print ncfile

   # Get grid data
   grid   = netCDF4.Dataset(ncfile) # opendap(ncfile) # when netCDF4 was not compiled with OPeNDAP

   # Get all data from lat and lon, but just a pointer to the topography
   # Because getting it all takes a bit too long
   lat = grid.variables['lat'][:]
   lon = grid.variables['lon'][:]

   # nonzero returns a tuple of idx per dimension
   # we're unpacking the tuple here so we can lookup max and min
   (latidx,) = np.logical_and(lat >= BB['lat'][0], lat < BB['lat'][1]).nonzero()
   (lonidx,) = np.logical_and(lon >= BB['lon'][0], lon < BB['lon'][1]).nonzero()

   assert lat.shape[0] == grid.variables['topo'].shape[0], 'We assume first dim is lat here'
   assert lon.shape[0] == grid.variables['topo'].shape[1], 'We assume 2nd dim is lon here'

   # get rid of the non used lat/lon now
   lat = lat[latidx]
   lon = lon[lonidx]
   # Get the topography data with the new indices
   topo = grid.variables['topo'][latidx.min():latidx.max(),
   title = grid.title # remember so we can close the file before plotting

   Lon,Lat = np.meshgrid(lon, lat)
   # make a pcolormesh (the fastest way to plot simple grids)
   mesh = plt.pcolormesh(Lon,Lat,topo)
   # some customizations (called on the axes and figure)
   mesh.axes.set_aspect(1/np.cos(np.pi*np.mean(BB['lat'])/180.)) # set aspect to match meters
   mesh.axes.set_xlabel('lon [deg]')
   mesh.axes.set_ylabel('lat [deg]')
   mesh.figure.colorbar(mesh) # use the mesh to make a colorbar
   # Save the figure
   mesh.figure.savefig('%s.png' % title)

Download the code of this python example (repos,manual download)

See also: OPeNDAP access with python, OPeNDAP subsetting with Matlab,OPeNDAP subsetting with R, PostgreSQL access with python