#!/usr/bin/env python

#-------------
# Load modules
#-------------
from netCDF4 import Dataset
import matplotlib.pyplot as plt            # pyplot module import
import numpy as np                         # Numpy import
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap   # basemap import
import datetime
import time
import sys

#----------------------
# Provide the file name
#----------------------
dirName  = '/discover/nobackup/jkouatch/gmiMetFields/MERRA/2x2.5/'
fileName = dirName+'MERRA300.prod.assim.20080316.2x2.5x72.nc'

#------------------
# Opening the file
#------------------
ncFid = Dataset(fileName, mode='r')

#----------------------
# Extracting dimensions
#----------------------
time = ncFid.variables['time'][:]
lev  = ncFid.variables['lev'][:]
lat  = ncFid.variables['lat'][:]
lon  = ncFid.variables['lon'][:]

#---------------------
# Extracting variables
#---------------------
slp = 0.01*ncFid.variables['SLP'][:]

#-------------
# Closing file
#-------------
ncFid.close()

nlat = lat.size - 1
nlon = lon.size - 1

#-------------------------------------------
# Pull out the first time record of the data
#-------------------------------------------
mySLP = slp[0,:,:]

#---------------------------------
# Create figure and axes instances
#---------------------------------
#fig = plt.figure(1,figsize=(15,8),dpi=50) # display 600x800
fig = plt.figure(1,figsize=(15,8),dpi=75)  # display 1024x768
ax = fig.add_axes([0.05,0.05,0.9,0.85])

#--------------------------------------------
# Create polar stereographic Basemap instance.
#--------------------------------------------

m = Basemap(projection='mill', 
            llcrnrlat=lat[0],
            urcrnrlat=lat[-1],
            llcrnrlon=lon[0],
            urcrnrlon=lon[-1] )

#----------------------------------------
# Draw the coastlines of continental area
#----------------------------------------
m.drawcoastlines(linewidth=1.25)

#--------------------------
# Fill the continental area
#--------------------------
m.fillcontinents(color='0.8')

#-----------------------------
# Draw parallels and meridians
#-----------------------------
m.drawparallels(np.arange(-80,81,20),labels=[1,1,0,0])
m.drawmeridians(np.arange(-180,180,60),labels=[0,0,0,1])

#-----------------------------
# Display the data on the axes
#-----------------------------
im = m.imshow(mySLP,
             interpolation='nearest',
             extent=[lon[0], lon[-1], lat[0], lat[-1]],
             cmap=plt.cm.jet)

plt.colorbar(orientation='hoirzontal',shrink=.8)

plt.title('Sea Level Pressure')
plt.savefig('fig_slp.png')
plt.show()
#plt.clf()