def normalize_colors(vmin, vmax, clip=False):
    """Helper to handle matplotlib API"""
    import matplotlib.pyplot as plt
    try:
        return plt.Normalize(vmin, vmax, clip=clip)
    except AttributeError:
        return plt.normalize(vmin, vmax, clip=clip)

def normalize_colors(vmin, vmax, clip=False):
    """Helper to handle matplotlib API"""
    import matplotlib.pyplot as plt
    if 'Normalize' in vars(plt):
        return plt.Normalize(vmin, vmax, clip=clip)
    else:
        return plt.normalize(vmin, vmax, clip=clip)

def plot_veg_types(yc, xc, Cv, baresoil):
    Projection_Parameters = projection

    labels = ['Evergreen Needleleaf', 'Evergreen Broadleaf',
              'Deciduous Needleleaf', 'Deciduous Broadleaf', 'Mixed Cover',
              'Woodland', 'Wooded Grasslands', 'Closed Shrublands',
              'Open Shrublands', 'Grasslands', 'Crop land', 'Bare Soil/Ice']

    fig = plt.figure(figsize=(10, 12))

    gs1 = gridspec.GridSpec(4, 3)

    for loc in xrange(11):
        ax = fig.add_subplot(gs1[loc])
        c = plot_map(ax, yc, xc, Cv[loc], Projection_Parameters, vmin=0,
                     cmap='Jet')
        ax.set_title(labels[loc])
    ax = fig.add_subplot(gs1[11])
    c = plot_map(ax, yc, xc, baresoil, Projection_Parameters, vmin=0,
                 cmap='Jet')
    ax.set_title(labels[11])

    sm = plt.cm.ScalarMappable(cmap='Jet', norm=plt.normalize(vmin=0, vmax=1))
    colorbar_ax = fig.add_axes([0.92, 0.1, 0.03, 0.8])
    sm._A = []
    plt.colorbar(sm, cax=colorbar_ax)
    fig.suptitle('Fraction of Vegetation Type', fontsize=20, fontweight='bold')
    fig.text(0.5, 0.93, 'Regional Arctic Climate Model',
             ha='center', fontsize=18)

    plt.show()

def buildHeatmap(dataDictionary,title,bartitle,barcolor,barmin,barmax):
    # Lambert Conformal map of USA lower 48 states
    geoMap = Basemap(llcrnrlon=-119, llcrnrlat=22, urcrnrlon=-64,
    urcrnrlat=49, projection='lcc', lat_1=33, lat_2=45,
    lon_0=-95, resolution='i', area_thresh=10000)
    # Uses a Shape File from the US Census Bureau
    # http://www2.census.gov/geo/tiger/GENZ2010/
    shapeInfo = geoMap.readshapefile('data/shape_files/gz_2010_us_040_00_500k','states',drawbounds=False)
    
    # Process based off an example from our professor's lecture slides,
    # James Bagrow - UVM Spring 2014
    # http://bagrow.com/dsv/
    # I took the base idea and modified it heavily for the needs 
    # of this project.
    # Color each state based on Population
    colors={}
    statenames=[]
    # Use the Green color map.
    # Others may be found here : http://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps
    colorMap = barcolor
    heatmapMinAmt = barmin; heatmapMaxAmt = barmax # Range for Heatmap (Valuewise)
    scalarMapping = plt.cm.ScalarMappable(cmap=colorMap,norm=plt.normalize(vmin=heatmapMinAmt, vmax=heatmapMaxAmt))
    
    for shapedict in geoMap.states_info:
        statename = shapedict['NAME']
        # Grab population info for each state
        try:
            population = float(dataDictionary[statename])
        except KeyError:
            population = 0.0
            
        # Define the color for this state.
        # Population / Max Color Amt.
        colors[statename] = colorMap((population-heatmapMinAmt)/(heatmapMaxAmt-heatmapMinAmt))
        statenames.append(statename)
    
    # Actually color each state in this structure.
    for nshape, seg in enumerate(geoMap.states):
        xx, yy = zip(*seg)
        color = rgb2hex(colors[statenames[nshape]]) 
        plt.fill(xx,yy,color,edgecolor=color)
    
    # Bound our map to the continental US.
    # otherwise, we'll have some weird display issues
    # with Alaska and Hawaii.
    geoMap.drawparallels(np.arange(25,65,20),labels=[0,0,0,0],zorder=-1,color="w")
    geoMap.drawmeridians(np.arange(-120,-40,20),labels=[0,0,0,0],zorder=-1,color="w")
    
    # Build the Color Bar at the bottom of the graph
    mm = plt.cm.ScalarMappable(cmap=colorMap)
    mm.set_array([heatmapMinAmt,heatmapMaxAmt])
    plt.colorbar(mm, label=bartitle,ticks=[0,5000,10000,15000,20000,25000],
    orientation="horizontal",fraction=0.05)
                
    # Adjust size of image
    plt.gcf().set_size_inches(12.0,8.0)
    plt.gca().axis("off")
    # Add Title and show.
    plt.title(title)
    plt.show()

