def get_era_winds(m, rawdatapath, start_year, end_year, start_month,  end_month, xyres):

	wind_data = Dataset(rawdatapath+'/WINDS/ERA/DATA/ERAI_WINDS_MONTHLY_1979-2014.nc', 'r')
	lats = wind_data.variables['latitude'][::-1]
	lons = wind_data.variables['longitude'][:]
	time = wind_data.variables['time'][:]/(24.*30.)
	time = time-time[0]
	time=time.reshape(time.shape[0]/12,12)
	u10 = wind_data.variables['u10'][:, ::-1, :]
	v10 = wind_data.variables['v10'][:, ::-1, :]
	u10=u10.reshape(u10.shape[0]/12, 12,u10.shape[1],u10.shape[2])
	v10=v10.reshape(v10.shape[0]/12, 12, v10.shape[1],v10.shape[2])

	u10_winter_mean= np.mean(u10[start_year-1979:end_year-1979+1,start_month:end_month+1], axis=tuple(range(0, 2)))
	v10_winter_mean= np.mean(v10[start_year-1979:end_year-1979+1,start_month:end_month+1], axis=tuple(range(0, 2)))


	u10_winter_meanS, lonsS = shiftgrid(180.,u10_winter_mean,lons,start=False)
	v10_winter_meanS, lonsS = shiftgrid(180.,v10_winter_mean,lons,start=False)

	u10_winter_meanSC, lonsSC = addcyclic(u10_winter_meanS, lonsS)
	v10_winter_meanSC, lonsSC = addcyclic(v10_winter_meanS, lonsS)

	xyres=100

	xvel,yvel,xptsW,yptsW = m.transform_vector(u10_winter_meanSC,v10_winter_meanSC,lonsSC,lats,xyres,xyres,returnxy=True,masked=True)
	wind_speed = sqrt((xvel**2) + (yvel**2))


	wind_speed = sqrt((xvel**2) + (yvel**2))
	return xptsW, yptsW, xvel, yvel, wind_speed

def region_plot(fname, fname2, time, output_dir, julday):

	ncfile = nc.Dataset(fname,'r')

	fig = plt.figure(figsize=(20.48,10.24))
	ax = plt.axes([0,0,1,1],frameon=False)
	ax.set_axis_off()

	# set up a Basemap
	latrng = (-90.0,90.0)
	lonrng = (0,360)

	m = Basemap(projection='cyl',llcrnrlat=latrng[0],urcrnrlat=latrng[1],
        llcrnrlon=lonrng[0],urcrnrlon=lonrng[1],resolution='h',
        area_thresh=800,fix_aspect=False)

	lats = ncfile.variables["yu_ocean"][:]
	lons = ncfile.variables["xu_ocean"][:]
	u = ncfile.variables["u"][time,0,:,:]

	ncfile2 = nc.Dataset(fname2,'r')
	v = ncfile2.variables["v"][time,0,:,:]

	lons2 = lons

	u = ma.array(u)
	u = ma.masked_where(u <= -10.0, u)
	u = ma.masked_where(u >= 10.0, u)
	u = ma.masked_where(u == 0.0, u)

	u, lons = addcyclic(u, lons)

	v = ma.array(v)
	v = ma.masked_where(v <= -10.0, v)
	v = ma.masked_where(v >= 10.0, v)
	v = ma.masked_where(v == 0.0, v)

	v, lons2 = addcyclic(v, lons2)

	field = ((u**2 + v**2)**0.5)
	field = np.log(field)

	x,y = m(*np.meshgrid(lons[:],lats[:]))

	# cm.Blues looks pretty good
	colorMap = mpl.cm.Blues
	m.drawmapboundary(fill_color='white')
	nlevs = 256 
	clevs = np.linspace(-4.0,0.5,num=nlevs)
	
	cs = m.contourf(x,y,field, levels=clevs, cmap=colorMap, extend='both', norm=mpl.colors.normalize(),antialiased=False )
	m.fillcontinents(color='black', lake_color='black', zorder=1)
	save_file = "%s/%f.png" % (output_dir,julday)
	plt.savefig(save_file)
	
	ncfile.close()

    def plot(self):

        self.data, lonwrap = addcyclic(self.data, self.lons)

        # Sort latitudes and data
        lat_idx = np.argsort(self.lats)
        self.lats = self.lats[lat_idx]
        self.data = self.data[lat_idx]

        data_lon_min = min(lonwrap)
        data_lon_max = max(lonwrap)
        data_lat_min = min(self.lats)
        data_lat_max = max(self.lats)

        new_lons = np.arange(data_lon_min - 1.0, data_lon_max + 1.0, 1.0)
        new_lats = np.arange(data_lat_min - 1.0, data_lat_max + 1.0, 1.0)

        x, y = self.m(*np.meshgrid(new_lons[:], new_lats[:]))

        # Two pass interpolation to deal with the mask.
        # First pass does bilinear, the next does nearest neighbour
        # interpolation.
        # It's not clear this is working, and the problem is likely
        # solved by ensuring the right mask is used!
        data_bl = interp(self.data, lonwrap[:], self.lats[:], x, y, order=1)
        data_nn = interp(self.data, lonwrap[:], self.lats[:], x, y, order=0)

        data_bl[data_nn.mask == 1] = data_nn[data_nn.mask == 1]

        if self.parameters.has_key("color_levels"):
            self.__print_custom_color_plot(x, y, data_bl)
        else:
            self.__print_cmap_plot(x, y, data_bl)

        return self.main_render

