Python

How to use the OPeNDAP server in Python.

The following explains the steps of loading data served by OPeNDAP using Python.

Download and Install Python
with required external modules:

Or download and install Anaconda (Python 2.7)
with update/install external modules in command prompt:

Open Python:

In [1]:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:		 from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset, date2index
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime

# Contours ssh for specified time and space
# Find the index for date Oct 3, 2004.
date = datetime(2004,10,3,0)
# Open dataset
dataset = Dataset('http://apdrc.soest.hawaii.edu/dods/public_data/NLOM/nlom_ssh')
timevar = dataset.variables['time']
timeindex = date2index(date,timevar) # find time index for desired date.
# Read lats and lons (representing centers of grid boxes).
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]
# Find the indexes for 165W-153W, 18N-24N
lat_bnds = [ 18 , 24 ]
lon_bnds = [ -165+360 , -154+360 ]
lat_inds = np.where((lats >= lat_bnds[0]) & (lats <= lat_bnds[1]))[0]
lon_inds = np.where((lons >= lon_bnds[0]) & (lons <= lon_bnds[1]))[0]
# Load the NLOM sea surface height from the above information
ssh = dataset.variables['ssh'][timeindex,lat_inds,lon_inds].squeeze()
# Create figure, axes instances.
fig = plt.figure()
ax = fig.add_axes([0.05,0.05,0.9,0.9])
# Create Basemap instance.
m = Basemap(llcrnrlon=-165.,llcrnrlat=18.,urcrnrlon=-154.,urcrnrlat=24.,projection='mill')
# Plot land outlines
m.drawcoastlines()
# Use subset of lats and lons
sublons, sublats = np.meshgrid(lons[lon_inds]-360,lats[lat_inds])
# Contour sea surface height
im1 = m.contourf(sublons, sublats, ssh, latlon=True)
cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
ax.set_title('NLOM 1/16 degree Sea Surface Height [cm] October 3, 2004')
plt.show()

Using pcolormesh to plot:

In [1]:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
 from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset, date2index
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime

date = datetime(2000,1,1,0)
# open dataset.
dataset = Dataset('http://apdrc.soest.hawaii.edu/dods/public_data/SCUD/current')
timevar = dataset.variables['time']
timeindex = date2index(date,timevar)
# read u current
uc = dataset.variables['uc'][timeindex,:].squeeze()
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]
lons, lats = np.meshgrid(lons,lats)
# create figure, axes instances.
fig = plt.figure()
ax = fig.add_axes([0.05,0.05,0.9,0.9])
# create Basemap instance.
m = Basemap(projection='mill',lon_0=180)
m.drawcoastlines()
im1 = m.pcolormesh(lons,lats,uc,latlon=True)
cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
ax.set_title('SCUD u current for %s'%date)
plt.show()

Using Argo from ERDDAP (through APDRC Dapper):
Reference: http://matplotlib.org/basemap/users/examples.html

In[1]:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:		 from netCDF4 import Dataset, num2date
import time, calendar, datetime, numpy
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import urllib, os

# data downloaded from the form at
# http://coastwatch.pfeg.noaa.gov/erddap/tabledap/apdrcArgoAll.html
filename, headers = urllib.urlretrieve('http://coastwatch.pfeg.noaa.gov/erddap/tabledap/apdrcArgoAll.nc? longitude,latitude,time&longitude>=0&longitude<=360&latitude>=-90&latitude<=90&time>=2017-01-01&time<=2017-01-02&distinct()')
dset = Dataset(filename)
lats = dset.variables['latitude'][:]
lons = dset.variables['longitude'][:]
time = dset.variables['time']
times = time[:]
t1 = times.min(); t2 = times.max()
date1 = num2date(t1, units=time.units)
date2 = num2date(t2, units=time.units)
dset.close()
os.remove(filename)
# draw map with markers for float locations
m = Basemap(projection='hammer',lon_0=180)
x, y = m(lons,lats)
m.drawmapboundary(fill_color='#99ffff')
m.fillcontinents(color='#cc9966',lake_color='#99ffff')
m.scatter(x,y,3,marker='o',color='k')
plt.title('Locations of %s ARGO floats active between %s and %s' %\
(len(lats),date1,date2),fontsize=12)
plt.show()

Using FNMOC Indian Ocean from ERDDAP (through APDRC FNMOC):

In[1]:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
      ...:
 from netCDF4 import Dataset, num2date
import time, calendar, datetime, numpy
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import urllib, os

filename, headers = urllib.urlretrieve('http://apdrc.soest.hawaii.edu/erddap/tabledap/fnmoc_temperature_indian.nc?longitude,latitude,time,temperature&longitude>=0&longitude<=360&latitude>=-90&latitude<=90&time>=2014-01-01&time<=2014-01-31&distinct()')
dset = Dataset(filename)
lats = dset.variables['latitude'][:]
lons = dset.variables['longitude'][:]
time = dset.variables['time']
temps = dset.variables['temperature'][:]
times = time[:]
t1 = times.min(); t2 = times.max()
date1 = num2date(t1, units=time.units)
date2 = num2date(t2, units=time.units)
dset.close()
os.remove(filename)
# draw map with markers for float locations
m = Basemap(llcrnrlon=17.5,llcrnrlat=-59.17,urcrnrlon=133.37,urcrnrlat=29.88,projection='mill')
x, y = m(lons,lats)
m.drawmapboundary(fill_color='#99ffff')
m.fillcontinents(color='#cc9966',lake_color='#99ffff')
cm = plt.cm.get_cmap('jet')
temp1 = temps.min(); temp2 = temps.max()
sc = m.scatter(x,y,c=temps, vmin=temp1, vmax=temp2, s=5, cmap=cm)
plt.colorbar(sc)
plt.title('FNMOC Indian Ocean temperature active between %s and %s' %\
(date1,date2),fontsize=12)
plt.show()