def showOverlapTable(modes_x, modes_y, **kwargs):
    """Show overlap table using :func:`~matplotlib.pyplot.pcolor`.  *modes_x*
    and *modes_y* are sets of normal modes, and correspond to x and y axes of
    the plot.  Note that mode indices are incremented by 1.  List of modes
    is assumed to contain a set of contiguous modes from the same model.

    Default arguments for :func:`~matplotlib.pyplot.pcolor`:

      * ``cmap=plt.cm.jet``
      * ``norm=plt.normalize(0, 1)``"""

    import matplotlib.pyplot as plt

    overlap = abs(calcOverlap(modes_y, modes_x))
    if overlap.ndim == 0:
        overlap = np.array([[overlap]])
    elif overlap.ndim == 1:
        overlap = overlap.reshape((modes_y.numModes(), modes_x.numModes()))

    cmap = kwargs.pop('cmap', plt.cm.jet)
    norm = kwargs.pop('norm', plt.normalize(0, 1))
    show = (plt.pcolor(overlap, cmap=cmap, norm=norm, **kwargs),
            plt.colorbar())
    x_range = np.arange(1, modes_x.numModes() + 1)
    plt.xticks(x_range-0.5, x_range)
    plt.xlabel(str(modes_x))
    y_range = np.arange(1, modes_y.numModes() + 1)
    plt.yticks(y_range-0.5, y_range)
    plt.ylabel(str(modes_y))
    plt.axis([0, modes_x.numModes(), 0, modes_y.numModes()])
    if SETTINGS['auto_show']:
        showFigure()
    return show

def discrete_calldata_colormesh(X, labels=None, colors='wbgrcmyk', states=None, ax=None, **kwargs):
    """
    Make a meshgrid from discrete calldata (e.g., genotypes).

    Parameters
    ----------

    X: array
        2-dimensional array of integers of shape (#variants, #samples)
    labels: sequence of strings
        Axis labels (e.g., sample IDs)
    colors: sequence
        Colors to use for different values of the array
    states: sequence
        Manually specify discrete calldata states (if not given will be determined from the data)
    ax: axes
        Axes on which to draw

    Remaining keyword arguments are passed to ax.pcolormesh.

    """

    # set up axes
    if ax is None:
        fig = plt.figure(figsize=(7, 5))
        ax = fig.add_subplot(111)

    # determine discrete states
    if states is None:
        states = np.unique(X)
    colors = colors[:max(states)-min(states)+1]  # only need as many colors as states

    # plotting defaults
    pltargs = {
        'cmap': ListedColormap(colors),
        'norm': plt.normalize(min(states), max(states)+1),
    }
    pltargs.update(kwargs)

    ax.pcolormesh(X.T, **pltargs)
    ax.set_xlim(0, X.shape[0])
    ax.set_ylim(0, X.shape[1])

    ax.set_yticks(np.arange(X.shape[1]) + .5)
    if labels is not None:
#        labels = ['%s [%s] ' % (s, i) for (i, s) in enumerate(labels)]
        ax.set_yticklabels(labels, rotation=0)

    return ax

    def coRes(self):
        layers = self.iface.legendInterface().layers()
        for ly in layers:
            if ly.name() == 'r1d_in_data':
                vl2 = ly
        print vl2.name()
        x=[]
        y=[]
        values =[]
        ab = self.dockwidget.txtAB2CoRes.text()
        print ab

        for f in vl2.getFeatures():
            #print f['ab2']
            if f['ab2'] == float(ab):
                geom = f.geometry()
                x.append(geom.asPoint().x())
                y.append(geom.asPoint().y())
                values.append(f['ra'])

        #Creating the output grid (100x100, in the example)
        #x = [10,60,40,70,10,50,20,70,30,60]
        #y = [10,20,30,30,40,50,60,70,80,90]
        #values = [1,2,2,3,4,6,7,7,8,10]
        print x
        print y
        #print xx

        txi = np.linspace(min(x), max(x), 10000)
        tyi = np.linspace(min(y), max(y), 10000)
        XI, YI = np.meshgrid(ti, ti)

        #Creating the interpolation function and populating the output matrix value
        rbf = Rbf(x, y, values, function='inverse')
        ZI = rbf(XI, YI)

        # Plotting the result
        n = plt.normalize(0.0, 1000.0)
        plt.subplot(1, 1, 1)
        plt.pcolor(XI, YI, ZI)
        plt.scatter(x, y, 1000, values)
        plt.title('RBF interpolation')
        plt.xlim(min(x), max(x))
        plt.ylim(min(y), max(y))
        plt.colorbar()

        plt.show()

