Recipe: Python program as OPeNDAP client

From GEOS-5
Jump to: navigation, search

Back to G5NR Data Access Guide.

Problem

By accessing the collection inst01hr_3d_T_Cv via the OPeNDAP server

http://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km

we want to read the surface temperature data inside the box bound by latitudes 25oN, 50oN and longitudes -130oW, -65oW for 2006/Sep/18, 9z, compute its min/max and plot the temperature data.

Solution

Here we use the netcdf4-python package. If the command

> python -c "import netCDF4"

does not return any error, we have the package installed.

Code

This code accesses the collection inst01hr_3d_T_Cv from the OPeNDAP server, reads a subset of the temperature data (all levels inside the bounding box specified above) and computes its max/min. It then plots the data at the surface (level=72).

NOTE:

  1. Instead of reading a downloaded NetCDF-4 file, we read an OPeNDAP URL.
  2. While in the downloaded file, the temperature variable appears in the uppercase (T), on the OPeNDAP server, this variable is in lowercase.
  3. Via the OPeNDAP URL, we now have access to all times for which data exists. The hourly inst files are available starting at 2005/15/15, 2200z. Our desired time, 2006/09/18, 0900z is then the 11772th step.
 1 #!/usr/bin/env python
 2 
 3 import sys
 4 import numpy as np
 5 import netCDF4 as nc4
 6 import matplotlib.pyplot as plt
 7 
 8 from mpl_toolkits.basemap import Basemap
 9 
10 rootgrp = nc4.Dataset('http://opendap.nccs.nasa.gov:9090/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv', 'r')
11 print "rootgrp.variables['t'].shape", rootgrp.variables['t'].shape
12 
13 # read subset of air temperature
14 print 'Reading T (subset)...',; sys.stdout.flush()
15 T = rootgrp.variables['t'][11771,:,229:280,99:230]
16 print 'done.'; sys.stdout.flush()
17 print 'T.shape:', T.shape
18 
19 # max/min
20 print 'max(T): %.4f' % np.max(T)
21 print 'min(T): %.4f' % np.min(T)
22 
23 # set up cylindrical map
24 m = Basemap(
25     projection='cyl',
26     llcrnrlat=25, urcrnrlat=50,
27     llcrnrlon=-130, urcrnrlon=-65,
28     resolution='c'
29     )
30 m.drawcoastlines(linewidth=1.5)
31 m.drawmapboundary()
32  
33 # plot contour
34 level = 71 # surface
35 X = np.arange(-130.0, -64.99, .5) # -65 is the last element
36 Y = np.arange(25.0, 50.01, .5) # 50 is the last element
37 cp = plt.contour(X, Y, T[level,:,:], 20, zorder=2)
38 plt.clabel(cp, inline=1, fontsize=9)
39 plt.title('Air temperature at the surface')
40 plt.show()

Output

Running this python script, we get the text output

rootgrp.variables['t'].shape (18288, 72, 361, 720)
Reading T (subset)... done.
T.shape: (72, 51, 131)
max(T): 305.6048
min(T): 191.6956

and the plot

BoundingBoxSurfaceTemp.png

Discussions

See Also

  1. File Spec: File:G5NR-Ganymed-7km FileSpec-ON6-V1.0.pdf
  2. Recipe: Python program to read data from downloaded file


No Warranty

Copyright