def free_energy(centers, A, levels=None, norm=None,\
                fmt='%.1f', method='linear', fill_value=np.nan,
                ax = None):
    r"""Make contourplot of alanine-dipeptide free energy.

    The scattered data is interpolated onto a regular grid 
    before plotting.

    Parameters
    ----------
    centers : (N, 2) ndarray 
        (phi, psi) coordinates of MSM discretization.
    A : (N, ) ndarray,
        Free energy.

    ax : optional matplotlib axis to plot to
    """
    X, Y=np.meshgrid(xcenters, ycenters)
    Z=griddata(centers, A, (X, Y), method=method, fill_value=fill_value)
    Z=Z-Z.min()
    if levels is None:
        levels=np.linspace(0.0, 50.0, 10)
    V=np.asarray(levels)
    if ax is None:
        fig=plt.figure()
        ax=fig.add_subplot(111)
    ax.set_xlim(-180.0, 180.0)
    ax.set_ylim(-180.0, 180.0)
    ax.set_xticks(np.linspace(-180.0, 180.0, 11))
    ax.set_yticks(np.linspace(-180.0, 180.0, 11))
    ax.set_xlabel(r"$\phi$")
    ax.set_ylabel(r"$\psi$")
    cs=ax.contour(X, Y, Z, V, norm=norm)
    plt.clabel(cs, inline=1, fmt=fmt)
    plt.grid()

def Pcolor(xs, ys, zs, pcolor=True, contour=False, **options):
    """Makes a pseudocolor plot.
    
    xs:
    ys:
    zs:
    pcolor: boolean, whether to make a pseudocolor plot
    contour: boolean, whether to make a contour plot
    options: keyword args passed to pyplot.pcolor and/or pyplot.contour
    """
    Underride(options, linewidth=3, cmap=matplotlib.cm.Blues)

    X, Y = np.meshgrid(xs, ys)
    Z = zs

    x_formatter = matplotlib.ticker.ScalarFormatter(useOffset=False)
    axes = pyplot.gca()
    axes.xaxis.set_major_formatter(x_formatter)

    if pcolor:
        pyplot.pcolormesh(X, Y, Z, **options)

    if contour:
        cs = pyplot.contour(X, Y, Z, **options)
        pyplot.clabel(cs, inline=1, fontsize=10)

def plot(variables,prev_vars,pltenv):

    cont_int = 10
    cont_smooth = 0.5
    thin = 10
    lvl = 2
    time = 0

    x = pltenv['x']
    y = pltenv['y']
    m = pltenv['map']


    bbox = dict(boxstyle="square",ec='None',fc=(1,1,1,0.75))

  
    T= variables['T_PL'][time][lvl]
    temp = T - 273.15 # convert to Celsius 
    levels = np.arange(-100,150,5)
    levels2 = np.arange(-50,50,1)

    P = m.contour(x,y,temp,levels=levels,colors='k')
    plt.clabel(P,inline=1,fontsize=10,fmt='%1.0f',inline_spacing=1)

    P = m.contour(x,y,temp,levels=[0],colors='r')
    plt.clabel(P,inline=1,fontsize=10,fmt='%1.0f',inline_spacing=1)


    m.contourf(x,y,temp, levels=levels2, extend='both')

def drawContour(img):
    contour_image = plt.contour(img, 5)

    plt.clabel(contour_image, inline=1, fontsize=10)
    plt.show()
    
    return contour_image