def showOverlapTable(rows, cols, **kwargs):
    """Show overlap table using :func:`~matplotlib.pyplot.pcolor`.  *rows* and 
    *cols* are sets of normal modes, and correspond to rows and columns of the 
    displayed overlap matrix.  Note that mode indices are incremented by 1.  
    List of modes should contain a set of contiguous modes from the same model. 
    
    Default arguments for :func:`~matplotlib.pyplot.pcolor`:
        
      * ``cmap=plt.cm.jet``
      * ``norm=plt.normalize(0, 1)``
    
    .. plot::
       :context:
       :include-source:
        
       plt.figure(figsize=(5,4))
       showOverlapTable( p38_pca[:6], p38_anm[:6] )
       plt.title('p38 PCA vs ANM')

    .. plot::
       :context:
       :nofigs:
        
       plt.close('all')"""
    
    import matplotlib.pyplot as plt
    
    overlap = abs(calcOverlap(rows, cols))
    if isinstance(rows, NMA):
        rows = rows[:]
    if isinstance(cols, NMA):
        cols = cols[:]
    cmap = kwargs.pop('cmap', plt.cm.jet)
    norm = kwargs.pop('norm', plt.normalize(0, 1))
    show = (plt.pcolor(overlap, cmap=cmap, norm=norm, **kwargs),
            plt.colorbar())
    x_range = np.arange(1, len(cols)+1)
    plt.xticks(x_range-0.5, x_range)
    plt.xlabel(str(cols))
    y_range = np.arange(1, len(rows)+1)
    plt.yticks(y_range-0.5, y_range)
    plt.ylabel(str(rows))
    plt.axis([0, len(cols), 0, len(rows)])
    return show

def quick_ortho_plot(x,y,z, bbox=[], loc=None, title='', norm=None):
    import matplotlib.pyplot as pp
    import xipy.volume_utils as vu
    from xipy.vis.single_slice_plot import SliceFigure
    # find or make the x, y, z plane extents
    if bbox:
        extents = vu.limits_to_extents(bbox)
    else:
        extents = [ reduce(lambda x,y: x+y, zip([0,0],p.shape[:2][::-1]))
                    for p in (x, y, z) ]
        
    if norm is None and len(x.shape) < 3:
        mx = max([ x.max(), y.max(), z.max() ])
        mn = min([ x.min(), y.min(), z.min() ])
        norm = pp.normalize(mn, mx)
    fig = pp.figure()

    loc = list(loc)
    zloc = loc[0],loc[1] if loc else None
    ax_z = fig.add_subplot(131)
    sf_z = SliceFigure(fig, extents[2])
    img_z = sf_z.spawn_image(z, extent=extents[2], loc=zloc,
                             interpolation='nearest', norm=norm)
    ax_z.set_title('plot z')

    yloc = loc[0],loc[2] if loc else None
    ax_y = fig.add_subplot(132)
    sf_y = SliceFigure(fig, extents[1])
    img_y = sf_y.spawn_image(y, extent=extents[1], loc=yloc,
                             interpolation='nearest', norm=norm)
    ax_y.set_title('plot y')

    xloc = loc[1],loc[2] if loc else None
    ax_x = fig.add_subplot(133)
    sf_x = SliceFigure(fig, extents[0])
    img_x = sf_x.spawn_image(x, extent=extents[0], loc=xloc,
                             interpolation='nearest', norm=norm)
    ax_x.set_title('plot x')
    if title:
        fig.text(.5, .05, title, ha='center')
    pp.colorbar(img_x.img)
    pp.show()
    return fig

def draw():
    global m, learn, z
    pyp.background(200,50)

    learn.collect_data()
    m = learn.get_mean()


    pyp.loadPixels()
    m = np.atleast_2d(m)
    norm = normalize(vmin=min(min(m)), vmax=max(max(m)))
    cmap = get_cmap('jet')
    m_normed = norm(m)
    rgba_data=cmap(m_normed)*255
    r = rgba_data[0,:,0].astype('uint32')
    g = rgba_data[0,:,1].astype('uint32')
    b = rgba_data[0,:,2].astype('uint32')
    a = rgba_data[0,:,3].astype('uint32')
    pyp.screen.pixels = a << 24 | r << 16 | g << 8 | b
    pyp.updatePixels()

def find_image_threshold(arr, percentile=90., debug=False):
    nbins = 200
    bsizes, bpts = np.histogram(arr.flatten(), bins=nbins)
    # heuristically, this should show up near the middle of the
    # second peak of the intensity histogram
    start_pt = np.abs(bpts - arr.max()/2.).argmin()
    db = np.diff(bsizes[:start_pt])
##     zcross = np.argwhere((db[:-1] < 0) & (db[1:] >= 0)).flatten()[0]
    bval = bsizes[1:start_pt-1][ (db[:-1] < 0) & (db[1:] >= 0) ].min()
    zcross = np.argwhere(bval==bsizes).flatten()[0]
    thresh = (bpts[zcross] + bpts[zcross+1])/2.
    # interpolate the percentile value from the bin edges
    bin_lo = int(percentile * nbins / 100.0)
    bin_hi = int(round(percentile * nbins / 100.0 + 0.5))
    p_hi = percentile - bin_lo # proportion of hi bin
    p_lo = bin_hi - percentile # proportion of lo bin
