''' NAME NetCDF with Python PURPOSE To demonstrate how to read and write data with NetCDF files using a NetCDF file from the NCEP/NCAR Reanalysis. Plotting using Matplotlib and Basemap is also shown. PROGRAMMER(S) Chris Slocum REVISION HISTORY 20140320 -- Initial version created and posted online 20140722 -- Added basic error handling to ncdump Thanks to K.-Michael Aye for highlighting the issue REFERENCES netcdf4-python -- http://code.google.com/p/netcdf4-python/ NCEP/NCAR Reanalysis -- Kalnay et al. 1996 http://dx.doi.org/10.1175/1520-0477(1996)077<0437:TNYRP>2.0.CO;2 ''' import os import errno import datetime as dt # Python standard library datetime module import numpy as np import netCDF4 from netCDF4 import Dataset # http://code.google.com/p/netcdf4-python/ import matplotlib.pyplot as plt import matplotlib as mpl from mpl_toolkits.basemap import Basemap, addcyclic, shiftgrid ''' from simplekml import (Kml, OverlayXY, ScreenXY, Units, RotationXY, AltitudeMode, Camera) def make_kml(llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat, figs, colorbar=None, **kw): """TODO: LatLon bbox, list of figs, optional colorbar figure, and several simplekml kw...""" kml = Kml() altitude = kw.pop('altitude', 2e7) roll = kw.pop('roll', 0) tilt = kw.pop('tilt', 0) altitudemode = kw.pop('altitudemode', AltitudeMode.relativetoground) camera = Camera(latitude=np.mean([urcrnrlat, llcrnrlat]), longitude=np.mean([urcrnrlon, llcrnrlon]), altitude=altitude, roll=roll, tilt=tilt, altitudemode=altitudemode) kml.document.camera = camera draworder = 0 for fig in figs: # NOTE: Overlays are limited to the same bbox. draworder += 1 ground = kml.newgroundoverlay(name='GroundOverlay') ground.draworder = draworder ground.visibility = kw.pop('visibility', 1) ground.name = kw.pop('name', 'overlay') ground.color = kw.pop('color', '9effffff') ground.atomauthor = kw.pop('author', 'ocefpaf') ground.latlonbox.rotation = kw.pop('rotation', 0) ground.description = kw.pop('description', 'Matplotlib figure') ground.gxaltitudemode = kw.pop('gxaltitudemode', 'clampToSeaFloor') ground.icon.href = fig ground.latlonbox.east = llcrnrlon ground.latlonbox.south = llcrnrlat ground.latlonbox.north = urcrnrlat ground.latlonbox.west = urcrnrlon if colorbar: # Options for colorbar are hard-coded (to avoid a big mess). screen = kml.newscreenoverlay(name='ScreenOverlay') screen.icon.href = colorbar screen.overlayxy = OverlayXY(x=0, y=0, xunits=Units.fraction, yunits=Units.fraction) screen.screenxy = ScreenXY(x=0.015, y=0.075, xunits=Units.fraction, yunits=Units.fraction) screen.rotationXY = RotationXY(x=0.5, y=0.5, xunits=Units.fraction, yunits=Units.fraction) screen.size.x = 0 screen.size.y = 0 screen.size.xunits = Units.fraction screen.size.yunits = Units.fraction screen.visibility = 1 kmzfile = kw.pop('kmzfile', 'overlay.kmz') #doc = kml.newdocument(kmzfile) kml.save("overlay.kml") kml.savekmz(kmzfile) def make_ground_overlay_kml(llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat, figs, colorbar=None, **kw): kml = Kml() ground = kml.newgroundoverlay(name='GroundOverlay') ground.icon.href = figs[0] ground.latlonbox.east = llcrnrlon ground.latlonbox.south = llcrnrlat ground.latlonbox.north = urcrnrlat ground.latlonbox.west = urcrnrlon kml.save("GroundOverlay.kml") def practice(request): map = Basemap(projection='merc' , resolution='i' , fix_aspect=True, llcrnrlon=119.0 , llcrnlat=21.8, urcrnlon=122.05 , urcrnrlat=25.4, lat_ts =20) ''' def ncdump(nc_fid, verb=True): ''' 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 def createFig(air, lons, lats, filename, depth_idx): path = filename+'/' fig = plt.figure(figsize=(20,20)) fig.subplots_adjust(left=None, bottom=None, right=None, top=None, \ wspace=None, hspace=None) axes = fig.add_subplot(1,1,1) ll_lat = 0 ll_lon = 100 ur_lat = 50 ur_lon = 150 m = Basemap(projection='merc', llcrnrlat=ll_lat,urcrnrlat=ur_lat,\ llcrnrlon=ll_lon,urcrnrlon=ur_lon, resolution='i') m.drawcoastlines() m.drawmapboundary(fill_color='aqua') air_cyclic, lons_cyclic = addcyclic(air[:, :], lons) # Create 2D lat/lon arrays for Basemap lon2d, lat2d = np.meshgrid(lons_cyclic, lats) # Transforms lat/lon into plotting coordinates for projection x, y = m(lon2d, lat2d) # Plot of air temperature with 11 contour intervals cs = m.contourf(x, y, air_cyclic, 21, cmap=plt.cm.Spectral_r) fig.tight_layout() fig.savefig(path+filename+'_'+depth_idx+'.png', transparent=True, format='png', pad_inches=0, bbox_inches='tight') if filename == "P": fig = plt.figure(figsize=(1.5, 4.0), facecolor=None, frameon=False) else : fig = plt.figure(figsize=(1.0, 4.0), facecolor=None, frameon=False) ax = fig.add_axes([0.05, 0.05, 0.2, 0.9]) cb = fig.colorbar(cs, cax=ax) #cb.set_label('Mean Dynamic Topography [m]', rotation=-90, color='k', labelpad=20) #if 'units' in nc_fid.variables[filename].ncattrs(): if filename == "currents": cb.set_label("%s (%s)" % ("", nc_fid.variables["U"].units), rotation=0, color='k', labelpad=15) elif 'units' in nc_fid.variables[filename].ncattrs() : cb.set_label("%s (%s)" % ("", nc_fid.variables[filename].units), rotation=0, color='k', labelpad=15) else : cb.set_label("%s (%s)" % ("", ""), rotation=0, color='k', labelpad=15) fig.savefig(path+filename+'_'+depth_idx+'_cbar.png', transparent=False, format='png') plt.close("all") def drawVector(U, V, lons, lats, filename, depth): path = filename+'/' ll_lat = 0 ll_lon = 100 ur_lat = 50 ur_lon = 150 fig = plt.figure(figsize=(20,20)) # KML friendly image. If using basemap try: `fix_aspect=False`. #ax = fig.add_axes([0, 0, 1, 1]) fig.subplots_adjust(left=None, bottom=None, right=None, top=None, \ wspace=None, hspace=None) ax = fig.add_subplot(1,1,1) ax.set_xlim(ll_lon, ur_lon) ax.set_ylim(ll_lat, ur_lat) ax.set_xbound(ll_lon, ur_lon) ax.set_ybound(ll_lat, ur_lat) u_cyclic, lons_cyclic = addcyclic(U[:, :], lons) #v_cyclic, lons_cyclic = addcyclic(V[:, :], lons) lon2d, lat2d = np.meshgrid(lons_cyclic, lats) #u2d, v2d = np.meshgrid(u_cyclic, v_cyclic) #u, v = m(u, v) m = Basemap(projection='merc', llcrnrlat=ll_lat,urcrnrlat=ur_lat,\ llcrnrlon=ll_lon,urcrnrlon=ur_lon, resolution='i') #m.drawcoastlines() x, y = m(lon2d, lat2d) m.drawmapboundary(fill_color='aqua') m.quiver(x[::10, ::10], y[::10, ::10], U[::10, ::10], V[::10, ::10], scale=240, scale_units='inches', units='inches') U_cyclic, lons_cyclic = addcyclic(U[:, :], lons) V_cyclic, lons_cyclic = addcyclic(V[:, :], lons) # Create 2D lat/lon arrays for Basemap lon2d, lat2d = np.meshgrid(lons_cyclic, lats) ''' #fig, ax = gearth_fig(llcrnrlon=ll_lon, llcrnrlat=ll_lat, urcrnrlon=ur_lon, urcrnrlat=ur_lat, pixels=pixels) #Q = plt.quiver(u, v) Q = ax.quiver(lon2d[::10, ::10], lat2d[::10, ::10], U_cyclic[::10, ::10], V_cyclic[::10, ::10], \ scale=90, units='x', angles='uv', scale_units='x') #Q = ax.quiver(lon[::10, ::10], lat[::10, ::10], u[::10, ::10], v[::10, ::10], scale=3) #ax.quiverkey(Q, 0.86, 0.45, 1, '1 m s$^{-1}$', labelpos='W') ax.set_axis_off() ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) ax.get_xaxis().get_label().set_visible(False) ax.get_yaxis().get_label().set_visible(False) ''' plt.axis('off') fig.tight_layout() #plt.tight_layout() fig.savefig(path+filename+'_'+depth+'_v.png', transparent=True, format='png', pad_inches=0, bbox_inches='tight') plt.close("all") def checkDir(): dir = ['P', 'S', 'T', 'currents', 'e_p', 'heat_up', 'swr_down', 'taux', 'tauy'] for d in dir : try: os.stat(d) except: os.mkdir(d) ''' # Plot of global temperature on our random day fig = plt.figure(figsize=(1, 16)) ax1 = fig.add_axes([0.05, 0.80, 0.9, 0.15]) norm = mpl.colors.Normalize(vmin=np.nanmin(a=air_cyclic), vmax=np.nanmax(a=air_cyclic)) cbar = mpl.colorbar.ColorbarBase(ax1, orientation='vertical', cmap=plt.cm.Spectral_r, norm=norm) cbar.set_label("%s (%s)" % ("var_desc", nc_fid.variables['U'].units)) #plt.title("%s on %s" % (nc_fid.variables['air'].var_desc, cur_time)) fig.savefig(filename+'_cbar.png', transparent=True, format='png', pad_inches=0, bbox_inches='tight') ''' checkDir() print np.__version__ print netCDF4.__version__ print mpl.__version__ nc_f = './member001_T1_00644.nc' # Your filename nc_fid = Dataset(nc_f, 'r') # Dataset is the class behavior to open the file # and create an instance of the ncCDF4 class nc_attrs, nc_dims, nc_vars = ncdump(nc_fid) # Extract data from NetCDF file lats32 = nc_fid.variables['Latitude'][:] # extract/copy the data lons32 = nc_fid.variables['Longitude'][:] lats_ = nc_fid.variables['Latitude'][:][1:487] lons_ = nc_fid.variables['Longitude'][:][1:405] depth = nc_fid.variables['Depth'][:] time32 = nc_fid.variables['Time'][:] air = nc_fid.variables['U'][0][:] # shape is time, lat, lon as shown above V = nc_fid.variables['V'][0][:] # shape is time, lat, lon as shown above P = nc_fid.variables['P'][0][:] # shape is time, lat, lon as shown above T = nc_fid.variables['T'][0][:] # shape is time, lat, lon as shown above S = nc_fid.variables['S'][0][:] # shape is time, lat, lon as shown above time = [np.float64(t)\ for t in time32] time_idx = 0 # some random day in 2012 # Python and the renalaysis are slightly off in time so this fixes that problem offset = dt.timedelta(hours=48) # List of all times in the file as datetime objects dt_time = [dt.date(1, 1, 1) + dt.timedelta(hours=t) - offset\ for t in time] cur_time = dt_time[time_idx] ''' for i in range(0, 25, 1): P = nc_fid.variables['P'][0][i] # shape is time, lat, lon as shown above T = nc_fid.variables['T'][0][i] # shape is time, lat, lon as shown above S = nc_fid.variables['S'][0][i] # shape is time, lat, lon as shown above U = nc_fid.variables['U'][0][i] # shape is time, lat, lon as shown above V = nc_fid.variables['V'][0][i] # shape is time, lat, lon as shown above currents = np.sqrt(U*U+V*V) createFig(P[1:487, 1:405], lons_, lats_, 'P', str(i)) createFig(T[1:487, 1:405], lons_, lats_, 'T', str(i)) createFig(S[1:487, 1:405], lons_, lats_, 'S', str(i)) createFig(currents[1:487, 1:405], lons_, lats_, 'currents', str(i)) drawVector(U[1:487, 1:405], V[1:487, 1:405], lons_, lats_, 'currents', str(i)) #createFig(P, lons32, lats32, 'P', str(i)) #createFig(T, lons32, lats32, 'T', str(i)) #createFig(S, lons32, lats32, 'S', str(i)) #createFig(currents, lons32, lats32, 'currents', str(i)) #drawVector(U, V, lons32, lats32, 'currents', str(i)) nc_f = './member001_fluxes_00644.nc' # Your filename nc_fid = Dataset(nc_f, 'r') # Dataset is the class behavior to open the file # and create an instance of the ncCDF4 class nc_attrs, nc_dims, nc_vars = ncdump(nc_fid) heat_up = nc_fid.variables['heat_up'] swr_down = nc_fid.variables['swr_down'] taux = nc_fid.variables['taux'] tauy = nc_fid.variables['tauy'] e_p = nc_fid.variables['e_p'] createFig(heat_up, lons_, lats_, 'heat_up', '0') createFig(swr_down, lons_, lats_, 'swr_down', '0') createFig(taux, lons_, lats_, 'taux', '0') createFig(tauy, lons_, lats_, 'tauy', '0') createFig(e_p, lons_, lats_, 'e_p', '0') ''' #plt.show() # Close original NetCDF file. nc_fid.close()