def doplot(nc_file=None,varname=None,vmin=None,vmax=None,
           title=None):
    lons=nc_file.variables['lon'][...]
    lats=nc_file.variables['lat'][...]
    vals=nc_file.variables[varname]
    vals.set_auto_maskandscale(True)
    vals=vals[...]
    vals,lons=shiftgrid(180.,vals,lons,start=False)
    vals,lons=addcyclic(vals,lons)
    fig,ax = plt.subplots(1,1)
    ax.cla()
    cmap=cm.RdBu_r
    cmap.set_over('y')
    cmap.set_under('k')
    the_norm=Normalize(vmin=vmin,vmax=vmax,clip=False)
    params=dict(projection='moll',lon_0= 0,resolution='c')
    m = Basemap(**params)
    x, y = m(*np.meshgrid(lons, lats))
    im=m.pcolormesh(x,y,vals,cmap=cmap,norm=the_norm,ax=ax)
    cb=m.colorbar(im,extend='both',location='bottom')
    m.drawcoastlines()
    m.drawparallels(np.arange(-90.,120.,30.))
    m.drawmeridians(np.arange(0.,420.,60.))
    ax.set_title(title)
    return fig,m,ax,vals,lons

def test_run_hgt():
  #long term mean
  fname_mean = '/data01/forONR/hgt.mon.1981-2010.ltm.nc'
  lev_500 = 5 #level of 500hPa
  month = 8 #august
  
  dataMean = netCDF4.Dataset(fname_mean,'r')
  hgtMean = dataMean.variables['hgt'][month,lev_500,:,:] #lat,lon
  #hgtMean += dataMean.variables['hgt'][month-1,lev_500,:,:]
  #hgtMean += dataMean.variables['hgt'][month-2,lev_500,:,:]
  #hgtMean /= 3.
  latMean = dataMean.variables['lat'][:]*np.pi/180. #stored in degrees!
  lonMean = dataMean.variables['lon'][:]*np.pi/180.
  dataMean.close()
  
  nlat = len(latMean)
  nlon = len(lonMean)
  var = np.zeros((nlat,nlon), dtype=float)
  
  #mpas files
  #fnames = []
  fnames = sorted(glob.glob('/arctic1/nick/cases/vduda/x4/x4.t.output.2006-08-*'))
  fnames.extend(sorted(glob.glob('/arctic1/nick/cases/vduda/x4.t.output.2006-08-15*')))
  #fnames = sorted(glob.glob('/arctic1/nick/cases/vduda/x4/*.output.2006-08*'))
  counter = 0
  for iFile, fpath in enumerate(fnames):
    data = output_data.open_netcdf_data(fpath)
    nEdgesOnCell = data.variables['nEdgesOnCell'][:];
    cellsOnCell = data.variables['cellsOnCell'][:]-1;
    latCell = data.variables['latCell'][:];
    lonCell = data.variables['lonCell'][:];
    nTimes = len(data.dimensions['Time'])
    #nTimes = 1 #3
    
    nCells = len(latCell)
    geop_mpas = np.zeros(nCells,dtype=float)
    
    for iTime in xrange(nTimes):
      #get average value for mesh over these times so can even anomaly different meshes together
      output_data.print_xtime(data, iTime)
      hgt_mpas = data.variables['height_500hPa'][iTime,:]
      #average here to avoid summing to large number over many times
      geop_mpas += mpas_hgt2geopotential(nCells, hgt_mpas, latCell)/nTimes
      
    llMap = makeMap_ll2mpas(latMean, lonMean, latCell, lonCell, nEdgesOnCell, cellsOnCell)
    var += calc_diff_field(hgtMean, geop_mpas, llMap, nlat, nlon)
    #print var
    counter = counter+1
    
    data.close()
    
  var /= counter
  #var = hgtMean
  # add wrap-around point in longitude.
  var, lonMean = addcyclic(var, lonMean)
  np.savez('tmp.dat',var,latMean,lonMean)
  lats,lons = np.meshgrid(latMean, lonMean)
  lats *= 180./np.pi; lons *= 180./np.pi
  plot_anomaly_ll(lats, lons, np.transpose(var))