##     print bin_hi, bin_lo, p_hi, p_lo
    pval = bpts[bin_lo] * p_lo + bpts[bin_hi] * p_hi
    if debug:
        import matplotlib as mpl
        import matplotlib.pyplot as pp
        f = pp.figure()
        ax = f.add_subplot(111)
        ax.hist(arr.flatten(), bins=nbins)
        l = mpl.lines.Line2D([thresh, thresh], [0, .25*bsizes.max()],
                             linewidth=2, color='r')
        ax.add_line(l)
        ax.xaxis.get_major_formatter().set_scientific(True)
        f = pp.figure()
        norm = pp.normalize(0, pval)
        ax = f.add_subplot(211)
        plot_arr = arr
        while len(plot_arr.shape) > 2:
            plot_arr = plot_arr[plot_arr.shape[0]/2]
        ax.imshow(plot_arr, norm=norm)
        ax = f.add_subplot(212)
        simple_mask = (plot_arr < thresh)
        ax.imshow(np.ma.masked_array(plot_arr, mask=simple_mask), norm=norm)
        pp.show()
    
    return thresh, pval

leg = ax.legend([l1, l2, l3], labels, ncol=3, frameon=False, fontsize=16, 
                bbox_to_anchor=[1.1, 0.11], handlelength=2, 
                handletextpad=1, columnspacing=2, title='Average Page Size')
 
# Customize legend title
# Set position to increase space between legend and labels
plt.setp(leg.get_title(), fontsize=20, alpha=a)
leg.get_title().set_position((0, 10))
# Customize transparency for legend labels
[plt.setp(label, alpha=a) for label in leg.get_texts()]

# Create a fake colorbar
ctb = LinearSegmentedColormap.from_list('custombar', customcmap, N=2048)
# Trick from http://stackoverflow.com/questions/8342549/
# matplotlib-add-colorbar-to-a-sequence-of-line-plots
sm = plt.cm.ScalarMappable(cmap=ctb, norm=plt.normalize(vmin=1, vmax=3))
# Fake up the array of the scalar mappable
sm._A = []
 
# Set colorbar, aspect ratio
cbar = plt.colorbar(sm, alpha=0.05, aspect=16, shrink=0.4)
cbar.solids.set_edgecolor("face")
# Remove colorbar container frame
cbar.outline.set_visible(False)
# Fontsize for colorbar ticklabels
cbar.ax.tick_params(labelsize=16)
# Customize colorbar tick labels
mytks = range(1,4,1)
cbar.set_ticks(mytks)
cbar.ax.set_yticklabels([str(i) for i in mytks], alpha=a)
 

def vectorize(hillshade_file, m_value_file):

    import matplotlib.pyplot as pp
    import numpy as np
    import matplotlib.colors as colors
    import matplotlib.cm as cmx
    from matplotlib import rcParams

    # get data
    hillshade, hillshade_header = read_flt(hillshade_file)
    m_values, m_values_header = read_flt(m_value_file)

    # handle plotting hillshades which are larger than the m_value raster
    # cannot cope with m_value raster larger than the hillshade
    corrected_x = 0
    corrected_y = 0
    if (hillshade_header[0] != m_values_header[0]) or (hillshade_header[1] != m_values_header[1]):
        corrected_x = (m_values_header[2] - hillshade_header[2]) / hillshade_header[4]
        corrected_y = (
            ((m_values_header[3] / m_values_header[4]) + m_values_header[1])
            - ((hillshade_header[3] / hillshade_header[4]) + hillshade_header[1])
        ) * -1

    # ignore nodata values
    hillshade = np.ma.masked_where(hillshade == -9999, hillshade)
    m_values = np.ma.masked_where(m_values == -9999, m_values)

    # fonts
    rcParams["font.family"] = "sans-serif"
    rcParams["font.sans-serif"] = ["arial"]
    rcParams["font.size"] = 12

    fig, ax = pp.subplots()

    ax.imshow(hillshade, vmin=0, vmax=255, cmap=cmx.gray)

    xlocs, xlabels = pp.xticks()
    ylocs, ylabels = pp.yticks()

    new_x_labels = np.linspace(
        hillshade_header[2], hillshade_header[2] + (hillshade_header[1] * hillshade_header[4]), len(xlocs)
    )
    new_y_labels = np.linspace(
        hillshade_header[3], hillshade_header[3] + (hillshade_header[0] * hillshade_header[4]), len(ylocs)
    )

    new_x_labels = [str(x).split(".")[0] for x in new_x_labels]  # get rid of decimal places in axis ticks
    new_y_labels = [str(y).split(".")[0] for y in new_y_labels][::-1]  # invert y axis
    pp.xticks(xlocs[1:-1], new_x_labels[1:-1], rotation=30)  # [1:-1] skips ticks where we have no data
    pp.yticks(ylocs[1:-1], new_y_labels[1:-1])

    pp.xlabel("Easting (m)")
    pp.ylabel("Northing (m)")

    # SET UP COLOURMAPS
    jet = pp.get_cmap("jet")

    m_MIN = np.min(m_values)
    m_MAX = np.max(m_values)
    cNorm_m_values = colors.Normalize(vmin=m_MIN, vmax=m_MAX)
    scalarMap_m_values = cmx.ScalarMappable(norm=cNorm_m_values, cmap=jet)

    for i in xrange(len(m_values)):
        for j in xrange(len(m_values[0])):
            if m_values[i][j] > 0:
                colorVal = scalarMap_m_values.to_rgba(m_values[i][j])
                pp.scatter(j + corrected_x, i + corrected_y, marker=".", color=colorVal, edgecolors="none")

    # Configure final plot
    sm = pp.cm.ScalarMappable(cmap=jet, norm=pp.normalize(vmin=m_MIN, vmax=m_MAX))
    sm._A = []
    cbar = pp.colorbar(sm)
    cbar.set_label("M Values")

    pp.show()

