def test_jpeg_alpha():
    Image = pytest.importorskip('PIL.Image')

    plt.figure(figsize=(1, 1), dpi=300)
    # Create an image that is all black, with a gradient from 0-1 in
    # the alpha channel from left to right.
    im = np.zeros((300, 300, 4), dtype=float)
    im[..., 3] = np.linspace(0.0, 1.0, 300)

    plt.figimage(im)

    buff = io.BytesIO()
    with rc_context({'savefig.facecolor': 'red'}):
        plt.savefig(buff, transparent=True, format='jpg', dpi=300)

    buff.seek(0)
    image = Image.open(buff)

    # If this fails, there will be only one color (all black). If this
    # is working, we should have all 256 shades of grey represented.
    num_colors = len(image.getcolors(256))
    assert 175 <= num_colors <= 185
    # The fully transparent part should be red.
    corner_pixel = image.getpixel((0, 0))
    assert corner_pixel == (254, 0, 0)

def plotClusters(websites, amount, clusters=8,xFactor=10, yFactor=10, myDpi=96, sampleSize=20):
    """
     We want to plot every image according to the appropriate point
     on the x-axis according to its cluster number. We want to plot
     each new member of a given cluster at a higher y position

     Inputs:
         imagePath is a string to where the images are located
         websites is a list of tuples of the form:
             (clusterNumber, websitename)
    """
    imagePath = getDataDir(amount, cut=True, big=False)
    clusterDict = [0 for n in xrange(clusters)]
    plt.figure(figsize=(800/myDpi, 800/myDpi), dpi=myDpi)
    for site in websites:
        clusterIndex, address = site
        try:
            yIndex = clusterDict[clusterIndex]
            if yIndex > sampleSize:
                continue
            image = mpimg.imread(imagePath+address)
            y = yIndex * yFactor
            clusterDict[clusterIndex]+=1
            plt.figimage(image, clusterIndex*xFactor, y)
        except IOError:
            # usually if we don't have the small cut image yet
            pass

    plt.show()

def plot_montage(images, ndx, ncol=None):

    N = len(ndx)
    
    if not ncol:
        ncol = int(sqrt(N))

    f = plt.figure(dpi=100)
    
    row = 0
    col = 0
    tot_height = 0
    for i in range(N):

        I = sp.array(images[ndx[i]], dtype='float')
        I = I / I.max()
        
        height = I.shape[0]
        width = I.shape[1]
        ax = plt.figimage(I, xo=width*col, yo=height*row, origin='lower')

        col += 1
        if col % ncol == 0:
            row += 1
            col = 0
            tot_height += height

    tot_height += height
    tot_width = width*ncol
    
    f.set_figheight(tot_height/100)
    f.set_figwidth(tot_width/100)

    return f