def region_plot(fname, time, output_dir, julday):

    #print "fname = %s, day = %d\n" % (fname, time);

	ncfile = nc.Dataset(fname,'r')

	fig = plt.figure(figsize=(20.48,10.24))
	ax = plt.axes([0,0,1,1],frameon=False)
	ax.set_axis_off()

	# set up a Basemap
	latrng = (-90.0,90.0)
	lonrng = (0,360)

	m = Basemap(projection='cyl',llcrnrlat=latrng[0],urcrnrlat=latrng[1],
        llcrnrlon=lonrng[0],urcrnrlon=lonrng[1],resolution='h',
        area_thresh=800,fix_aspect=False)

	# for mom output
	lats = ncfile.variables["yt_ocean"][:]
	lons = ncfile.variables["xt_ocean"][:]
	field = ncfile.variables["eta_t"][time,:,:]

	aveNcfile = nc.Dataset('ocean_eta_annual.nc','r')
	average = aveNcfile.variables["eta_t"][0,:,:]

	field = field - average

	field = ma.array(field)
	field = ma.masked_where(field == -1.0e10, field)
	field = ma.masked_where(field<= -10,field)
	field = ma.masked_where(field>= 100, field)

	field, lons = addcyclic(field, lons)

	#field_min = np.min(field) 
	#field_max = np.max(field)

	x,y = m(*np.meshgrid(lons[:],lats[:]))

	colorMap = colormaps.make_colormap({0.:'#191970',  0.2:'#448BB9',  0.35:'#1e90ff', 0.5:'w', 0.6:'#F8BC4E', 0.7:'#ffa500', 0.8:'#ff8c00',0.9:'#F56D2E',  1.0:'#CC4000'}) 

	m.drawmapboundary(fill_color='white')
	nlevs = 256
	clevs = np.linspace(-0.8,0.8,num=nlevs)
	
	cs = m.contourf(x,y,field, levels=clevs, cmap=colorMap, extend='both', norm=mpl.colors.normalize(),antialiased=False )
	m.fillcontinents(color='black', lake_color='black', zorder=1)
	save_file = "%s/%f.png" % (output_dir,julday)
	plt.savefig(save_file)
	
	#plt.show()
	ncfile.close()

def project_SkyMap_Mollweide( skymap , Pnorm=1 , Qnorm=1
                              , Ppath=None , Qpath=None ):

    Pw , lonsw = addcyclic( skymap.P , skymap.sky.lons )
    Qw , lonsw = addcyclic( skymap.Q , skymap.sky.lons )
    meshlon , meshlat = pylab.meshgrid( lonsw , skymap.sky.lats )

    projection = Basemap( projection='moll' , lon_0=180 , resolution='c' )