#.........这里部分代码省略.........
        if channel_id[i]==0:
            plt.plot(chi[i], m_mean[i], "o", markersize=10, color='black', markeredgecolor = 'black',alpha=0.7)  
        else:
            plt.plot(chi[i], m_mean[i], "o", markersize=8, color=colorVal, markeredgecolor = colorVal,alpha = 0.7)
                 
           # plt.plot(chi[i], m_mean[i] + m_standard_error[i], ".", linewidth=0.25, color='k')
           # plt.plot(chi[i], m_mean[i] - m_standard_error[i], ".", linewidth=0.25, color='k')
 
    rect1 =plt.Rectangle((21.4577,0), 18.2625, 20, color='red', alpha= 0.1)
    rect2 =plt.Rectangle((69.0736,0), 29.3623, 20, color='green', alpha= 0.1)
    ax.add_artist(rect1)
    ax.add_artist(rect2)

           
    # Configure final plot
    ax.spines['top'].set_linewidth(2.5)
    ax.spines['left'].set_linewidth(2.5)
    ax.spines['right'].set_linewidth(2.5)
    ax.spines['bottom'].set_linewidth(2.5) 
    ax.tick_params(axis='both', width=2.5)    
    plt.xlabel('$\chi$ (m)', fontsize = axis_size)
    plt.ylabel('Gradient in $\chi$ space', fontsize = axis_size)
    plt.title('$A_0$: '+str(A_0)+' m$^2$, and $m/n$: '+str(m_over_n), fontsize = label_size)

    #=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    # LONGITUDINAL PROFILES WITH COLOUR SCALE GIVING CHI-m VALUE
    plt.figure(3, facecolor='white',figsize=(10,7.5))  
    ax = plt.subplot(1,1,1)    
    for i in range (0,len(channel_id)):
            colorVal = scalarMap_m_values.to_rgba(m_mean[i])
            plt.plot((flow_dist[i]-minfd)/1000, elevation[i], "o", markersize=6, color=colorVal, markeredgecolor = colorVal,alpha = 0.5)       
                  
    # Configure final plot
    sm = plt.cm.ScalarMappable(cmap=hot,norm=plt.normalize(vmin=np.min(m_mean), vmax=np.max(m_mean)))
    sm._A = []
    
    cbar = plt.colorbar(sm,orientation='horizontal',use_gridspec=True)
    cbar.set_label('Gradient in $\chi$ space', fontsize = axis_size)
    plt.xlabel('Distance upstream (km)', fontsize = axis_size)
    plt.ylabel('Elevation (m)', fontsize = axis_size)
    plt.title('$A_0$: '+str(A_0)+' m$^2$, and $m/n$: '+str(m_over_n), fontsize = label_size)
    
    ax.spines['top'].set_linewidth(2.5)
    ax.spines['left'].set_linewidth(2.5)
    ax.spines['right'].set_linewidth(2.5)
    ax.spines['bottom'].set_linewidth(2.5) 
    ax.tick_params(axis='both', width=2.5)   
    
    cbar.ax.spines['top'].set_linewidth(2.5)
    cbar.ax.spines['left'].set_linewidth(2.5)
    cbar.ax.spines['right'].set_linewidth(2.5)
    cbar.ax.spines['bottom'].set_linewidth(2.5) 
    cbar.ax.tick_params(axis='both', width=2.5)       
    
    plt.tight_layout()

    #=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    # BASIC CHI-PLOT WITH EACH CHANNEL LABELLED WITH A DIFFERENT COLOUR
    plt.figure(1, facecolor='white',figsize=(10,7.5))
    ax = plt.subplot(1,1,1)

    
    data_pointer = 0    # points to data element in chi and elev vectors
    for i in range (0,len(segment_lengths)):
        chi_seg = np.zeros(segment_lengths[i])
        elev_seg = np.zeros(segment_lengths[i])