def plot_map(res, bins, cmap, cropName, fileName, technologyName, variableName, unitLabel, techFolder):
	
	figTitle = cropName + ' ' + technologyName + ' ' + variableName + ' (' + unitLabel + ')' 

	plt.close()
	plt.figtext(0.025,0.92, figTitle, clip_on = 'True', size = 'large', weight = 'semibold'); 
	plt.figtext(0.025,0.86, 'Spatially disaggregated production statistics of circa 2005 using the Spatial Production Allocation Model (SPAM).', clip_on = 'True', size = 'x-small', stretch = 'semi-condensed', weight = 'medium');
	plt.figtext(0.025,0.835,'Values are for 5 arc-minute grid cells.', clip_on = 'True', size = 'x-small', stretch = 'semi-condensed', weight = 'medium');

	plt.figtext(0.025,0.072, 'You, L., U. Wood-Sichra, S. Fritz, Z. Guo, L. See, and J. Koo. 2014', clip_on = 'True', size = 'xx-small', stretch = 'semi-condensed', weight = 'medium')
	plt.figtext(0.025,0.051, 'Spatial Production Allocation Model (SPAM) 2005 Version 2.0.', clip_on = 'True', size = 'xx-small', stretch = 'semi-condensed', weight = 'medium')
	plt.figtext(0.025,0.030, '03.10.2015. Available from http://mapspam.info', clip_on = 'True', size = 'xx-small', stretch = 'semi-condensed', weight = 'medium');

	plt.figimage(logos, 4100, 200)

	map = Basemap(projection='merc',resolution='i', epsg=4326, lat_0 = 0, lon_0 = 20, llcrnrlon=-160, llcrnrlat=-70, urcrnrlon=200, urcrnrlat=90)
	map.drawlsmask(land_color='#fffff0', lakes = False, zorder = 1)
	shp_coast = map.readshapefile(baseShpLoc + '/ne_50m_coastline/ne_50m_coastline', 'scalerank', drawbounds=True, linewidth=0.1, color = '#828282', zorder = 7)
	shp_rivers = map.readshapefile(baseShpLoc + '/ne_110m_rivers_lake_centerlines/ne_110m_rivers_lake_centerlines', 'scalerank', drawbounds=True, color='#e8f8ff', linewidth=0.1, zorder = 5)
	shp_lakes = map.readshapefile(baseShpLoc + '/ne_110m_lakes/ne_110m_lakes', 'scalerank', drawbounds=True, linewidth=0.1, color='#e8f8ff', zorder=4)

	paths = []
	for line in shp_lakes[4]._paths:
		paths.append(matplotlib.path.Path(line.vertices, codes=line.codes))
	coll_lakes = matplotlib.collections.PathCollection(paths, linewidths=0, facecolors='#e8f8ff', zorder=3)

	cs = map.pcolormesh(res.lons, res.lats, res.d2, cmap=cmap, norm=BoundaryNorm(bins, 256, clip=True), zorder = 6)

	map.drawcountries(linewidth=0.1, color='#828282', zorder = 8)

	ax = plt.gca()
	ax.add_collection(coll_lakes)

	cbar = map.colorbar(cs,location='bottom', pad='3%')

	if (variableName == 'Production' or variableName == 'Yield'):
		labelSize = 8
	else :
		labelSize = 7
	labels = [item.get_text() for item in cbar.ax.get_xticklabels()]
	labels[0] = '1'; labels[labelSize - 1] = labels[labelSize - 1] + ' <'; labels[labelSize] = ''
	cbar.ax.set_xticklabels(labels)

	plt.tight_layout(h_pad=0.9, w_pad = 0.9)
	
	outputFile = 'spam2005v2r0_' + techFolder + '_' + fileName + '_' + technologyName.lower()
	plt.savefig(outputFolder + '/' + techFolder + '/' + outputFile + '.png', format='png', dpi=900)

def plot_digits(X, ndx, ncol=50, width=IM_WIDTH, cmap=cm.gray):
    """
    plot a montage of the examples specified in ndx (as rows of X)
    """

    row = 0
    col = 0
    for i in range(ndx.shape[0]):
        
        plt.figimage(reshape_digit(X, ndx[i]),
                     xo=width*col, yo=width*row,
                     origin='upper', cmap=cmap)

        col += 1
        if col % ncol == 0:
            row += 1
            col = 0

def draw_mol_slider(mol, legend, matching, symbol, color, values, index):
    import matplotlib.pyplot as plt
    from matplotlib.widgets import Slider
    if op.isfile('/home/flo/temp.png'):
        os.remove('/home/flo/temp.png')
    subimg = MolToImage(mol, (200, 200), legend=legend, highlightAtoms=matching, symbol=symbol, background=color,
                        mol_adder=MolDrawing.AddMol)
    # Then we will integrate a slider showing the relative ranking of the molecule
    rank = float(len(values[:index + 1])) / len(values)
    alpha_axis = plt.axes([0.37, 0.59, 0.14, 0.018], axisbg='grey')
    Slider(alpha_axis, '1(+)', 0, 1, valinit=rank, facecolor='w', valfmt='(-)%.0f')
    # Manually set to look good on the svg: [0.37, 0.59, 0.14, 0.018] for alpha_axis and 180, 70 for subimg
    # Manually set to look good on the png: [0.37, 0.59, 0.14, 0.018] for alpha_axis and 247, 172 for subimg
    # Now we need to combine subimg and slider in one image
    plt.figimage(subimg, 247, 172)
    plt.show()
    plt.savefig('/home/flo/temp.png')  # , bbox_inches='tight')

