#This script is used to make misc. maps from NCEP data,
#illustrating synoptic eddy activity,relative humidity patterns,
#and variations in static stability.
#
#It uses the Nio package (distributed as PyNio from the same
#place as PyNgl) for reading the data files.


#ToDo:  
#       *Include the full time-dependent dataset so other
#        days can be looked at and movies can be made. Make
#        Quicktime animations.
#
#       *Perhaps the computation of the seasonal cycle map
#        in Chapter 7 should be merged into this code. 

Chapter = "9"
Figure = "**FIG"

from ClimateUtilities import *
import Nio

NCEPpath = '/Users/rtp1/bigdata/ClimateBook/' #Edit this to appropriate data directory
#Note that the NCEP netcdf files are not currently in the WorkbookDatasets.
#They should be moved there.


#This is a specialized plotting routine for contour plots
#on map projections, with present Earth continent outlines
#It is not very extensively used in the courseware, so it
#hasn't been put in ClimateUtilities.  If you want to
#make plots of various aspects of present Earth climate
#based on netCDF data, you can use this routine as a model
#This doesn't have the features built in to allow saving the map
#to a file. Needs a lot of prettying up.
def map(v,**kwargs):
    rw = Dummy()
    #Pick the color map. "gsdtol" is greyscale
    if field == 'T':
        rw.wkColorMap = "temp1"
    elif field == 'rh':
        rw.wkColorMap = "hotres"
    else:
        rw.wkColorMap = "temp1"
    w = Ngl.open_wks('x11','Climate Workbook',rw)
    #
    r = Dummy()
    r.cnFillOn = True
    r.cnLinesOn = False
    r.cnLineLabelsOn = False # Set whether we want line labels
    
    r.pmLabelBarDisplayMode = "Always"
    r.mpProjection      ="CylindricalEquidistant" #"Mollweide"
    r.mpGridAndLimbOn   = True
    r.mpGridLatSpacingF = 30.
    r.mpGridLonSpacingF = 30.
    r.mpFillOn = False

    #The following tweaks the contour levels
    #field is a global that lets us select which
    #levels we want
    if field == 'T':
        r.cnLevelSelectionMode = "ManualLevels"
        r.cnMinLevelValF       = 240
        r.cnMaxLevelValF       = 290
        r.cnLevelSpacingF      = 5
    elif field == 'rh':
        r.cnLevelSelectionMode = "ManualLevels"
        r.cnMinLevelValF       = 10.
        r.cnMaxLevelValF       = 100.
        r.cnLevelSpacingF      = 5
    else:
        pass
    
    #I think the following only apply for the Mollweide case
    #r.vpXF = .1
    #r.vpYF = .9
    #r.vpWidthF = .8
    #r.vpHeightF = .8

    try:
        #Look for axis data in the array (Just for cdms variables)
        r.sfXArray= v.getLongitude()[:]
        r.sfYArray = v.getLatitude()[:]
    except:
        #If it's not there, it has to be put in as an argument
        r.sfXArray  = kwargs['x']
        r.sfYArray  = kwargs['y']

    #Now make the plot
    r.nglDraw = False
    r.nglFrame = False
    plot = Ngl.contour_map(w,v,r)#cdms vbls need v._data
    #Now draw the plot
    r.nglDraw = True
    Ngl.panel(w,[plot],[1,1],r)
    return plotObj(w,plot,rw) #So user can delete or save plot

#Open the file and get the temperature data
f = Nio.open_file(NCEPpath + 'T700_2003.nc')

#f.variables.keys() will return a list of the variables in
#the file if you don't already know.  From previous experience,
#I know that the temperature in this file is called 'T700'

istart = 178
T = f.variables['T700'][istart:istart + 50]
# NCEP daily data is already in Kelvins

#Now make a map of a sequence of days
lat = f.variables['lat'][:]
lon = f.variables['lon'][:].copy()
#w = contour(diff,x=lon,y=lat) #Use simple contour plotter
lon[-1]=360. #Kludge to eliminate the bar at the greenwich meridian
field = 'T' #Selects the nice contour levels for temperature
for i in range(10):
    w = map(T[4*i],x=lon,y=lat) #Use fancy mapmaker
    w.save('map%d'%i)

#Now do a relative humidity map
field = "rh"
f = Nio.open_file(NCEPpath + 'rh600Jan1_2001.nc')
rh = f.variables['rh600Jan1']
w = map(rh[:,:],x=lon,y=lat)