#    projection = Basemap( projection='ortho' , lat_0=60 , lon_0=180 , resolution='c' )
    projection.drawmapboundary()
    x , y = projection( meshlon , meshlat )

    if Ppath == None :
        pass
    else :
        fig = pyplot.figure( figsize = (8,8) )
        ax = fig.add_axes( [ 0.05 , 0.15 , 0.8 , 0.8 ] )
        projection.contourf( x , y , Pnorm*Pw , 30 )
        projection.drawmeridians( numpy.arange(0,360,30) )
        projection.drawparallels( numpy.arange(-90,90,30) , labels=[1,0,0,0] )
        pos = ax.get_position()
        l , b , w , h = pos.bounds
        cax = pyplot.axes( [ l+w+0.03 , b , 0.04 , h ] )
        pyplot.colorbar( cax=cax , orientation='vertical' , format='%.1e' )
        pylab.savefig( Ppath )

    if Qpath == None :
        pass
    else :
        fig = pyplot.figure( figsize = (8,8) )
        ax = fig.add_axes( [ 0.05 , 0.15 , 0.8 , 0.8 ] )
        projection.contourf( x , y , Qnorm*Qw , 30 )
        projection.drawmeridians( numpy.arange(0,360,30) )
        projection.drawparallels( numpy.arange(-90,90,30) , labels=[1,0,0,0] )
        pos = ax.get_position()
        l , b , w , h = pos.bounds
        cax = pyplot.axes( [ l+w+0.03 , b , 0.04 , h ] )
        pyplot.colorbar( cax=cax , orientation='vertical' )
        pylab.savefig( Qpath )
    return

 def plotMap(self,figdir=''):
     if not self.gotpix:
         self.getPixels()
     lats = self.sky.lats
     lons = self.sky.lons
     pixpw , lonsw = addcyclic( self.pixp , lons )
     pixqw , lonsw = addcyclic( self.pixq , lons )
     meshlon, meshlat = pylab.meshgrid( lonsw , lats )
     projection = Basemap(projection='moll',lon_0=180,resolution='c')
     x , y = projection( meshlon , meshlat )
     projection.contourf( x,y, pixpw , 30)
     projection.drawmeridians(numpy.arange(0,360,30))
     projection.drawparallels(numpy.arange(-90,90,30))
     projection.drawmapboundary()
     figname = self.outfileprefix + '-' + self.maptype + 'real.png'
     pylab.savefig( ( figdir + figname ) )
     projection.contourf( x,y, pixqw , 30)
     projection.drawmeridians(numpy.arange(0,360,30))
     projection.drawparallels(numpy.arange(-90,90,30))
     projection.drawmapboundary()
     figname = self.outfileprefix + '-' + self.maptype + 'imag.png'
     pylab.savefig( ( figdir + figname ) )

def shade_coord(xin,yin,datain=None,lon0=None):
	from mpl_toolkits.basemap import Basemap, shiftgrid, addcyclic
	from mcmath import edges

	# invert lat coords if they are descending instead of ascending
	if yin[-1] < yin[0]:
		yin = yin[::-1]
		if datain is not None:
			ydimloc = np.where(yin.size == np.array(datain.shape))[0]
			if len(ydimloc) == 0:
				raise MCPlotError("no dimension in 'datain' matches length of 'yin'")
			y_ind = []
			for i in xrange(datain.ndim):
				y_ind.append(slice(None))
			y_ind[ydimloc[0]] = slice(None,None,-1)
			datain = datain[y_ind]
	yedg = edges(yin)
	
	# convert xin to -180:180 range; roll and addcyclic as needed to have lon0 as central lon
	xedg = edges(xin)
	xspan = xedg[-1] - xedg[0]
	if ((xspan < 365.) & (xspan > 355.)):
		xin = np.where(xin >= 180,xin-360,xin)
		roll_ind = np.where(xin < np.roll(xin,1))[0]
		if len(roll_ind) == 0:
			raise MCPlotError("can't find pivot between +180 and -180 in xin")
		xin = np.roll(xin,-roll_ind)
		if datain is not None:
			xdimloc = np.where(xin.size == np.array(datain.shape))[-1]
			if len(xdimloc) == 0:
				raise MCPlotError("no dimension in 'datain' matches length of 'xin'")
			datain = np.roll(datain,-roll_ind,axis=xdimloc)

			if (lon0 is not None):
				(datain,xin) = addcyclic(datain,xin)
				(datain,xin) = shiftgrid(lon0-180,datain,xin)

			xedg = edges(xin)

			## Fix for keeping the squares near the edges of global maps filled just to the edge (map boundary)
			if (lon0 is not None):
				if xedg[0] < (lon0-180):
					xedg[0] = lon0-180
				if xedg[-1] > (lon0+180):
					xedg[-1] = lon0+180

			return [xedg,yedg,datain,xin,yin]
		else:
			xedg = edges(xin)

	return [xedg,yedg,datain,xin,yin]