def _mock_worker(exact, estimated, param, nologo):
    """Plot the exact and estimated values of a parameter for the mock analysis

    Parameters
    ----------
    exact: Table column
        Exact values of the parameter.
    estimated: Table column
        Estimated values of the parameter.
    param: string
        Name of the parameter
    nologo: boolean
        Do not add the logo when set to true.

    """

    range_exact = np.linspace(np.min(exact), np.max(exact), 100)

    # We compute the linear regression
    if (np.min(exact) < np.max(exact)):
        slope, intercept, r_value, p_value, std_err = stats.linregress(exact,
                                                                       estimated)
    else:
        slope = 0.0
        intercept = 1.0
        r_value = 0.0

    plt.errorbar(exact, estimated, marker='.', label=param, color='k',
                 linestyle='None', capsize=0.)
    plt.plot(range_exact, range_exact, color='r', label='1-to-1')
    plt.plot(range_exact, slope * range_exact + intercept, color='b',
             label='exact-fit $r^2$ = {:.2f}'.format(r_value**2))
    plt.xlabel('Exact')
    plt.ylabel('Estimated')
    plt.title(param)
    plt.legend(loc='best', fancybox=True, framealpha=0.5, numpoints=1)
    plt.minorticks_on()
    if nologo is False:
        image = plt.imread(pkg_resources.resource_filename(__name__,
                                                           "data/CIGALE.png"))
        plt.figimage(image, 510, 55, origin='upper', zorder=10, alpha=1)

    plt.tight_layout()
    plt.savefig(OUT_DIR + 'mock_{}.pdf'.format(param))

    plt.close()

def hydro_plot(view, hydro_code, image_size, figname):
    """
    view: the (physical) region to plot [xmin, xmax, ymin, ymax]
    hydro_code: hydrodynamics code in which the gas to be plotted is defined
    image_size: size of the output image in pixels (x, y)
    """
    if not HAS_MATPLOTLIB:
        return
    shape = (image_size[0], image_size[1], 1)
    size = image_size[0] * image_size[1]
    axis_lengths = [0.0, 0.0, 0.0] | units.m
    axis_lengths[0] = view[1] - view[0]
    axis_lengths[1] = view[3] - view[2]
    grid = Grid.create(shape, axis_lengths)
    grid.x += view[0]
    grid.y += view[2]
    speed = grid.z.reshape(size) * (0 | 1/units.s)
    rho, rhovx, rhovy, rhovz, rhoe = hydro_code.get_hydro_state_at_point(
            grid.x.reshape(size),
            grid.y.reshape(size),
            grid.z.reshape(size), speed, speed, speed)

    min_v = 800.0 | units.km / units.s
    max_v = 3000.0 | units.km / units.s
    min_rho = 3.0e-9 | units.g / units.cm**3
    max_rho = 1.0e-5 | units.g / units.cm**3
    min_E = 1.0e11 | units.J / units.kg
    max_E = 1.0e13 | units.J / units.kg

    v_sqr = (rhovx**2 + rhovy**2 + rhovz**2) / rho**2
    E = rhoe / rho
    log_v = numpy.log((v_sqr / min_v**2)) / numpy.log((max_v**2 / min_v**2))
    log_rho = numpy.log((rho / min_rho)) / numpy.log((max_rho / min_rho))
    log_E = numpy.log((E / min_E)) / numpy.log((max_E / min_E))

    red = numpy.minimum(numpy.ones_like(rho.number), numpy.maximum(
        numpy.zeros_like(rho.number), log_rho)).reshape(shape)
    green = numpy.minimum(numpy.ones_like(rho.number), numpy.maximum(
        numpy.zeros_like(rho.number), log_v)).reshape(shape)
    blue = numpy.minimum(numpy.ones_like(rho.number), numpy.maximum(
        numpy.zeros_like(rho.number), log_E)).reshape(shape)
    alpha = numpy.minimum(
            numpy.ones_like(log_v),
            numpy.maximum(
                numpy.zeros_like(log_v),
                numpy.log((rho / (10*min_rho)))
                )
            ).reshape(shape)

    rgba = numpy.concatenate((red, green, blue, alpha), axis=2)

    pyplot.figure(figsize=(image_size[0]/100.0, image_size[1]/100.0), dpi=100)
    im = pyplot.figimage(rgba, origin='lower')

    pyplot.savefig(figname, transparent=True, dpi=100,
                   facecolor='k', edgecolor='k')
    print "\nHydroplot was saved to: ", figname
    pyplot.close()

	def make_image(self):
		rawdata = load(self.binfile)
		rotated = rot90(rawdata)
		rotated = rot90(rotated)
		rotated = rot90(rotated)
		rotated = reshape(rotated,(1,32,32))

		rawshape = shape(rawdata)
		
		if self.taking_sky == True:
			self.skycount += rotated
		
		#print 'rawdata shape ', shape(rawdata)
		newcounts = reshape(rawdata,(1,1024))
		self.counts = append(self.counts,newcounts,axis=0)
		
		self.rotated_counts = append(self.rotated_counts, rotated, axis=0)
		
		#print 'shape of counts', shape(self.counts)
		tf = self.image_time
		self.int_time = self.ui.int_time_spinBox.value()
		ti = tf-self.int_time
		if ti<0:
			ti = 0
		if tf==0 or tf==ti:
			image_counts = newcounts
		else:
			image_counts = sum(self.counts[ti:tf],axis=0)
			image_counts = reshape(image_counts,(1,1024))

		if self.sky_subtraction == True:
			image_counts = image_counts-reshape(self.skyrate*(tf-ti),(1,1024))

		indices = sort(image_counts)
		
		brightest = self.ui.brightpix.value()
		self.vmax = indices[0,-1*(brightest)]

		photon_count = reshape(image_counts,rawshape)
		photon_count = rot90(photon_count)
		photon_count = rot90(photon_count)
		photon_count = rot90(photon_count)
		
		fig = plt.figure(figsize=(.32,.32), dpi=100, frameon=False)
		im = plt.figimage(photon_count, cmap='gray', vmax = self.vmax)
		#im = plt.figimage(rawdata, cmap='gray')
		plt.savefig("Arcons_frame.png", pad_inches=0)
		print "Generated image ",tf
		self.image_time+=1
		if self.taking_sky == True:
			if self.image_time == self.skytime:
				self.taking_sky = False
				self.skyrate = self.skycount/self.skytime