def map_dict(export_fn,states_happiness):
    # Lambert Conformal map of lower 48 states.
    m = Basemap(llcrnrlon=-119,llcrnrlat=22,
                    urcrnrlon=-64, urcrnrlat=49,
                                projection='lcc', lat_1=33,lat_2=45,lon_0=-95)
    
    # laod state boundaries.
    # data from U.S Census Bureau
    # http://www.census.gov/geo/www/cob/st1990.html
    shp_info = m.readshapefile('gz_2010_us_040_00_500k/gz_2010_us_040_00_500k','states',drawbounds=False)
    # This loads three files:
    #   gz_2010_us_040_00_500k.dbf
    #   gz_2010_us_040_00_500k.shp
    #   gz_2010_us_040_00_500k.shx
    
    max_score = -sys.maxint - 1
    min_score = sys.maxint+1
    happiest_state = ''
    saddest_state = ''
    for state in states_happiness:
          score = states_happiness[state]
          if (score > max_score):
                happiest_state = state
                max_score = score
          if (score < min_score):
                saddest_state = state
                min_score = score
    
    for state in states_happiness:
        states_happiness[state] -= min_score
    
    min_score = 0
    max_score = max_score-min_score
    
    # choose a color for each state based on happiness score.
    colors={}
    statenames=[]
    cmap = plt.cm.Blues_r # use 'hot' colormap
    vmin = min_score; vmax = max_score # set range.
    sm = plt.cm.ScalarMappable(cmap=cmap, 
                                   norm=plt.normalize(vmin=vmin, vmax=vmax))
    
    for shapedict in m.states_info:
        statename = shapedict['NAME']
        #print statesNames[statename], ' ', statename
        try:
            #score = states_happiness[statesNames[statename]]
            score = states_happiness[statename]
        except KeyError:
            score = 0.0
            
        # calling colormap with value between 0 and 1 returns
        # rgba value. Invert color range (hot colors are
        # high
        # population), take sqrt root to spread out
        # colors more.
        try:
            colors[statename] = cmap((score-vmin)/(vmax-vmin))[:3]
        except KeyError:
            colors[statename] = cmap(0)[:3]
        statenames.append(statename)
        
    # cycle through state names,
    # color each one.
    for nshape,seg in enumerate(m.states):
        xx,yy = zip(*seg)
        if statenames[nshape] != 'District of Columbia' and \
        statenames[nshape] != "Puerto Rico":
            color = rgb2hex(colors[statenames[nshape]]) 
            plt.fill(xx,yy,color,edgecolor='black')
    
    # draw
    # meridians
    # and
    # parallels.
    m.drawparallels(np.arange(25,65,20),   labels=[0,0,0,0],
                        zorder=-1,color="w")
    m.drawmeridians(np.arange(-120,-40,20),labels=[0,0,0,0],
                        zorder=-1,color="w")
     
    # set  up colorbar:
    mm = plt.cm.ScalarMappable(cmap=cmap)
    mm.set_array([0,1])
    #plt.colorbar(mm, label="Happiness",
    #                              orientation="horizontal", fraction=0.05)
    plt.colorbar(mm,orientation="horizontal", fraction=0.05)
     
    #plt.title('Filling State Polygons by Population Density')
    plt.gca().axis("off")
    plt.savefig(export_fn)
    plt.close('all')