def region_plot(fname, time, output_dir, julday):

	ncfile = nc.Dataset(fname,'r')

	fig = plt.figure(figsize=(20,10))
	ax = plt.axes([0,0,1,1],frameon=False)
	ax.set_axis_off()

	# set up a Basemap
	latrng = (-90.0,90.0)
	lonrng = (0,360)

	m = Basemap(projection='cyl',llcrnrlat=latrng[0],urcrnrlat=latrng[1],
        llcrnrlon=lonrng[0],urcrnrlon=lonrng[1],resolution='h',
        area_thresh=800,fix_aspect=False)

	# for ofam salinity 
	lats = ncfile.variables["yt_ocean"][:]
	lons = ncfile.variables["xt_ocean"][:]
	field = ncfile.variables["salt"][time,0,:,:]

	aveNcfile = nc.Dataset('ocean_salt_annual.nc','r')
	average = aveNcfile.variables["salt"][0,0,:,:]
	
	field = field - average

	field = ma.array(field)
	field = ma.masked_where(field == -1.0e10, field)
	field = ma.masked_where(field<= -10,field)
	field = ma.masked_where(field>= 100, field)

	field, lons = addcyclic(field, lons)

	x,y = m(*np.meshgrid(lons[:],lats[:]))

	colorMap =  colormaps.make_colormap({0.:'#010233', 0.3:'#01034B',0.4:'#02035F',0.42:'#28297B', 0.45:'#3B3D8E', 0.5:'w', 0.55:'#B32222',0.58:'#731515',0.6:'#600202',0.7:'#4B0202', 1.0:'#330101'})  # red to blue

	m.drawmapboundary(fill_color='white')
	nlevs = 256 
	clevs = np.linspace(-2,2,num=nlevs)
	
	cs = m.contourf(x,y,field, levels=clevs, cmap=colorMap, extend='both', norm=mpl.colors.normalize(),antialiased=False )
	m.fillcontinents(color='black', lake_color='black', zorder=1)
	save_file = "%s/%f.png" % (output_dir,julday)

	plt.savefig(save_file)
	ncfile.close()

def test_plot():
  fname_mean = '/data01/forONR/hgt.mon.1981-2010.ltm.nc'  
  dataMean = netCDF4.Dataset(fname_mean,'r')
  latMean = dataMean.variables['lat'][:]
  lonMean = dataMean.variables['lon'][:]
  dataMean.close()
  
  nlat = len(latMean)
  nlon = len(lonMean)
  var = np.zeros((nlat,nlon), dtype=float)
  for iLat in xrange(nlat):
    for iLon in xrange(nlon):
      var[iLat, iLon] = lonMean[iLon] #latMean[iLat]
  
  var, lonMean = addcyclic(var, lonMean)
  lats,lons = np.meshgrid(latMean, lonMean)
  plot_anomaly_ll(lats, lons, np.transpose(var))