def compare_files(dir_path1, dir_path2, name):
    result = None
    
    path1 = os.path.join(dir_path1, name)
    path2 = os.path.join(dir_path2, name)
    
    if not os.path.isfile(path1) or not os.path.isfile(path2):
        result = Difference(name, 'missing file')
    else:
        image1 = plt.imread(path1)
        image2 = plt.imread(path2)
        if image1.shape != image2.shape:
            result = Difference(name, "shape: %s -> %s" % (image1.shape, image2.shape))
        else:
            diff = (image1 != image2).any(axis=2)
            if diff.any():
                diff_path = os.path.join(DIFF_DIR, name)
                diff_path = os.path.realpath(diff_path)
                plt.figure(figsize=reversed(diff.shape), dpi=1)
                plt.figimage(diff, cmap=cm.gray)
                plt.savefig(diff_path, dpi=1)
                result = Difference(name, "diff: %s" % diff_path)
    return result

def testDivisonOfGrid():
	topo = topology.Topology()
	topoFile="topology_200.txt"
	try:
		topo.importFromFile(topoFile)
	except: 
		print "could not read topology file",topoFile
	
	habNM = simulation.makeNeighConnMap(topo)
	sects = simulation.divideGridToSection(habNM,GRID_SIZE,GRID_SIZE)
	print "sects ",sects
	k=0
	for habs in sects.values():
		k=k+1
		X = pylab.ones((10,10,3))*float(k*0.4)
		print "habs" ,habs
		for hab in habs:
			hab=hab-1
			xo=int(hab%GRID_SIZE)*2
			yo=int(hab/GRID_SIZE)*2
			print "X ",xo
			print "Y ",yo
			plt.figimage(X, xo, yo, origin='lower')
	plt.show()