def plot_fgmax_grid():

    fg = fgmax_tools.FGmaxGrid()
    fg.read_input_data('fgmax_grid1.txt')
    fg.read_output()

    #clines_zeta = [0.01] + list(numpy.linspace(0.05,0.3,6)) + [0.5,1.0,10.0]
    clines_zeta = [0.001] + list(numpy.linspace(0.05,0.25,10))
    colors = geoplot.discrete_cmap_1(clines_zeta)
    plt.figure(1)
    plt.clf()
    zeta = numpy.where(fg.B>0, fg.h, fg.h+fg.B)   # surface elevation in ocean
    plt.contourf(fg.X,fg.Y,zeta,clines_zeta,colors=colors)
    plt.colorbar()
    plt.contour(fg.X,fg.Y,fg.B,[0.],colors='k')  # coastline

    # plot arrival time contours and label:
    arrival_t = fg.arrival_time/3600.  # arrival time in hours
    #clines_t = numpy.linspace(0,8,17)  # hours
    clines_t = numpy.linspace(0,2,5)  # hours
    #clines_t_label = clines_t[::2]  # which ones to label 
    clines_t_label = clines_t[::1]  # which ones to label 
    clines_t_colors = ([.5,.5,.5],)
    con_t = plt.contour(fg.X,fg.Y,arrival_t, clines_t,colors=clines_t_colors) 
    plt.clabel(con_t, clines_t_label)

    # fix axes:
    plt.ticklabel_format(format='plain',useOffset=False)
    plt.xticks(rotation=20)
    plt.gca().set_aspect(1./numpy.cos(fg.Y.mean()*numpy.pi/180.))
    plt.title("Maximum amplitude / arrival times (hrs)")

def animation_plot(i, pressure, wind, direction, step=24):
    """
    Function to update the animation frame.

    """
    # Clear figure to refresh colorbar
    plt.gcf().clf()

    ax = plt.axes(projection=ccrs.Mercator())
    ax.set_extent([-10.5, 3.8, 48.3, 60.5], crs=ccrs.Geodetic())

    contour_wind = iplt.contourf(wind[i][::10, ::10], cmap='YlGnBu',
                                 levels=range(0, 31, 5))
    contour_press = iplt.contour(pressure[i][::10, ::10], colors='white',
                                 linewidth=1.25, levels=range(938, 1064, 4))
    plt.clabel(contour_press, inline=1, fontsize=14, fmt='%i', colors='white')
    quiver_plot(wind[i], direction[i], step)
    plt.gca().coastlines(resolution='50m')

    time_points = pressure[i].coord('time').points
    time = str(pressure[i].coord('time').units.num2date(time_points)[0])
    plt.gca().set_title(time)

    colorbar = plt.colorbar(contour_wind)
    colorbar.ax.set_ylabel('Wind speed ({})'.format(str(wind[i].units)))

    def test(self):
        xlim=[-2.0,2.0]
        ylim=[-2.0,2.0]
        
        resol = 0.025
        sample_x = np.arange(-2.0,2.0,resol)
        sample_y = np.arange(-2.0,2.0,resol)
        sample_X, sample_Y = np.meshgrid(sample_x, sample_y)
        sample_XX = np.hstack([sample_X.reshape(sample_X.size,1),sample_Y.reshape(sample_Y.size,1)])

        sample_Z1 = np.exp(self.gmm_0.score(sample_XX))
        sample_Z2 = np.exp(self.gmm_1.score(sample_XX))

        plt.figure()
        ax1 = plt.subplot(121, aspect='equal')    
        CS = plt.contour(sample_X, sample_Y, sample_Z1.reshape(sample_X.shape))
        plt.clabel(CS, inline=1, fontsize=10)

        ax2 = plt.subplot(122, aspect='equal')    
        CS = plt.contour(sample_X, sample_Y, sample_Z2.reshape(sample_X.shape))
        plt.clabel(CS, inline=1, fontsize=10)

        ax1.set_xlim(xlim)
        ax1.set_ylim(ylim)

        ax2.set_xlim(xlim)
        ax2.set_ylim(ylim)

        plt.show()

        #print gmm_1.get_params(deep=True)
        print self.gmm_0.weights_
        print self.gmm_0.means_
        print self.gmm_0.covars_