def region_plot(fname, time, output_dir, julday):

	ncfile = nc.Dataset(fname,'r')

	fig = plt.figure(figsize=(20.48,10.24))
	ax = plt.axes([0,0,1,1],frameon=False)
	ax.set_axis_off()

	# set up a Basemap
	latrng = (-90.0,90.0)
	lonrng = (0,360)

	m = Basemap(projection='cyl',llcrnrlat=latrng[0],urcrnrlat=latrng[1],
        llcrnrlon=lonrng[0],urcrnrlon=lonrng[1],resolution='h',
        area_thresh=800,fix_aspect=False)

	# for salinity
	lats = ncfile.variables["yt_ocean"][:]
	lons = ncfile.variables["xt_ocean"][:]
	field = ncfile.variables["salt"][time,0,:,:]

	field = ma.array(field)
	field = ma.masked_where(field == -1.0e10, field)
	field = ma.masked_where(field<= -10,field)
	field = ma.masked_where(field>= 100, field)

	field, lons = addcyclic(field, lons)

	#field_min = np.min(field) 
	#field_max = np.max(field)
	
	x,y = m(*np.meshgrid(lons[:],lats[:]))

	# do one with jet colormap
	colorMap = mpl.cm.jet
	
	m.drawmapboundary(fill_color='white')
	nlevs = 256
	clevs = np.linspace(28,38,num=nlevs)
	
	cs = m.contourf(x,y,field, levels=clevs, cmap=colorMap, extend='both', norm=mpl.colors.normalize(),antialiased=False )
	m.fillcontinents(color='black', lake_color='black', zorder=1)
	save_file = "%s/%f.png" % (output_dir,julday)
	plt.savefig(save_file)
	
	ncfile.close()

def expand_longitude(
    test_lon,orig_lon,orig_data):
    """
    Makes a cyclic array using basemap and expands the longitudes to give the buffer aroud the edge needed 
    for the area weighted re-gridding.
    """
    import numpy as np
    from mpl_toolkits.basemap import addcyclic, shiftgrid
    import copy
    # shift the grid of the emissions so it fits with the test longitude grid
    # first, find the value of the orig long grid nearest to the test long grid
    if len(orig_data.shape) != 2:
        raise Exception, 'Data array must be two-dimensional'
    idx = (np.abs(orig_lon-test_lon[0])).argmin()
    start = orig_lon[idx]
    orig_data_temp,orig_lon_temp=addcyclic(orig_data,orig_lon)
    orig_data_temp,orig_lon_temp=shiftgrid(start,orig_data_temp,orig_lon_temp)
    test_dlon = test_lon[1] - test_lon[0]
    orig_dlon = orig_lon[1] - orig_lon[0]
    extra_cells=np.int(test_dlon/orig_dlon)+2
    # NOTE: DATA AND LOGITUDE HAVE BEEN MADE CYCLIC, ORIGIN IS THAT CLOSEST IN ORIG TO THE START OF TEST
    new_lon   = np.linspace(orig_lon_temp[0] - extra_cells*orig_dlon,\
                              orig_lon_temp[-2]+(extra_cells)*orig_dlon,\
                              len(orig_lon)+2*extra_cells)
    lon_index = np.arange(-1*extra_cells,len(new_lon)-extra_cells,1,dtype=np.int)
    new_data = np.zeros((orig_data.shape[0],len(new_lon)))
    i = 0
    while i < len(new_lon):
        # ignore cyclic grid here for length, since orig_data_temp[0] == orig_data_temp[1], but
        # we don't want to replicate that
        new_data[:,i] = copy.deepcopy(orig_data_temp[:,lon_index[i]%orig_data.shape[1]])
        i+=1
    del idx
    del i
    del start
    del orig_data_temp
    del orig_lon_temp
    del test_dlon
    del orig_dlon
    del test_lon
    del orig_lon
    del orig_data
    del extra_cells
    del lon_index
    return new_lon,new_data

 def __transform_lons(self,bbox,lons,var):
     """ Take Bounding box longitudes and transform them for sensible
     Basemap plotting.
     
     bbox: Bounding Box ... values will change in there
     lons: numpy array of longitudes
     var: numpy array of field to be plotted
     """
     
     to_360 = lambda x: ( x % 360. )
     
     # Put everything into a range greater 0 right away
     lon_min = to_360(bbox.lon_min)
     lon_max = to_360(bbox.lon_max)
     lons = to_360(lons)
     
     # If lon_min is greater than lon_max after transformation
     if lon_min > lon_max:
         lon_max += 360.
     
     # If both are the same i.e. 360 degree print
     if lon_min == lon_max:
         lon_max += 360.
     
     # If lon_max greater than 360, remap longitudes to that range
     if lon_max > 360.:
         idx = lons < lon_min
         lons[idx] = lons[idx] + 360.
     
     bbox.lon_min = lon_min
     bbox.lon_max = lon_max
     
     # Ensure longitudes are ascending, and shuffle field with same indices
     idx = np.argsort(lons)
     lons = lons[idx]
     var = var[:,idx]
     var, lons = addcyclic(var,lons)
     
     return bbox,lons,var

    def _meshgrid(self, data_array):
        """
        given some data, returns the np.meshgrid associated
        lons and lats + addcyclic longitude if proj is 'spstere'
        or 'npstere'
        """

        latitudes = self.dset_domain['latitudes'].data
        longitudes = self.dset_domain['longitudes'].data

        if self.proj in ['npstere','spstere','moll']:
            """
            if north polar or south polar stereographic projection
            we take care of the discontinuity at Greenwich using the
            `addcyclic` function of basemap
            """
            data_array, lon = addcyclic(data_array, longitudes)
            lons, lats = np.meshgrid(lon, latitudes)
        else:
            lons, lats = np.meshgrid(longitudes, latitudes)

        return lons, lats, data_array