def test_jpeg_alpha():
    plt.figure(figsize=(1, 1), dpi=300)
    # Create an image that is all black, with a gradient from 0-1 in
    # the alpha channel from left to right.
    im = np.zeros((300, 300, 4), dtype=float)
    im[..., 3] = np.linspace(0.0, 1.0, 300)

    plt.figimage(im)

    buff = io.BytesIO()
    with rc_context({'savefig.facecolor': 'red'}):
        plt.savefig(buff, transparent=True, format='jpg', dpi=300)

    buff.seek(0)
    image = Image.open(buff)

    # If this fails, there will be only one color (all black). If this
    # is working, we should have all 256 shades of grey represented.
    print("num colors: ", len(image.getcolors(256)))
    assert len(image.getcolors(256)) >= 175 and len(image.getcolors(256)) <= 185
    # The fully transparent part should be red, not white or black
    # or anything else
    print("corner pixel: ", image.getpixel((0, 0)))
    assert image.getpixel((0, 0)) == (254, 0, 0)

def saveFig(w,h,data,output,cmapData=None, minData=None, maxData=None):

    # --- Save image
    fig = plt.figure(figsize=((w/100)+1, (h/100)+1), dpi=100)
    cax = plt.figimage(data, vmin=minData, vmax=maxData, cmap=cmapData)
    fig.savefig(output)#, bbox_inches=extent)
    plt.close()

    # Load and resave image
    im = Image.open(output)
    (widthN, heightN) = im.size
    logger.info("Detected size: ",widthN,heightN, "targeted", w, h)
    im2 = im.crop(((widthN-w),
                   (heightN-h),
                   w,h))
    im.close()
    im2.save(output)

    def loadIm(self):
        logger.debug("Complex load")
        fig = plt.figure(figsize=((self.width/100.0), (self.height/100.0)), dpi=100)
        data = convertImage(self.pixelsHDR, self.height,
                            self.width, self.inverse)

        if(self.inverse):
            pRefLum = [0.0 if lum(p)==0 else 1.0/lum(p) for p in self.pixelsHDR]
        else:
            pRefLum = [lum(p) for p in self.pixelsHDR]
        pRefLum.sort()

        maxData = pRefLum[-1]
        minData = pRefLum[0]

        logger.info("Find the min/max: ",str(minData),str(maxData))


        if(self.pMax != -1):
            maxData = pRefLum[int(len(pRefLum)*self.pMax)]
        if(self.pMin != -1):
            minData = pRefLum[int(len(pRefLum)*self.pMin)]

        if self.minV != -10000.005454:
            minData = self.minV
        if self.maxV != 10000.005454:
            maxData = self.maxV

        logger.info("Used min/max: ",str(minData),str(maxData))

        # --- Save the figure
        cax = plt.figimage(data, vmin=minData, vmax=maxData, cmap=self.cmap)
        fig.savefig(wk + os.path.sep + self.output)#, bbox_inches=extent)

        self.im = Image.open(wk + os.path.sep + self.output)
        (widthN, heightN) = self.im.size
        self.im = self.im.crop(((widthN-self.width),
                                (heightN-self.height),
                                self.width,
                                self.height))

