# plot 2D 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


def ncdump(nc_mod, 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_mod : 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_mod.variables[key].dtype)
            for ncattr in nc_mod.variables[key].ncattrs():
                print '\t\t%s:' % ncattr,\
                      repr(nc_mod.variables[key].getncattr(ncattr))
        except KeyError:
            print "\t\tWARNING: %s does not contain variable attributes" % key

    # NetCDF global attributes
    nc_attrs = nc_mod.ncattrs()
    if verb:
        print "NetCDF Global Attributes:"
        for nc_attr in nc_attrs:
            print '\t%s:' % nc_attr, repr(nc_mod.getncattr(nc_attr))
    nc_dims = [dim for dim in nc_mod.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_mod.dimensions[dim])
            print_ncattr(dim)
    # Variable information.
    nc_vars = [var for var in nc_mod.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_mod.variables[var].dimensions
                print "\t\tsize:", nc_mod.variables[var].size
                print_ncattr(var)
    return nc_attrs, nc_dims, nc_vars

#############################
pdf = PdfPages('out.pdf')

nc_r = './ref_data.nc'  
nc_ref = Dataset(nc_r, 'r') 
nc_attrs, nc_dims, nc_vars = ncdump(nc_ref)
var_name='VAR_NAME'
ref = nc_ref.variables[var_name][:]  # shape is time, lat, lon as shown above

nc_f = './my_data.nc'  
nc_mod = Dataset(nc_f, 'r') 
nc_attrs, nc_dims, nc_vars = ncdump(nc_mod)

# Extract data from NetCDF file
if ('longitude' in nc_vars):
    lons = nc_mod.variables['longitude'][:]
else :
    lons = nc_mod.variables['lon'][:]

if ('latitude' in nc_vars):
    lats = nc_mod.variables['latitude'][:]  # extract/copy the data
else :
    lats = nc_mod.variables['lat'][:]  # extract/copy the data

if ('lev' in nc_vars):
     levs = nc_mod.variables['lev'][:]  # extract/copy the data
else :
  if ('level' in nc_vars):
     levs = nc_mod.variables['level'][:]  # extract/copy the data
  else :
     levs = nc_mod.variables['plev'][:]  # extract/copy the data

if ('time' in nc_vars):
   time = nc_mod.variables['time'][:]
else :
   time = nc_mod.variables['leadtime'][:]

var_name='VAR_NAME'
data = nc_mod.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)
lon2d, lat2d = np.meshgrid(lons, lats)
x, y = m(lon2d, lat2d)

time_idx = [TIME_INDEX]
input_lev=LEVELS

for i in range(len(time_idx)):
  for j in range(len(levs)):
   if (input_lev == levs[j]):
     data1 = np.subtract(data[time_idx[i],j, :, :], ref[time_idx[i],j, :, :])
     if ( 'ta' in nc_vars ):
        con_levels = np.arange(-10.0, 10, 1)
     else:
       if ( 'ua' in nc_vars or 'va' in nc_vars ):
          con_levels = np.arange(-5.0, 5, 0.5)
       else:
         con_levels = np.arange(-0.01, 0.01, 0.001)

     cs = m.contourf(x, y, data1, con_levels, cmap=plt.cm.Spectral_r, extend='both')
     m.drawcoastlines()
     m.drawmapboundary()
     m.drawparallels(np.arange(-90.,120.,30.), labels=[1,0,0,0])
     m.drawmeridians(np.arange(0.,360.,60.), labels=[0,0,0,1])
     cbar = plt.colorbar(cs, orientation='horizontal', shrink=0.6)
     cbar.set_label("%s at %s hPa (%s)" % (var_name, levs[j], nc_mod.variables[var_name].units))
     if ( nc_mod.variables['leadtime'].units == 'hours' ) :
       plt.title("%s on %s lead days : \n%s" % (nc_mod.variables[var_name].long_name, time[time_idx[i]]/24.0, "MODEL-REF"))
     else:
       plt.title("%s on %s lead days : \n%s" % (nc_mod.variables[var_name].long_name, time[time_idx[i]], "MODEL-REF"))

     plt.contour(cs,cs.levels,colors='k',linewidths=0.5)
     pdf.savefig(fig)
     plt.clf()

#fig.savefig('foo.png',dpi = 200)
#pdf.savefig(fig)

# Close original NetCDF file.
pdf.close()
nc_mod.close()