lats, uwnd, vwnd = order_latdim(lats, uwnd, vwnd)

# Create a VectorWind instance to handle the computation of streamfunction and
# velocity potential.
w = VectorWind(uwnd, vwnd)

# Compute the streamfunction and velocity potential. Also use the bundled
# tools to re-shape the outputs to the 4D shape of the wind components as they
# were read off files.
sf, vp = w.sfvp()
sf = recover_data(sf, uwnd_info)
vp = recover_data(vp, uwnd_info)

# Pick out the field for December and add a cyclic point (the cyclic point is
# for plotting purposes).
sf_dec, lons_c = addcyclic(sf[11], lons)
vp_dec, lons_c = addcyclic(vp[11], lons)

# Plot streamfunction.
m = Basemap(projection='cyl', resolution='c', llcrnrlon=0, llcrnrlat=-90,
            urcrnrlon=360.01, urcrnrlat=90)
x, y = m(*np.meshgrid(lons_c, lats))
clevs = [-120, -100, -80, -60, -40, -20, 0, 20, 40, 60, 80, 100, 120]
m.contourf(x, y, sf_dec*1e-06, clevs, cmap=plt.cm.RdBu_r,
           extend='both')
m.drawcoastlines()
m.drawparallels((-90, -60, -30, 0, 30, 60, 90), labels=[1,0,0,0])
m.drawmeridians((0, 60, 120, 180, 240, 300, 360), labels=[0,0,0,1])
plt.colorbar(orientation='horizontal')
plt.title('Streamfunction ($10^6$m$^2$s$^{-1}$)', fontsize=16)

grid = AxesGrid(fig, [0.05,0.01,0.9,0.9],
                nrows_ncols=(3, 2),
		axes_pad=0.25,
		cbar_mode='single',
		cbar_pad=0.3,
		cbar_size=0.1,
                cbar_location='top',
                share_all=True,
		)

# create Basemap instance for Orthographic projection.
m = Basemap(lon_0=-90,lat_0=60,projection='ortho')
# add wrap-around point in longitude.
t2m = np.zeros((ntimes,nlats,nlons+1),np.float32)
for nt in range(ntimes):
    t2m[nt,:,:], lons = addcyclic(t2mvar[nt,:,:], lons1)
# convert to celsius.
t2m = t2m-273.15
# contour levels
clevs = np.arange(-30,30.1,2.)
lons, lats = np.meshgrid(lons, lats)
x, y = m(lons, lats)
# make subplots.
for nt,fcsthr in enumerate(fcsthrs):
    ax = grid[nt]
    m.ax = ax
    cs = m.contourf(x,y,t2m[nt,:,:],clevs,cmap=plt.cm.jet,extend='both')
    m.drawcoastlines(linewidth=0.5)
    m.drawcountries()
    m.drawparallels(np.arange(-80,81,20))
    m.drawmeridians(np.arange(0,360,20))

print verifdates
lats = latitudes[:]
nlats = len(lats)
lons1 = longitudes[:]
nlons = len(lons1)

# unpack 2-meter temp forecast data.