def make_png(infile, outdir):
    """
    """
    # Get DN
    sarmp = SARMapper(infile)
    sarim = SARImage(sarmp)
    orig_spacing = sarim.get_original_spacing()
    # n = np.ceil(sarim.get_info('number_of_samples')/600.)
    # spacing = orig_spacing[1]*n
    # spacing = orig_spacing[1]*10
    im = abs(sarim.get_values('digital_number', step=[10, 10]))
    final_spacing = orig_spacing*[10, 10]
    sar_size = sarim.get_full_extent()[2:4]
    extent = [sar_size[0]/2, sar_size[1]/2, sar_size[0]/2+1, sar_size[1]/2+1]
    lon = sarim.get_values('lon', extent=extent)
    lat = sarim.get_values('lat', extent=extent)
    inc = sarim.get_values('incidence', extent=extent)
    heading = sarim.get_info('platform_heading')
    north = 90+heading # north dir in image
    im_num = sarim.get_info('image_number')
    dist = (np.array(im.shape)-1)*final_spacing/1000.
    if sarim.get_info('pass') == 'Descending':
        im = im[::-1, ::-1]
        north = north+180
    # ind = im.nonzero()
    # ind = np.where(im > 50)
    # im2 = im[ind]
    # vmin = im2.min()
    # vmax = im2.max()
    imsh = im.shape
    im2 = im[imsh[0]*0.1:imsh[0]*0.9, imsh[1]*0.1:imsh[1]*0.9]
    vmin = im2.min()
    vmax = im2.max()
    # Make figure
    dpi = 100
    imsize = np.array(im.shape[::-1])
    margin = np.array(((900-imsize[0])/2, 60))
    figsize = np.array((900, imsize[1]+2*margin[1]))
    imsizeN = imsize.astype('float')/figsize
    marginN = margin.astype('float')/figsize
    #print imsize, margin, figsize
    fig = plt.figure(figsize=figsize.astype('float')/dpi, dpi=dpi)
    # imax = fig.gca()
    # imax.set_position([marginN[0], marginN[1], imsizeN[0], imsizeN[1]])
    imax = fig.add_axes([marginN[0], marginN[1], imsizeN[0], imsizeN[1]])
    imax.set_xlim(0, dist[1])
    imax.set_ylim(0, dist[0])
    plt.imshow(im, origin='lower', interpolation='nearest', vmin=vmin,
                vmax=vmax, cmap=cm.get_cmap('Greys_r'),
                extent=[0, dist[1], 0, dist[0]], aspect='auto')
    imax.set_xlabel('range distance [km]')
    imax.set_ylabel('azimuth distance [km]')
    tit = '#%03i / lon=%.2f / lat=%.2f / inc=%.2f' % (im_num, lon, lat, inc)
    imax.set_title(tit)
    #imax.set_frame_on(False)
    #imax.set_axis_off()
    # Add colorbar
    cbax = fig.add_axes([1-0.75*marginN[0], .25, 20./figsize[0], .50])
    plt.colorbar(label='digital number', cax=cbax, format='%.1e')
    meanstr = r'$\mu$=%.2f' % im2.mean()
    cbax.text(0.5, -0.025, meanstr, ha='center', va='top')
    stdstr = r'$\sigma$=%.2f' % im2.std()
    cbax.text(0.5, -0.1, stdstr, ha='center', va='top')
    minstr = 'min=%.2f' % im2.min()
    cbax.text(0.5, 1.025, minstr, ha='center', va='bottom')
    maxstr = 'max=%.2f' % im2.max()
    cbax.text(0.5, 1.1, maxstr, ha='center', va='bottom')
    # Add north
    cpsizeN = (margin[0]-margin[1])/figsize.astype('float')
    cpax = fig.add_axes([0., 0.5-cpsizeN[1]/2, cpsizeN[0], cpsizeN[1]])
    plt.arrow(.5, .5, .5*np.cos(north*np.pi/180),
              .5*np.sin(north*np.pi/180), facecolor='black',
              width=0.01, length_includes_head=True,
              head_width=0.1, head_length=0.1)
    plt.annotate('North', (.5, .5), ha='center', va='top')
    cpax.set_axis_off()
    # Add Logos
    python_logo = '/local/home/gilles/Documents/logo/python-powered-w-70x28.png'
    #python_logo = '/home/cercache/users/gguitton/Documents/logo/python-powered-w-70x28.png'
    logo = imread(python_logo)
    plt.figimage(logo, 5, figsize[1]-logo.shape[0]-5)
    odl_logo = '/local/home/gilles/Documents/logo/oceandatalab-85x32.png'
    #odl_logo = '/home/cercache/users/gguitton/Documents/logo/oceandatalab-85x32.png'
    logo = imread(odl_logo)
    plt.figimage(logo, 5, 5)
    # Save as PNG
    infileid = os.path.basename(os.path.dirname(os.path.dirname(infile)))
    outdir2 = os.path.join(outdir, infileid)
    if os.path.exists(outdir2) == False:
        os.makedirs(outdir2)
    outfile = os.path.join(outdir2, os.path.basename(infile).replace('.tiff','-digital_number.png'))
    plt.savefig(outfile, dpi=dpi)#, bbox_inches='tight')
    plt.close()
    cmd = 'convert '+outfile+' -colors 256 '+outfile
    os.system(cmd)

"""
This illustrates placing images directly in the figure, with no axes.

"""
import numpy as np
import matplotlib
import matplotlib.cm as cm
import matplotlib.pyplot as plt


fig = plt.figure()
Z = np.arange(10000.0)
Z.shape = 100, 100
Z[:, 50:] = 1.