def contour(x,y,z, linewidth = 2, labels = None):
    """
    Plots contours for non-evenly spaced data.
    x,y,z must be 1d arrays.
    lines = # of contour lines (default 18 )
    linewidth = line width of lines (default 2 )
    """
    
    assert x.shape[0] == y.shape[0] == z.shape[0], "arrays x,y,z must be the same size"
    
    #make a grid that surrounds x,y support
    xi = np.linspace(x.min(),x.max(),100)
    yi = np.linspace(y.min(),y.max(),100)
    # grid the data.
    zi = griddata((x, y), z, (xi[None,:], yi[:,None]), method='cubic')
    # contour the gridded data, plotting dots at the randomly spaced data points.
    plt.figure()
    CS = plt.contour(xi,yi,zi,linewidth=2)
    plt.clabel(CS, inline=1, fontsize=10)
    
    if labels:
        plt.xlabel(labels[0])
        plt.ylabel(labels[1])
    # plot data points.
    plt.scatter(x,y,c=z,s=60, alpha = 0.7, edgecolors = "none")
    plt.xlim(x.min(),x.max())
    plt.ylim(y.min(),y.max())
    plt.show()

    def draw(self):

        filename = self.filename
        file = open(os.getcwd() + "\\" + filename, 'r')
        lines = csv.reader(file)
        #
        data = []
        x = []
        y = []
        z = []
        for line in lines:
            try:
                data.append(line)
            except Exception as e:
                print e
                pass
        # print data
        for i in range(1, len(data)):
            try:
                x.append(float(data[i][0]))
                y.append(float(data[i][1]))
                z.append(float(data[i][3]))
            finally:
                pass

        xx = np.array(x)
        yy = np.array(y)
        zz = np.array(z)
        # print np.min(xx)

        tx = np.linspace(np.min(xx), np.max(xx), 100)
        ty = np.linspace(np.min(yy), np.max(yy), 100)

        XI, YI = np.meshgrid(tx, ty)

        rbf = interpolate.Rbf(xx, yy, zz, epsilon=2)
        ZI = rbf(XI, YI)

        #

        plt.gca().set_aspect(1.0)

        font = font_manager.FontProperties(family='times new roman', style='italic', size=16)

        cs = plt.contour(XI, YI, ZI, colors="black")
        plt.clabel(cs, cs.levels, inline=True, fontsize=10, prop=font)

        plt.subplot(1, 1, 1)
        plt.pcolor(XI, YI, ZI, cmap=cm.jet)
        plt.scatter(xx, yy, 100, zz, cmap=cm.jet)


        plt.title('interpolation example')
        plt.xlim(int(xx.min()), int(xx.max()))
        plt.ylim(int(yy.min()), int(yy.max()))
        plt.colorbar()
        plt.savefig("interpolation.jpg")
        #plt.show()

        return ZI, XI, YI

def main():
    # Lets relax a 1D array of values. 
    N = 20
    a = 0.0
    b = 1.0
    u = np.zeros(20)
    u[0] = a
    u[-1] = b
    x = np.arange(0.0, float(len(u)), 1.0) / float(len(u) - 1)
    u = relax1D(u)
    print "It took %d iterations to reach the desired tolerance." % (len(u) - 1)
    animate1Dframes(x, u)
    plt.plot(x, u[-1], "-sk")
    plt.savefig("Relaxed1D.png")
    plt.close()

    # Now lets relax a 2D array of values. Such that the top and bottom edges are set to zero
    # and the left and right edges are set to one. 
    u2D = np.zeros((N, N))
    u2D[0,:]  = a
    u2D[-1,:] = a
    u2D[:,0]  = b
    u2D[:,-1] = b
    u2D = relax2D(u2D)
    print "It took %d iterations to reach the desired tolerance." % (len(u2D) - 1)
    animate2Dframes(u2D)
    u2D_cf = plt.contour(u2D[-1], levels=np.arange(0, 1, 0.1))
    plt.clabel(u2D_cf, colors='k')
    plt.savefig("Relaxed2D.png")
    plt.close()

def plot_pressuremap(data, 
                    title='pressure pattern', 
                    sub_title='ploted in birdhouse'):
  """
  plots pressure data
  :param data: 2D or 3D array of pressure data. if data == 3D a mean will be calculated
  :param title: string for title
  :param sub_title: string for sub_title
  """
  from numpy import squeeze, mean
  
  d = squeeze(data)

  if len(d.shape)==3:
    d = mean(d, axis=0)
  if len(d.shape)!=2:
    logger.error('data are not in shape for map display')

  co = plt.contour(d, lw=2, c='black')
  cf = plt.contourf(d)

  plt.colorbar(cf)
  plt.clabel(co, inline=1) # fontsize=10
  
  plt.title(title)
  plt.annotate(sub_title, (0,0), (0, -30), xycoords='axes fraction',
               textcoords='offset points', va='top')

  ip, image = mkstemp(dir='.',suffix='.png')
  plt.savefig(image)
  plt.close()

  return image