t2mvar = data.variables['tmp2m']
t2min = t2mvar[0:ntimes,:,:]
t2m = np.zeros((ntimes,nlats,nlons+1),t2min.dtype)
# create Basemap instance for Orthographic projection.
m = Basemap(lon_0=-90,lat_0=60,projection='ortho')
# add wrap-around point in longitude.
for nt in range(ntimes):
    t2m[nt,:,:], lons = addcyclic(t2min[nt,:,:], lons1)
# convert to celsius.
t2m = t2m-273.15
# contour levels
clevs = np.arange(-30,30.1,2.)
lons, lats = np.meshgrid(lons, lats)
x, y = m(lons, lats)
# create figure.
fig=plt.figure(figsize=(6,8))
# make subplots.
for nt,fcsthr in enumerate(fcsthrs):
    ax = fig.add_subplot(321+nt)
    cs = m.contourf(x,y,t2m[nt,:,:],clevs,cmap=plt.cm.jet,extend='both')
    m.drawcoastlines(linewidth=0.5)
    m.drawcountries()
    m.drawparallels(np.arange(-80,81,20))

# Compute components of rossby wave source: absolute vorticity, divergence,
# irrotational (divergent) wind components, gradients of absolute vorticity.
eta = w.absolutevorticity()
div = w.divergence()
uchi, vchi = w.irrotationalcomponent()
etax, etay = w.gradient(eta)

# Combine the components to form the Rossby wave source term. Re-shape the
# Rossby wave source array to the 4D shape of the wind components as they were
# read off files.
S = -eta * div - uchi * etax + vchi * etay
S = recover_data(S, uwnd_info)

# Pick out the field for December at 200 hPa and add a cyclic point (the
# cyclic point is for plotting purposes).
S_200, lons_c = addcyclic(S[11, 9], lons)

# Plot Rossby wave source.
m = Basemap(projection='cyl', resolution='c', llcrnrlon=0, llcrnrlat=-90,
        urcrnrlon=360.01, urcrnrlat=90)
x, y = m(*np.meshgrid(lons_c, lats))
clevs = [-30, -25, -20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30]
m.contourf(x, y, S_200*1e11, clevs, cmap=plt.cm.RdBu_r,
        extend='both')
m.drawcoastlines()
m.drawparallels((-90, -60, -30, 0, 30, 60, 90), labels=[1,0,0,0])
m.drawmeridians((0, 60, 120, 180, 240, 300, 360), labels=[0,0,0,1])
plt.colorbar(orientation='horizontal')
plt.title('Rossby Wave Source ($10^{-11}$s$^{-1}$)', fontsize=16)
plt.show()

lats, uwnd, vwnd = order_latdim(lats, uwnd, vwnd)

# Create a VectorWind instance to handle the computation of streamfunction and
# velocity potential.
w = VectorWind(uwnd, vwnd)

# Compute the streamfunction and velocity potential. Also use the bundled
# tools to re-shape the outputs to the 4D shape of the wind components as they
# were read off files.
sf, vp = w.sfvp()
sf = recover_data(sf, uwnd_info)
vp = recover_data(vp, uwnd_info)

# Pick out the field for December at 200 hPa and add a cyclic point (the
# cyclic point is for plotting purposes).
sf_200, lons_c = addcyclic(sf[11, 9], lons)
vp_200, lons_c = addcyclic(vp[11, 9], lons)

# Plot streamfunction.
m = Basemap(projection='cyl', resolution='c', llcrnrlon=0, llcrnrlat=-90,
        urcrnrlon=360.01, urcrnrlat=90)
x, y = m(*np.meshgrid(lons_c, lats))
clevs = [-120, -100, -80, -60, -40, -20, 0, 20, 40, 60, 80, 100, 120]
m.contourf(x, y, sf_200*1e-06, clevs, cmap=plt.cm.RdBu_r,
        extend='both')
m.drawcoastlines()
m.drawparallels((-90, -60, -30, 0, 30, 60, 90), labels=[1,0,0,0])
m.drawmeridians((0, 60, 120, 180, 240, 300, 360), labels=[0,0,0,1])
plt.colorbar(orientation='horizontal')
plt.title('Streamfunction ($10^6$m$^2$s$^{-1}$)', fontsize=16)