im1 = plt.figimage(Z, xo=50, yo=0, origin='lower')
im2 = plt.figimage(Z, xo=100, yo=100, alpha=.8, origin='lower')

plt.show()

    p = reader.PelFile(location+"%04i.pel"%run)
    image = np.sum(p.make3d(),axis=2)
    return (p.header._asdict(),image)

if __name__=="__main__":
    plt.figure(figsize=(4,4))
    g = gaingen()
    gain = g.next()
    old = -1000
    d={gain:{}}
    x = np.arange(512)

    for i in range(1,75):
        (h,im) = extract(i)
        plt.clf()
        plt.figimage(im,vmin=0,vmax=200,cmap="spectral")
        plt.savefig(location+"%04i.png"%i,dpi=128)
        plt.clf()
        plt.fill_between(x,np.sum(im,axis=0))
        plt.savefig(location+"x%04i.png"%i,dpi=128)
        plt.clf()
        plt.fill_between(x,np.sum(im,axis=1))
        plt.savefig(location+"y%04i.png"%i,dpi=128)
        if h[gain] < old: #we've moved on to the next gain
            print (gain,h[gain],old)
            gain = g.next()
            if not d.has_key(gain):
                d[gain] = {}
        d[gain][h[gain]] = i
        old = h[gain]
    with open(location+"gains.html","w") as outfile:

		return np.concatenate((np.arange(L), np.ones(2*L+1, dtype=np.int64)*L,
			np.arange(L-1,-L,-1), np.ones(2*L+1, dtype=np.int64)*-L, np.arange(1-L,0)))

	count = 1
	arr = np.ones((layers*2-1,layers*2-1), dtype=np.uint32)
	for layer in range(1, layers):

		# Generate the Y indices for the layer and offset based on the center coordinates
		sequenceY = getsequence(layer) + layers - 1

		# Offset X indices from the generated Y
		sequenceX = np.roll(sequenceY, -2*layer)

		# Assign incremental values to the indices in order
		for x, y in zip(sequenceX, sequenceY):
			count += 1
			arr[x,y] = count

	return arr

if __name__ == '__main__':
	# Testing
	import matplotlib.pyplot as plt
	layers = 100
	a = spiral(layers)
	for x in range(layers*2-1):
		for y in range(layers*2-1):
			a[x,y] = isprime(a[x,y])
	plt.figimage(a*255)
	plt.show()

ax = plt.subplot(2,2,4)
ax.imshow(small_im, interpolation = 'nearest')
plt.subplots_adjust(left = 0.24, wspace = 0.2, hspace = 0.1, \
    bottom = 0.05, top = 0.86)

#Label the rows and columns of the table
fig.text(0.03, 0.645, 'Big Image\nScaled Down', ha = 'left')
fig.text(0.03, 0.225, 'Small Image\nBlown Up', ha = 'left')
fig.text(0.383, 0.90, "Interpolation = 'none'", ha = 'center')
fig.text(0.75, 0.90, "Interpolation = 'nearest'", ha = 'center')

#If you were going to run this example on your local machine, you
#would save the figure as a PNG, save the same figure as a PDF, and
#then compare them.  The following code would suffice.
txt = fig.text(0.452, 0.95, 'Saved as a PNG', fontsize = 18)
# plt.savefig('None_vs_nearest-png.png')
# txt.set_text('Saved as a PDF')
# plt.savefig('None_vs_nearest-pdf.pdf')

#Here, however, we need to display the PDF on a webpage, which means
#the PDF must be converted into an image.  For the purposes of this
#example, the 'Nearest_vs_none-pdf.pdf' has been pre-converted into
#'Nearest_vs_none-pdf.png' at 80 dpi.  We simply need to load and
#display it.
pdf_im_path = cbook.get_sample_data('None_vs_nearest-pdf.png')
pdf_im = plt.imread(pdf_im_path)
fig2 = plt.figure(figsize = [8.0, 7.5])
plt.figimage(pdf_im)

plt.show()