#.........这里部分代码省略.........
  optParser.add_option("--log-widths", dest="logWidths", action="store_true",
                         default=False, help="If set, widths are log-scaled")
  optParser.add_option("--min-color-value", dest="colorMin", 
                         type="float", default=None, help="If set, defines the minimum edge color value")
  optParser.add_option("--max-color-value", dest="colorMax", 
                         type="float", default=None, help="If set, defines the maximum edge color value")
  optParser.add_option("--min-width-value", dest="widthMin", 
                         type="float", default=None, help="If set, defines the minimum edge width value")
  optParser.add_option("--max-width-value", dest="widthMax", 
                         type="float", default=None, help="If set, defines the maximum edge width value")
  optParser.add_option("-v", "--verbose", dest="verbose", action="store_true",
                         default=False, help="If set, the script says what it's doing")
  # standard plot options
  helpers.addInteractionOptions(optParser)
  helpers.addPlotOptions(optParser)
  # parse
  options, remaining_args = optParser.parse_args(args=args)

  if options.net==None: 
    print "Error: a network to load must be given."
    return 1
  if options.verbose: print "Reading network from '%s'" % options.net
  net = sumolib.net.readNet(options.net)
  
  if options.measures==None: 
    print "Error: a dump file must be given."
    return 1
  
  times = []
  hc = None
  if options.dumps.split(",")[0]!="":  
    if options.verbose: print "Reading colors from '%s'" % options.dumps.split(",")[0]
    hc = WeightsReader(options.measures.split(",")[0])
    sumolib.output.parse_sax(options.dumps.split(",")[0], hc)
    times = hc._edge2value
  
  hw = None
  if options.dumps.split(",")[1]!="":  
    if options.verbose: print "Reading widths from '%s'" % options.dumps.split(",")[1]
    hw = WeightsReader(options.measures.split(",")[1])
    sumolib.output.parse_sax(options.dumps.split(",")[1], hw)
    times = hw._edge2value

  for t in times: 
    colors = {}
    maxColorValue = None
    minColorValue = None
    for e in net._id2edge:
      if hc and t in hc._edge2value and e in hc._edge2value[t]: 
        if options.colorMax!=None and hc._edge2value[t][e]>options.colorMax: hc._edge2value[t][e] = options.colorMax
        if options.colorMin!=None and hc._edge2value[t][e]<options.colorMin: hc._edge2value[t][e] = options.colorMin
        if maxColorValue==None or maxColorValue<hc._edge2value[t][e]: maxColorValue = hc._edge2value[t][e] 
        if minColorValue==None or minColorValue>hc._edge2value[t][e]: minColorValue = hc._edge2value[t][e]
        colors[e] = hc._edge2value[t][e] 
    if options.colorMax!=None: maxColorValue = options.colorMax 
    if options.colorMin!=None: minColorValue = options.colorMin 
    if options.logColors: 
      helpers.logNormalise(colors, maxColorValue)
    else:
      helpers.linNormalise(colors, minColorValue, maxColorValue)
    for e in colors:
      colors[e] = helpers.getColor(options, colors[e], 1.)
    if options.verbose: print "Color values are between %s and %s" % (minColorValue, maxColorValue)
  
    widths = {}
    maxWidthValue = None
    minWidthValue = None
    for e in net._id2edge:
      if hw and t in hw._edge2value and e in hw._edge2value[t]:
        v = abs(hw._edge2value[t][e]) 
        if options.widthMax!=None and v>options.widthMax: v = options.widthMax
        if options.widthMin!=None and v<options.widthMin: v = options.widthMin
        if not maxWidthValue or maxWidthValue<v: maxWidthValue = v 
        if not minWidthValue or minWidthValue>v: minWidthValue = v
        widths[e] = v
    if options.widthMax!=None: maxWidthValue = options.widthMax 
    if options.widthMin!=None: minWidthValue = options.widthMin 
    if options.logWidths: 
      helpers.logNormalise(widths, options.colorMax)
    else:
      helpers.linNormalise(widths, minWidthValue, maxWidthValue)
    for e in widths:
      widths[e] = options.minWidth + widths[e] * (options.maxWidth-options.minWidth)
    if options.verbose: print "Width values are between %s and %s" % (minWidthValue, maxWidthValue)
  
    fig, ax = helpers.openFigure(options)
    ax.set_aspect("equal", None, 'C')
    helpers.plotNet(net, colors, widths, options)
    
    
    # drawing the legend, at least for the colors
    sm = matplotlib.cm.ScalarMappable(cmap=get_cmap(options.colormap), norm=plt.normalize(vmin=minColorValue, vmax=maxColorValue))
    # "fake up the array of the scalar mappable. Urgh..." (pelson, http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots)
    sm._A = []
    plt.colorbar(sm)
    
    options.nolegend = True
    helpers.closeFigure(fig, ax, options)
    return 0 # !!! hack 
  return 0

    f = liblas.file.File(sys.argv[1],mode='r')

    lgrid = LidarGrid (f)

    print lgrid.minz
    print lgrid.maxz
    
    print "Adding data to grid"
    for p in f:
        if p.return_number == 1:
            lgrid.add (p.x, p.y, p.z)
    #   print "Add: "+str(p.x)+", "+str(p.y)+", "+str(p.z)

    print "Data added to grid."

    

    plt.subplot(111)
    n = plt.normalize(75, 155)
#img = plt.imshow(lgrid.grid, extent=(0, lgrid.xsize, lgrid.ysize, 0), cmap=plt.cm.Greys_r)
    img = plt.imshow(lgrid.grid, norm=n, interpolation='nearest', cmap=plt.cm.jet)
#    img.set_clim(lgrid.minz, 120)
    plt.colorbar()
    plt.show()

    print "Done."

#for x in range (0, lidar_grid.xsize):
#   for y in range (0, lidar_grid.ysize):
#       print "Get: "+str(lidar_grid.get (x, y))

