# plot 3D fields import datetime as dt # Python standard library datetime module import numpy as np from netCDF4 import Dataset # http://code.google.com/p/netcdf4-python/ import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap, addcyclic, shiftgrid from matplotlib.backends.backend_pdf import PdfPages from mpl_toolkits.axes_grid1 import make_axes_locatable def ncdump(nc_fid, verb=False): ''' ncdump outputs dimensions, variables and their attribute information. The information is similar to that of NCAR's ncdump utility. ncdump requires a valid instance of Dataset. Parameters ---------- nc_fid : netCDF4.Dataset A netCDF4 dateset object verb : Boolean whether or not nc_attrs, nc_dims, and nc_vars are printed Returns ------- nc_attrs : list A Python list of the NetCDF file global attributes nc_dims : list A Python list of the NetCDF file dimensions nc_vars : list A Python list of the NetCDF file variables ''' def print_ncattr(key): """ Prints the NetCDF file attributes for a given key Parameters ---------- key : unicode a valid netCDF4.Dataset.variables key """ try: print "\t\ttype:", repr(nc_fid.variables[key].dtype) for ncattr in nc_fid.variables[key].ncattrs(): print '\t\t%s:' % ncattr,\ repr(nc_fid.variables[key].getncattr(ncattr)) except KeyError: print "\t\tWARNING: %s does not contain variable attributes" % key # NetCDF global attributes nc_attrs = nc_fid.ncattrs() if verb: print "NetCDF Global Attributes:" for nc_attr in nc_attrs: print '\t%s:' % nc_attr, repr(nc_fid.getncattr(nc_attr)) nc_dims = [dim for dim in nc_fid.dimensions] # list of nc dimensions # Dimension shape information. if verb: print "NetCDF dimension information:" for dim in nc_dims: print "\tName:", dim print "\t\tsize:", len(nc_fid.dimensions[dim]) print_ncattr(dim) # Variable information. nc_vars = [var for var in nc_fid.variables] # list of nc variables if verb: print "NetCDF variable information:" for var in nc_vars: if var not in nc_dims: print '\tName:', var print "\t\tdimensions:", nc_fid.variables[var].dimensions print "\t\tsize:", nc_fid.variables[var].size print_ncattr(var) return nc_attrs, nc_dims, nc_vars ############################# pdf = PdfPages('out.pdf') nc_m = './mask.nc' nc_mask = Dataset(nc_m, 'r') omask = nc_mask.variables["sftlf"][:] # shape is time, lat, lon as shown above mask=omask[0, :, :] nc_r = './data2.nc' nc_mod2 = Dataset(nc_r, 'r') nc_attrs, nc_dims, nc_vars = ncdump(nc_mod2) nc_f = './data1.nc' nc_mod1 = Dataset(nc_f, 'r') nc_attrs, nc_dims, nc_vars = ncdump(nc_mod1) # Extract data from NetCDF file if ('longitude' in nc_vars): lons = nc_mod1.variables['longitude'][:] else : if ('rlon' in nc_vars): lons = nc_mod1.variables['rlon'][:] else : lons = nc_mod1.variables['lon'][:] if ('latitude' in nc_vars): lats = nc_mod1.variables['latitude'][:] # extract/copy the data else : if ('rlat' in nc_vars): lats = nc_mod1.variables['rlat'][:] # extract/copy the data else : lats = nc_mod1.variables['lat'][:] # extract/copy the data if ('time' in nc_vars): time = nc_mod1.variables['time'][:] ntime='time' else : time = nc_mod1.variables['leadtime'][:] ntime='leadtime' if ('lev' in nc_vars): levs = nc_mod1.variables['lev'][:] # extract/copy the data else : if ('depth' in nc_vars): levs = nc_mod1.variables['depth'][:] # extract/copy the data else : levs = nc_mod1.variables['level'][:] # extract/copy the data var_name='VAR_NAME' data1 = nc_mod1.variables[var_name][:] # shape is time, lat, lon as shown above data2 = nc_mod2.variables[var_name][:] # shape is time, lat, lon as shown above # Plot of global temperature on our random day fig = plt.figure() fig.subplots_adjust(left=0.1, right=0.9, bottom=0.1, top=0.9) m = Basemap(projection='cyl', llcrnrlat=-90, urcrnrlat=90,\ llcrnrlon=0, urcrnrlon=360, resolution='c', lon_0=0) if (lons[1] < 0 ): m = Basemap(projection='cyl', llcrnrlat=-90, urcrnrlat=90,\ llcrnrlon=-280, urcrnrlon=80, resolution='c', lon_0=-279.5) lon2d, lat2d = np.meshgrid(lons, lats) x, y = m(lon2d, lat2d) time_idx = [TIME_INDEX] input_lev=levs[0] for i in range(len(time_idx)): if ('wo' not in nc_vars): d1=data1[time_idx[i],0, :, :] d2=data2[time_idx[i],0, :, :] mask1=np.where(mask==100.0,np.nan,1) d1=d1*mask1 d2=d2*mask1 else : d1=data1[time_idx[i],0, :, :]*100000 d2=data2[time_idx[i],0, :, :]*100000 cs = plt.subplot(211,axisbg='black') cs = m.contourf(x, y, d1, 21, cmap=plt.cm.Spectral_r, extend='both') levels = cs.levels m.drawparallels(np.arange(-90.,120.,30.), labels=[1,0,0,0], fontsize=7) m.drawmeridians(np.arange(0.,360.,60.), labels=[0,0,0,1], fontsize=7) if ( nc_mod1.variables[ntime].units == 'hours' ) : plt.title("%s on %s lead days :\n %s" % (nc_mod1.variables[var_name].long_name, time[time_idx[i]]/24.0, "MODEL TYPE1")) else: plt.title("%s on %s lead days :\n %s" % (nc_mod1.variables[var_name].long_name, time[time_idx[i]], "MODEL TYPE1")) cbar = plt.colorbar(cs, orientation='vertical', shrink=0.5) cbar.ax.tick_params(labelsize=6) if ('wo' not in nc_vars): cbar.set_label("%s at %s m (%s)" % (var_name, input_lev, nc_mod1.variables[var_name].units)) else : cbar.set_label("%s at %s m (1E-5 %s)" % (var_name, input_lev, nc_mod1.variables[var_name].units)) cf = plt.subplot(212,axisbg='black') cf = m.contourf(x, y, d2, levels, cmap=plt.cm.Spectral_r, extend='both') m.drawparallels(np.arange(-90.,120.,30.), labels=[1,0,0,0],fontsize=7) m.drawmeridians(np.arange(0.,360.,60.), labels=[0,0,0,1],fontsize=7) cbar = plt.colorbar(cf, orientation='vertical', shrink=0.5) cbar.ax.tick_params(labelsize=6) if ('wo' not in nc_vars): cbar.set_label("%s at %s m (%s)" % (var_name, input_lev, nc_mod2.variables[var_name].units)) else : cbar.set_label("%s at %s m (1E-5 %s)" % (var_name, input_lev, nc_mod2.variables[var_name].units)) plt.title("MODEL TYPE2") pdf.savefig(fig) plt.clf() # Close original NetCDF file. pdf.close() nc_mod1.close() nc_mod2.close()