def plot_area(dataset, request, area_prediction):
    bot_lat, top_lat, left_lon, right_lon = request['predict_area']

    plt.figure(figsize=(20, 15))

    # Plot stations
    sp = plt.scatter(
        dataset.longitude, dataset.latitude,
        s=10, c=dataset['TT2m_error'],
        edgecolor='face',
        vmin=-5, vmax=5,
    )

    # Contours
    contour_handle = plt.contour(
        area_prediction,
        np.arange(-5, 5, 1),
        antialiased=True,
        extent=(left_lon, right_lon, top_lat, bot_lat),
        zorder=999,
        alpha=0.5
    )
    plt.clabel(contour_handle, fontsize=11)

    # Color bar
    cb = plt.colorbar(sp)
    cb.set_label('Temperature error')
    plt.show()

def main():
    dir='/Users/ph290/Public/mo_data/ostia/' # on ELD140
    filename = dir + '*.nc'

    cube = iris.load_cube(filename,'sea_surface_temperature',callback=my_callback)
    #reads in data using a special callback, because it is a nasty netcdf file
    cube.data=cube.data-273.15


    sst_mean = cube.collapsed('time', iris.analysis.MEAN)
    #average all 12 months together
    sst_stdev=cube.collapsed('time', iris.analysis.STD_DEV)

    caribbean = iris.Constraint(
                                    longitude=lambda v: 260 <= v <= 320,
                                    latitude=lambda v: 0 <= v <= 40,
                                    name='sea_surface_temperature'
                                    )

    caribbean_sst_mean = sst_mean.extract(caribbean)
    caribbean_sst_stdev = sst_stdev.extract(caribbean)
    #extract the Caribbean region

    fig = plt.figure()
    ax = plt.axes(projection=ccrs.PlateCarree())
    data=caribbean_sst_mean.data
    data2=caribbean_sst_stdev.data
    lons = caribbean_sst_mean.coord('longitude').points
    lats = caribbean_sst_mean.coord('latitude').points
    lo = np.floor(data.min())
    hi = np.ceil(data.max())
    levels = np.linspace(lo,hi,100)
    lo2 = np.floor(data2.min())
    hi2 = np.ceil(data2.max())
    levels2 = np.linspace(lo2,5,10)
    cube_label = 'latitude: %s' % caribbean_sst_mean.coord('latitude').points
    contour=plt.contourf(lons, lats, data,transform=ccrs.PlateCarree(),levels=levels,xlabel=cube_label)
    #filled contour the annually averaged temperature
    contour2=plt.contour(lons, lats, data2,transform=ccrs.PlateCarree(),levels=levels2,colors='k')
    #contour the standard deviations
    plt.clabel(contour2, inline=0.5, fontsize=12,fmt='%1.1f' )
    ax.add_feature(cartopy.feature.LAND)
    ax.coastlines()
    ax.add_feature(cartopy.feature.RIVERS)
    ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
    #ax.add_feature(cartopy.feature.LAKES, alpha=0.5)
    cbar = plt.colorbar(contour, ticks=np.linspace(lo,hi,7), orientation='horizontal')
    cbar.set_label('Sea Surface Temperature ($^\circ$C)')
    # enable axis ticks
    ax.xaxis.set_visible(True)
    ax.yaxis.set_visible(True)
    # fix 10-degree increments and don't clip range
    plt.locator_params(steps=(1,10), tight=False)
    # add gridlines
    plt.grid(True)
    # add axis labels
    plt.xlabel('longitude')
    plt.ylabel('latitude')
    #plt.show()
    plt.savefig('/home/h04/hador/public_html/twiki_figures/carib_sst_and_stdev.png')