#.........这里部分代码省略.........
    col = []
        
    with open(tree_file, 'r') as f:
        lines = f.readlines()
        
    for q,line in enumerate(lines):
        if q > 0: #skip first line
            channel_id.append(float(line.split()[0]))
            M_chi_value.append(float(line.split()[11]))
            row.append(float(line.split()[4]))
            col.append(float(line.split()[5]))

    #get bounding box & pad to 10% of the dimension
    x_max = max(col)
    x_min = min(col)
    y_max = max(row)
    y_min = min(row) 
    
    pad_x = (x_max - x_min) * 0.1
    pad_y = (x_max - y_min) * 0.1
    
    if (pad_y > pad_x):
        pad_x = pad_y
    else:
        pad_y = pad_x
    
    x_max += pad_x
    x_min -= pad_x
    y_max += pad_y
    y_min -= pad_y 
    
    fig = pp.figure(1, facecolor='white',figsize=(10,7.5))
    ax = fig.add_subplot(1,1,1)
    ax.imshow(hillshade, vmin=0, vmax=255, cmap=cmx.gray)
    
    # now get the tick marks    
    n_target_tics = 5
    xlocs,ylocs,new_x_labels,new_y_labels = format_ticks_for_UTM_imshow(hillshade_header,x_max,x_min,y_max,y_min,n_target_tics)  
    pp.xticks(xlocs, new_x_labels, rotation=60)  #[1:-1] skips ticks where we have no data
    pp.yticks(ylocs, new_y_labels) 
    
    for line in ax.get_xticklines():
        line.set_marker(mpllines.TICKDOWN)
        #line.set_markeredgewidth(3)

    for line in ax.get_yticklines():
        line.set_marker(mpllines.TICKLEFT)
        #line.set_markeredgewidth(3)  
 
    pp.xlim(x_min,x_max)    
    pp.ylim(y_max,y_min)  
   
    pp.xlabel('Easting (m)',fontsize = axis_size)
    pp.ylabel('Northing (m)',fontsize = axis_size)    
              
    # channel ID
    M_chi_value_MIN = np.min(M_chi_value)
    M_chi_value_MAX = np.max(M_chi_value)
    cNorm_M_chi_value  = colors.Normalize(vmin=M_chi_value_MIN, vmax=M_chi_value_MAX)  # the max number of channel segs is the 'top' colour
    hot = pp.get_cmap('RdYlBu_r')
    scalarMap_M_chi_value = cmx.ScalarMappable(norm=cNorm_M_chi_value, cmap=hot) 
    
    
    for a,i in enumerate(M_chi_value):
        #print "a: " +str(a)+" i: " +str(i)
        if channel_id[a] != 0:     
            # plot other stream segments
            colorVal = scalarMap_M_chi_value.to_rgba(i) # this gets the distinct colour for this segment
            pp.scatter(col[a], row[a], 30,marker=".", color=colorVal,edgecolors=colorVal) 

    for a,i in enumerate(M_chi_value):
        if channel_id[a] == 0:
            # plot trunk stream in black
            colorVal = scalarMap_M_chi_value.to_rgba(i)
            pp.scatter(col[a], row[a], 40,marker=".", color=colorVal,edgecolors=colorVal)

    sm = pp.cm.ScalarMappable(cmap=hot, norm=pp.normalize(vmin=min(M_chi_value), vmax=max(M_chi_value)))
    sm._A = []
    


    
    ax.spines['top'].set_linewidth(2.5)
    ax.spines['left'].set_linewidth(2.5)
    ax.spines['right'].set_linewidth(2.5)
    ax.spines['bottom'].set_linewidth(2.5) 
    ax.tick_params(axis='both', width=2.5)      
    
    from mpl_toolkits.axes_grid1 import make_axes_locatable
    divider = make_axes_locatable(pp.gca())
    cax = divider.append_axes("right", "5%", pad="3%")
    pp.colorbar(sm, cax=cax).set_label('$M_{\chi}$',fontsize=axis_size) 
    cax.tick_params(labelsize=label_size) 
    
    #pp.xlim(x_min,x_max)    
    #pp.ylim(y_max,y_min) 

    pp.tight_layout()
        
    pp.show()

import numpy as np
from scipy.interpolate import Rbf
import matplotlib.pyplot as plt
from matplotlib import cm

# 2-d tests - setup scattered data
x = np.random.rand(100)*4.0-2.0
y = np.random.rand(100)*4.0-2.0
z = x*np.exp(-x**2-y**2)
ti = np.linspace(-2.0, 2.0, 100)
XI, YI = np.meshgrid(ti, ti)

# use RBF
rbf = Rbf(x, y, z, epsilon=2)
ZI = rbf(XI, YI)

# plot the result
n = plt.normalize(-2., 2.)
plt.subplot(1, 1, 1)
plt.pcolor(XI, YI, ZI, cmap=cm.jet)
plt.scatter(x, y, 100, z, cmap=cm.jet)
plt.title('RBF interpolation - multiquadrics')
plt.xlim(-2, 2)
plt.ylim(-2, 2)
plt.colorbar()
plt.show()

def main(args=None):
  """The main function; parses options and plots"""
  ## ---------- build and read options ----------
  from optparse import OptionParser
  optParser = OptionParser()
  optParser.add_option("-n", "--net", dest="net", metavar="FILE",
                         help="Defines the network to read")
  optParser.add_option("--edge-width", dest="defaultWidth", 
                         type="float", default=1, help="Defines the edge width")
  optParser.add_option("--edge-color", dest="defaultColor", 
                         default='k', help="Defines the edge color")
  optParser.add_option("--minV", dest="minV", 
                         type="float", default=None, help="Define the minimum value boundary")
  optParser.add_option("--maxV", dest="maxV", 
                         type="float", default=None, help="Define the maximum value boundary")
  optParser.add_option("-v", "--verbose", dest="verbose", action="store_true",
                         default=False, help="If set, the script says what it's doing")
  # standard plot options
  helpers.addInteractionOptions(optParser)
  helpers.addPlotOptions(optParser)
  # parse
  options, remaining_args = optParser.parse_args(args=args)

  if options.net==None: 
    print "Error: a network to load must be given."
    return 1
  if options.verbose: print "Reading network from '%s'" % options.net
  net = sumolib.net.readNet(options.net)

  speeds = {}
  minV = None
  maxV = None
  for e in net._id2edge:
    v = net._id2edge[e]._speed
    if minV==None or minV>v:
        minV = v