Python

# Some general python stuff: Data Import
code import numpy as N import Scientific.IO.NetCDF as S
 * 1) NCEP reanalysis data, from Johnny Lin's AMS short course http://www.johnny-lin.com/ams2011/sc/gen/files/
 * 2) Here's the key reading lines. Scientific.IO.NetCDF module must be in the [|CDAT] distribution of Python.
 * 1) Module imports:

fileobj = S.NetCDFFile('air.mon.mean.nc', mode='r') T = fileobj.variables['air'].getValue lat = fileobj.variables['lat'].getValue lat = fileobj.variables['lon'].getValue fileobj.close
 * 1) netCDF Reading of air temperature: open, read, close

code
 * 1) ??? Do we need to unpack the T array (short integers), or did it happen automatically?
 * Q: Anybody know how to get at the datafile with OpenDAP?
 * 1) Is it as simple as replacing 'air.mon.mean.nc'
 * 2) with 'http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface_flux/variable.averagingperiod.statistic.nc'
 * 3) or do more modules need to be imported? (antigravity)
 * A: Should be simple, try importing pydap, get it from http://pydap.org/.

To download files from within a python program, do something like: (works just like wget or curl in Linux/UNIX) code format="python" from urllib import urlretrieve urlretrieve('ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/surface/uwnd.mon.ltm.nc','uwnd.mon.ltm.nc')
 * 1) Module import:
 * 1) Download file into current working directory:
 * 2) syntax is urlretrieve('url','filename')
 * 1) syntax is urlretrieve('url','filename')

code To access a data file through OpenDAP, use pydap module: code format="python" from pydap.client import open_url ncfile=open_url('http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface_gauss/prate.mon.ltm.nc')
 * 1) Module import

code To see contents of an object, just type the object name: code >>> ncfile {'lat': , 'lon': , 'time': , 'prate': {'prate': , 'time': , 'lat': , 'lon': }}

code Notice that ncfile is a dictionary object that contains {'lat', 'lon', 'time', 'prate'{'prate', 'time', 'lat', 'lon'}}, where the 4th element 'prate' is a dictionary object itself. This is because variables contain dimensions as their attributes. To get the precipitation rate field, we type: code format="python" >>> a=ncfile['prate']['prate'] # This line just accesses the object descriptor, not the data itself >>> a  >>> a.shape (12, 94, 192) >>> a=ncfile['prate']['prate'][:,:,:] # append [:,:,:] to actually download data into temp memory

code

Joni Lum code code format="python" import numpy as N # Numeric Python from pupynere import netcdf_file # NetCDF libraries from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as PL # Plotting libraries for varstr in ['vwnd','air','precip']: ncfile=netcdf_file(varstr+'.mon.ltm.nc','r') print ncfile.title if varstr=='vwnd': lat=ncfile.variables['lat'].getValue # Get latitude lon=ncfile.variables['lon'].getValue # Get longitude time=ncfile.variables['time'].getValue # Get time time=time/24 # in days v=ncfile.variables[varstr].getValue # Get zonal wind v_x=N.mean(v,axis=2) # Zonal wind mean v_90=v[:,:,N.argmin(abs(lon-90))] # 90W img=PL.contourf(time,lat,N.transpose(v_x)) PL.contour(time,lat,N.transpose(v_x),colors='black') PL.colorbar(img) PL.title('Zonally Averaged Meridonal Velocity {m/s}') PL.savefig('v_x.png',format='png') PL.clf PL.figure img=PL.contourf(time,lat,N.transpose(v_90)) PL.contour(time,lat,N.transpose(v_90),colors='black') PL.colorbar(img) PL.title('Meridional Velocity at 90W {m/s}') PL.savefig('V_90.png',format='png') elif varstr=='air': lat=ncfile.variables['lat'].getValue # Get latitude T=ncfile.variables[varstr].getValue # Get zonal wind dlon=N.cos(N.deg2rad(lat)) # Compute degree longitude size dlon=N.reshape(dlon,(1,73)) # Reshape to 3-D array dlon=N.tile(dlon,(12,1)) # Tile the time column T_x=N.mean(T,axis=2) # Zonal Mean T_xy=N.mean(T_x*dlon,axis=1)/N.mean(dlon) # Meridional mean T_xyt=N.mean(T_xy) # Time space mean T_stdev=N.mean((T-T_xyt)**2.,axis=2) # STD part 1 T_stdev=N.sqrt(N.mean(T_stdev*dlon)/N.mean(dlon)) # STD part 2 PL.clf img=PL.plot(T_xy) PL.title('Average Air Temp Over 12 Months') PL.savefig('Averge_Air_T.png',format='png') elif varstr=='precip': lat=ncfile.variables['lat'].getValue # Get latitude P=ncfile.variables[varstr].getValue # Get precip. P_t=N.mean(P,axis=0) # Take time mean P_t=N.reshape(P_t,(1,72,144)) # Reshape to 3-D P_t=N.tile(P_t,(12,1,1)) # Tile the time column P_stdev=100.*N.sqrt(N.mean((P-P_t)**2.,axis=0))/P_t # Standard deviation in % P_stdev=N.mean(P_stdev,axis=0) # STD time mean PL.clf img=PL.contourf(lon,lat,P_stdev) PL.colorbar(img,orientation='horizontal') img=PL.contour(lon,lat,P_stdev,colors='black') PL.title('Precipitation STD [%]') img=Basemap(llcrnrlon=0,llcrnrlat=-88.75,urcrnrlon=357.5,urcrnrlat=88.75) img.drawcoastlines PL.savefig('precip.png',format='png') P_x=N.mean(P,axis=2) # Zonal Mean P_90=P[:,:,N.argmin(abs(lon-90))] # 90W PL.clf img=PL.contourf(time,lat,N.transpose(P_x)) PL.contour(time,lat,N.transpose(P_x),colors='black') PL.colorbar(img) PL.title('Zonally Averaged Precipitation [mm/ day]') PL.savefig('precip_x.png',format='png') PL.clf img=PL.contourf(time,lat,N.transpose(P_90)) PL.colorbar(img) PL.contour(time,lat,N.transpose(P_90),colors='black') PL.title('Precipitation at 90W [mm/day]') PL.savefig('P_90.png',dpi=100) code
 * 1) !/usr/bin/env python
 * 2) MPO581 Homework #1
 * 3) First load modules we need:
 * 1) First load modules we need:
 * 1) First load modules we need:
 * 1) Load data files:
 * 1) Load data files:
 * 1) Figure 1
 * 1) Figure 1
 * 1) Figure 2
 * 1) Figure 2
 * 1) Figure Average Air Temp over 12 months
 * 1) Figure Average Air Temp over 12 months
 * 1) Figure 5
 * 1) Figure 5
 * 1) Figure 3
 * 1) Figure 3
 * 1) Figure 4
 * 1) Figure 4