def eigenvector(centers, ev, levels=None, norm=None,\
                    fmt='%.e', method='linear', fill_value=np.nan):
    r"""Make contourplot of alanine-dipeptide stationary distribution

    The scattered data is interpolated onto a regular grid before
    plotting.

    Parameters
    ----------
    centers : (N, 2) ndarray 
        (phi, psi) coordinates of MSM discretization.
    ev : (N, ) ndarray,
        Right eigenvector
    
    """  
    X, Y=np.meshgrid(xcenters, ycenters)
    Z=griddata(centers, ev, (X, Y), method=method, fill_value=fill_value)
    if levels is None:
        levels=np.linspace(-0.1, 0.1, 21)
    V=np.asarray(levels)    
    fig=plt.figure()
    ax=fig.add_subplot(111)
    ax.set_xlim(-180.0, 180.0)
    ax.set_ylim(-180.0, 180.0)
    ax.set_xticks(np.linspace(-180.0, 180.0, 11))
    ax.set_yticks(np.linspace(-180.0, 180.0, 11))
    ax.set_xlabel(r"$\phi$")
    ax.set_ylabel(r"$\psi$")
    cs=ax.contour(X, Y, Z, V, norm=norm)
    plt.clabel(cs, inline=1, fmt=fmt)
    plt.grid()

def __make_result_sps(session):
    sdf_extension = session['SPS']
    # if isinstance(sdf_extension, SDFExtension) == False:
    #    return sdf_extension
    scales = sdf_extension.actual_scales
    zs = sdf_extension.result_datasets
    x, y = numpy.meshgrid(scales[1].data, scales[0].data)
    response_string = ''
    for zds in zs:
        if zds.name != 'summary.batteryPower' and zds.name != 'summarySystemSPS.mechanicalPower':
            continue
        print zds.data
        z = numpy.fabs(zds.data)
        fig = plt.contourf(x, y, z, 10, alpha=.75, cmap=cm.coolwarm)
        C = plt.contour(x, y, z, 10, colors='black', linewidth=.5)
        plt.clabel(C, inline=1, fontsize=10)
        # TODO: Unit hackcode
        plt.xlabel(scales[1].name + ' (' + str(scales[1].unit) + ')')
        plt.ylabel(scales[0].name + ' (' + str(scales[0].unit) + ')')
        plt.title(zds.name)
        plt.colorbar(fig, shrink=0.5, aspect=5)
        plt.savefig('temp.svg')

        with open('temp.svg', 'rb') as f:
            data = f.read().encode('base64')
        os.remove('temp.svg')
        response_string += '<img src="data:image/svg+xml;base64,{0}"><br/>'.format(data)
        plt.close('all')
    return response_string

def plot_TS(temp, psal, depth, lon, lat, svec, tvec, density, title, m, figname):
    '''
    Create the T-S diagram
    '''
    logger = logging.getLogger(__name__)
    fig = plt.figure(figsize=(15, 15))
    rcParams.update({'font.size': 18})
    plt.scatter(psal, temp, s=5, c=depth, vmin=10., vmax=1000.,
               edgecolor='None', cmap=plt.cm.plasma)
    cbar = plt.colorbar(extend='max')
    plt.xlabel('Salinity', fontsize=18)
    plt.ylabel('Temperature\n($^{\circ}$C)', rotation=0, ha='right', fontsize=18)
    cont = plt.contour(svec, tvec, density, levels=np.arange(22., 32., 1.), 
                       colors='.65', linestyles='dashed', lineswidth=0.5)
    plt.clabel(cont,inline=True, fmt='%1.1f')
    plt.xlim(smin, smax)
    plt.ylim(tmin, tmax)
    cbar.set_label('Depth\n(m)', rotation=0, ha='left')
    plt.grid(color="0.6")

    # Add an inset showing the positions of the platform
    inset=plt.axes([0.135, 0.625, 0.3, 0.35])
    lon2plot, lat2plot = m(lon, lat)
    m.drawmapboundary(color='w')
    m.plot(lon2plot, lat2plot, 'ro', ms=2, markeredgecolor='r')
    #m.drawcoastlines(linewidth=0.25)
    m.drawlsmask(land_color='0.4', ocean_color='0.9', lakes=False)
    plt.title(title, fontsize=20)
    plt.savefig(figname, dpi=150)
    # plt.show()
    plt.close()

   def publish(self):
      alldata = []
      
      for i,val in enumerate(self.status['runs']):
         print val
         alldata.append(post_bout.read(path=val))
     
      alldata = np.array(alldata)
      print "alldata.shape: ", alldata.shape
      print  alldata[0]['Ni']['ave']['amp'].shape
      data = alldata[0]['Ni']['ave']['amp'] #Nt long array
      pp = PdfPages('output.pdf')
      plt.figure()
      cs = plt.plot(data)
      plt.title('amplitude ')
      plt.savefig(pp, format='pdf')

      gamma = alldata[2]['Ni']['modes'][0]['gamma'] #single value
      phase = alldata[2]['Ni']['modes'][0]['phase'] #Nt x Nx array 

      plt.figure()
      cs = plt.contour(phase)
      plt.clabel(cs, inline=1, fontsize=10)
      plt.title('phase')
      plt.savefig(pp, format='pdf')

      plt.close() 
      pp.close()
   
      explain = alldata[0]['meta']
      print explain
      
      allpickle = open('allpickled.pkl','wb')
      pickle.dump(alldata,allpickle)
      allpickle.close()