def make_png(infile, outdir, vmax=None):
    """
    """
    # Get Sea Surface Roughness
    sarim = sarimage(infile)
    spacing_ra = np.ceil(sarim.get_info('number_of_samples')/600.)
    spacing_ra_m = sarim.pixels2meters(spacing_ra)[1]
    spacing = np.round(sarim.meters2pixels(spacing_ra_m))
    spacing_m = sarim.pixels2meters(spacing)
    ssr = sarim.get_data('roughness', spacing=spacing)
    lon = sarim.get_data('lon', midrange=True, midazimuth=True)[0, 0]
    lat = sarim.get_data('lat', midrange=True, midazimuth=True)[0, 0]
    inc = sarim.get_data('incidence', midrange=True, midazimuth=True)[0, 0]
    heading = sarim.get_info('platform_heading')
    north = 90 + heading # north dir in image
    im_num = sarim.get_info('image_number')
    dist = (np.array(ssr.shape)-1)*spacing_m/1000.
    # vmin = ssr[ind].min()
    # vmax = ssr[ind].max()
    vmin = scoreatpercentile(ssr[25:-25, 25:-25], 0.1)
    if vmax is None:
        vmax = scoreatpercentile(ssr[25:-25, 25:-25], 99.9)
    if sarim.get_info('pass') == 'Descending':
        ssr = ssr[::-1, ::-1]
        north = north+180
    # Make figure
    dpi = 100
    imsize = np.array(ssr.shape[::-1])
    margin = np.array(((900-imsize[0])/2, 60))
    figsize = np.array((900, imsize[1]+2*margin[1]))
    imsizeN = imsize.astype('float')/figsize
    marginN = margin.astype('float')/figsize
    #print imsize, margin, figsize
    fig = plt.figure(figsize=figsize.astype('float')/dpi, dpi=dpi)
    # imax = fig.gca()
    # imax.set_position([marginN[0], marginN[1], imsizeN[0], imsizeN[1]])
    imax = fig.add_axes([marginN[0], marginN[1], imsizeN[0], imsizeN[1]])
    imax.set_xlim(0, dist[1])
    imax.set_ylim(0, dist[0])
    plt.imshow(ssr, origin='lower', interpolation='nearest', vmin=vmin,
                vmax=vmax, cmap=cm.get_cmap('Greys_r'),
                extent=[0, dist[1], 0, dist[0]], aspect='auto')
    imax.set_xlabel('range distance [km]')
    imax.set_ylabel('azimuth distance [km]')
    tit = '#%03i / lon=%.2f / lat=%.2f / inc=%.2f' % (im_num, lon, lat, inc)
    imax.set_title(tit)
    #imax.set_frame_on(False)
    #imax.set_axis_off()
    # Add colorbar
    cbax = fig.add_axes([1-0.75*marginN[0], .25, 20./figsize[0], .50])
    plt.colorbar(label='sea surface roughness', cax=cbax, format='%.1e')
    meanstr = r'$\mu$=%.1f' % ssr.mean()
    cbax.text(0.5, -0.025, meanstr, ha='center', va='top')
    # Add north
    cpsizeN = (margin[0]-margin[1])/figsize.astype('float')
    cpax = fig.add_axes([0., 0.5-cpsizeN[1]/2, cpsizeN[0], cpsizeN[1]])
    plt.arrow(.5, .5, .5*np.cos(north*np.pi/180),
              .5*np.sin(north*np.pi/180), facecolor='black',
              width=0.01, length_includes_head=True,
              head_width=0.1, head_length=0.1)
    plt.annotate('North', (.5, .5), ha='center', va='top')
    cpax.set_axis_off()
    # Add Logos
    python_logo = os.path.join(LOGO_PATH, 'python-powered-w-70x28.png')
    logo = imread(python_logo)
    plt.figimage(logo, 5, figsize[1]-logo.shape[0]-5)
    odl_logo = os.path.join(LOGO_PATH, 'oceandatalab-85x32.png')
    logo = imread(odl_logo)
    plt.figimage(logo, 5, 5)
    # Save as PNG
    outbase = os.path.splitext(os.path.basename(infile))[0]+'-roughness.png'
    outfile = os.path.join(outdir, outbase)
    plt.savefig(outfile, dpi=dpi)#, bbox_inches='tight')
    plt.close()
    cmd = 'convert '+outfile+' -colors 256 '+outfile
    os.system(cmd)