def plotSNR(R, snr, T0, fname = '/Users/jmrv/Desktop/snr.pdf'):
	pp = PP(fname)
	plt.clf()

	T = T0*10**(np.arange(-2, 1.5, 0.1))
	scale = np.sqrt(T/T0)

	snrGrid = np.zeros((len(R), len(T)))
	for i in range(len(snr)):
		snrGrid[i,:] = scale*snr[i]

	levels = np.array([0.5, 1, 2, 3, 5, 10, 15, 20])/10
	c = plt.contour(T/1e3, R*1000, snrGrid, levels, colors='k')
	plt.clabel(c, fontsize=9, inline=1)
	plt.xscale('log')
	plt.yscale('log')
	plt.title('Signal-to-noise of Sc signal in the X-ray band')
	plt.ylabel('Instrument resolution (eV, FWHM)')
	plt.xlabel('Exposure (ks)')
	#plt.ticklabel_format(style='plain', axis='both')
	
	pp.savefig()
	pp.close()

	print '\nplotted '+fname

def kden_plot(X,Y, im, support, pvals=(.1,.333,.5,.666,.9),
              contours=True, imshow=True, cbar=False, psupport=True,oned=True,
              fignum=None,cfontsize=10):
    if fignum is None:
        fignum=plt.gcf().number
        plt.clf()
    
    dx=np.abs(X[0,0]-X[1,0])
    dy=np.abs(Y[0,0]-Y[0,1])
    dxdy=dx*dy
    if contours:
        
        levelsdict={x:'{:3.0f}%'.format(p*100)
                    for x,p in find_prob_contours(im,dxdy,pvals=pvals)}
        cplt=plt.contour(X,Y,im,colors='k', levels=levelsdict.keys())
        plt.clabel(cplt, inline=1, fontsize=10, use_clabeltext=True,
                   fmt=levelsdict)
    if imshow:
        plt.imshow(np.rot90(im),aspect='auto',
                   extent=[X.min(), X.max(), Y.min(), Y.max()],
                   cmap=plt.cm.gist_earth_r)
    if cbar:
        plt.colorbar()

    if psupport:
        plt.plot(support[0],support[1],'r,')

    if oned:
        plt.figure(fignum+1)
        plt.plot(X[:,0],((im*dy).sum(1)))

def plot_image(file_in, num=0, xlab='', ylab='', map_range=[0, 0, 0, 0],
               color_loc=False, cm='', fig_title='', cont=''):
    """
    plot a given 2d image
    """
    from matplotlib.pyplot import imshow, colorbar, figure
    from matplotlib.pyplot import xlabel, ylabel
    from matplotlib.pyplot import title, contour, clabel
    figure(num)

    if fig_title != '':
        title(fig_title)

    if xlab != '':
        xlabel(xlab)

    if ylab != '':
        ylabel(ylab)

    if all([v == 0 for v in map_range]):
        map_range = [0, len(file_in), 0, len(file_in)]

    if cm == '':
        imshow(file_in, extent=map_range, origin='lower', interpolation='None')
    else:
        imshow(file_in, extent=map_range, origin='lower',
               interpolation='None', cmap=cm)
    if color_loc:
        colorbar()
    if cont != '':
        CS = contour(file_in, cont, extent=map_range, origin='lower',
                     hold='on', colors=('k', ))
        clabel(CS, inline=1, fmt='%2.1f', fontsize=10)