def transform(W):
    global N
    print W[20]
    t = 20


    MM = N
    N = 356

    #bakk = W[t]
    arr = np.zeros((N+1,N+1)).astype(np.float32)

    # interpolation in order to keep baking in 32 array elements (for stability)
    bakk = interp1d(range(MM+1), W[t], kind='cubic')

    for i in range(-N/2,N/2+1):
        for j in range(-N/2,N/2+1):
            i2 = i
            j2 = j
            r = np.sqrt(i2*i2+j2*j2).astype(np.float32)
            if(r < N and r >= 0):
                arr[N/2-i,N/2-j] = bakk((N-r)*MM/N)

    if(False):
        I2 = Image.frombuffer('L',(N+1,N+1), (arr).astype(np.uint8),'raw','L',0,1)
        imgplot = plt.imshow(I2)
        plt.colorbar()
        gx, gy = np.gradient(arr)
        pylab.quiver(gx,gy)
        pylab.show()

        plt.show()
    #print "Saving Image res.png..."
    #I2.save('res.png')
    return arr

def getvectorfield( fitsa, fitsb, graph=False):
    mysex = sexparser.sextractorresult(fitsa)
    mysex.loaddata()
    testsex = sexparser.sextractorresult(fitsb)
    testsex.loaddata()

    flann = pyflann.FLANN()
    dataset = buildvec( mysex.objects["X_IMAGE"], mysex.objects["Y_IMAGE"] )
    testset = buildvec( testsex.objects["X_IMAGE"], testsex.objects["Y_IMAGE"] )

    ids, dists = flann.nn(dataset,testset,1,algorithm="linear")

    
    f=numpy.vectorize( lambda x:  x > 0 and x < 100 )
    cond = numpy.where( f(dists) )

    if graph:
	pylab.hist(dists,bins=100,range=(0,10),histtype="bar")
	pylab.hist(dists[cond],bins=100,range=(0,10))
	pylab.show()
    
    diffvec = dataset[ids[cond]]-testset[cond]

    x,y=numpy.transpose(testset)
    q,u=numpy.transpose(diffvec)

    if graph:
	pylab.quiver(x,y,q,u)
	pylab.xlim(0,2048)
	pylab.ylim(0,2048)
	pylab.show()

    return q.mean(), u.mean(), q.std(), u.std()

def draw_MAP_residuals(objectsA, objectsB, P, scaled='no'):
    from pyBA.distortion import compute_displacements, compute_residual
    from numpy import array
    
    # Compute displacements between frames for tie objects
    xobs, yobs, vxobs, vyobs, sxobs, syobs = compute_displacements(objectsA, objectsB)

    # Compute residual
    dx, dy = compute_residual(objectsA, objectsB, P)

    # Draw residuals
    fig = figure(figsize=(16,16))
    ax = fig.add_subplot(111, aspect='equal')
    if scaled is 'yes':
        # Allow relative scaling of arrows
        quiver(xobs,yobs,dx,dy)
    else:
        # Show residuals in absolute size (often very tiny), with uncertainties

        # Also plot error ellipses
        ellipses = array([ Bivarg( mu = array([xobs[i] + dx[i], yobs[i] + dy[i]]),
                                   sigma = objectsA[i].sigma + objectsB[i].sigma )
                           for i in range(len(objectsA)) ])
        draw_objects(ellipses, replot='yes')

        # Residuals
        quiver(xobs,yobs,dx,dy,color='r', angles='xy', scale_units='xy', scale=1)
    ax.autoscale(enable=None, axis='both', tight=True)
    show()

def plotSurface(pt, td, winds, map, stride, title, file_name):
    pylab.figure()
    pylab.axes((0.05, 0.025, 0.9, 0.9))

    u, v = winds

    nx, ny = pt.shape
    gs_x, gs_y = goshen_3km_gs
    xs, ys = np.meshgrid(gs_x * np.arange(nx), gs_y * np.arange(ny))   

    data_thin = tuple([ slice(None, None, stride) ] * 2)

    td_cmap = matplotlib.cm.get_cmap('Greens')
    td_cmap.set_under('#ffffff')
    pylab.contourf(xs, ys, td, levels=np.arange(40, 80, 5), cmap=td_cmap)
    pylab.colorbar()
    CS = pylab.contour(xs, ys, pt, colors='r', linestyles='-', linewidths=1.5, levels=np.arange(288, 324, 4))
    pylab.clabel(CS, inline_spacing=0, fmt="%d K", fontsize='x-small')
    pylab.quiver(xs[data_thin], ys[data_thin], u[data_thin], v[data_thin])

    drawPolitical(map, scale_len=75)

    pylab.suptitle(title)
    pylab.savefig(file_name)
    pylab.close()
    return

 def vectorField(self,plt,X,Y,U,V,title):
     if self.plotFields:
         # Create mask corresponding to 0 fluid velocity values inside obstructions
         M = zeros([X.quiverLength,Y.quiverLength],dtype='bool')
         M = (U.quiver == 0)
         
         # Mask the obstructions in the fluid velocity vector field
         U.quiver = ma.masked_array(U.quiver,mask=M)
         V.quiver = ma.masked_array(V.quiver,mask=M)
         
         # Build and scale the plot
         fig = plt.figure()        
         ax = fig.add_subplot(111)
         ax.set_ylim(Y.minValue, Y.maxValue)
         ax.set_xlim(X.minValue, X.maxValue)
         pos1 = ax.get_position()  
         pos2 = [pos1.x0+((pos1.width - (pos1.width*self.scaleX))*.5), pos1.y0+((pos1.height - (pos1.height*self.scaleY))*.5),  pos1.width*self.scaleX, pos1.height*self.scaleY] 
         ax.set_position(pos2)
         title = title + ' Vector Field'        
         plt.title(title)
         plt.xlabel("X [m]")
         plt.ylabel("Y [m]")
         plt.grid()
         
         quiver(X.quiver,Y.quiver,U.quiver,V.quiver)        
     
     
 
             

def plot_stiffness_field(k_cart,plottitle=''):
    n_points = 20
    ang_step = math.radians(360)/n_points
    x_list = []
    y_list = []
    u_list = []
    v_list = []
    k_cart = k_cart[0:2,0:2]

    for i in range(n_points):
        ang = i*ang_step
        for r in [0.5,1.,1.5]:
            dx = r*math.cos(ang)
            dy = r*math.sin(ang)
            dtau = -k_cart*np.matrix([dx,dy]).T
            x_list.append(dx)
            y_list.append(dy)
            u_list.append(dtau[0,0])
            v_list.append(dtau[1,0])
    
    pl.figure()
#    mpu.plot_circle(0,0,1.0,0.,math.radians(360))
    mpu.plot_radii(0,0,1.5,0.,math.radians(360),interval=ang_step,color='r')
    pl.quiver(x_list,y_list,u_list,v_list,width=0.002,color='k',scale=None)
    pl.axis('equal')
    pl.title(plottitle)

def main():
    # Create the grid
    x = arange(-100, 101)
    y = arange(-100, 101)

    # Create the meshgrid
    Y, X = meshgrid(x, y)
    A = 1
    B = 2
    V = 6*pi / 201
    W = 4*pi / 201
    F = A*sin(V*X) + B*cos(W*Y)
    Fx = V*A*cos(V*X)
    Fy = W*B*-sin(W*Y)

    # Show the images
    show_image(F)
    show_image(Fx)
    show_image(Fy)

    # Create the grid for the quivers
    xs = arange(-100, 101, 10)
    ys = arange(-100, 101, 10)

    # Here we determine the direction of the quivers
    Ys, Xs = meshgrid(ys, xs)
    FFx = V*A*cos(V*Xs)
    FFy = W*B*-sin(W*Ys)

    # Draw the quivers and the image
    clf()
    imshow(F, cmap=cm.gray, extent=(-100, 100, -100, 100))
    quiver(ys, xs, -FFy, FFx, color='red')
    show()

        def doSubplot(multiplier=1.0, layout=(-1, -1)):
            if time_sec < 16200:
                xs, ys = xs_1, ys_1
                domain_bounds = bounds_1sthalf
                grid = grid_1
            else:
                xs, ys = xs_2, ys_2
                domain_bounds = bounds_2ndhalf
                grid = grid_2

            try:
                mo = ARPSModelObsFile("%s/%s/KCYS%03dan%06d" % (base_path, exp, min_ens, time_sec))
            except AssertionError:
                mo = ARPSModelObsFile("%s/%s/KCYS%03dan%06d" % (base_path, exp, min_ens, time_sec), mpi_config=(2, 12))
            except:
                print "Can't load reflectivity ..."
                mo = {'Z':np.zeros((1, 255, 255), dtype=np.float32)}

            pylab.contour(xs, ys, wind[exp]['w'][wdt][domain_bounds], levels=np.arange(2, 102, 2), styles='-', colors='k')
            pylab.contour(xs, ys, wind[exp]['w'][wdt][domain_bounds], levels=np.arange(-100, 0, 2), styles='--', colors='k')

            pylab.quiver(xs[thin], ys[thin], wind[exp]['u'][wdt][domain_bounds][thin], wind[exp]['v'][wdt][domain_bounds][thin])

            pylab.contourf(xs, ys, mo['Z'][0][domain_bounds], levels=np.arange(10, 85, 5), cmap=NWSRef, zorder=-10)

            grid.drawPolitical(scale_len=10)

            row, col = layout
            if col == 1:
                pylab.text(-0.075, 0.5, exp_names[exp], transform=pylab.gca().transAxes, rotation=90, ha='center', va='center', size=12 * multiplier)

def main(argv=None):
  if argv is None:
    argv=sys.argv
  iresfile=argv[1];
  try:
    ct=float(argv[2]);
  except:
    ct=0.0
  ires=adore.res2dict(iresfile);
  A=ires['fine_coreg']['results'];
  A=A[A[:,5]>ct,:]
  az=int(ires['fine_coreg']['Initial offsets'].split(',')[0]);
  rg=int(ires['fine_coreg']['Initial offsets'].split(',')[1]);
#figure1
  pylab.figure()
  pylab.title('Fine Correlation Results')
  pylab.quiver(A[:,2], A[:,1], A[:,4], A[:,3], A[:,5]);
  pylab.colorbar()
#figure2

  pylab.figure()
  pylab.title('Fine Correlation Deviations')
  pylab.quiver(A[:,2], A[:,1], A[:,4]-rg, A[:,3]-az, A[:,5]);
  pylab.colorbar()
  pylab.show()

def quiver_image(X, Y, U, V):
    pylab.figure(1)
    pylab.quiver(X, Y, U, V)
    canvas = pylab.get_current_fig_manager().canvas
    canvas.draw()
    pil_image = Image.fromstring('RGB', canvas.get_width_height(), canvas.tostring_rgb())
    return pil_image

def whiskerplot(shear,dRA=1.,dDEC=1.,scale=5, combine=1, offset=(0,0) ):
    if combine>1:
        s = (combine*int(shear.shape[0]/combine),
             combine*int(shear.shape[1]/combine))
        shear = shear[0:s[0]:combine, 0:s[1]:combine] \
                + shear[1:s[0]:combine, 0:s[1]:combine] \
                + shear[0:s[0]:combine, 1:s[1]:combine] \
                + shear[1:s[0]:combine, 1:s[1]:combine]
        shear *= 0.25

        dRA *= combine
        dDEC *= combine

    
    theta = shear**0.5
    RA = offset[0] + np.arange(shear.shape[0])*dRA
    DEC = offset[1] + np.arange(shear.shape[1])*dDEC

    pylab.quiver(RA,DEC,
                 theta.real.T,theta.imag.T,
                 pivot = 'middle',
                 headwidth = 0,
                 headlength = 0,
                 headaxislength = 0,
                 scale=scale)
    pylab.xlim(0,shear.shape[0]*dRA)
    pylab.ylim(0,shear.shape[1]*dDEC)
    pylab.xlabel('RA (arcmin)')
    pylab.ylabel('DEC (arcmin)')

def ode_slopes_2states(func, radii_list, theta_deg_list, time_list):
    """
    Plot field of arrows indicating derivatives of the state
    :param func:
    :param radii_list:
    :param theta_deg_list:
    :param time_list:
    :return:
    """
    # to convert angle unit from degree to radian
    deg2rad = np.pi / 180
    # list of angles in radian
    theta_rad_list = tuple([(theta_deg * deg2rad) for theta_deg in theta_deg_list])
    # radii x angles mesh grid
    radii_mesh, theta_rad_mesh = np.meshgrid(radii_list, theta_rad_list)
    # polar coordinate to cartesian coordinate
    y = radii_mesh * np.cos(theta_rad_mesh), radii_mesh * np.sin(theta_rad_mesh)
    # derivatives of state at each point
    y_dot = func(y, time_list)
    # color
    color_mesh = np.sqrt(y_dot[0] * y_dot[0] + y_dot[1] * y_dot[1])
    # 1
    pylab.figure()
    pylab.quiver(y[0], y[1], y_dot[0], y_dot[1], color_mesh, angles='xy')
    l, r, b, t = pylab.axis()
    x_span, y2_mesh = r - l, t - b
    pylab.axis([l - 0.05 * x_span, r + 0.05 * x_span, b - 0.05 * y2_mesh, t + 0.05 * y2_mesh])
    pylab.axis('equal')
    pylab.grid()

def plotForceField(mouse,points,delta=0.,cm=0.):
    N=51
    sz=12
    rng=np.linspace(-sz,sz,N)
    R=np.zeros((N,N,2))
    for x in range(rng.size):
        for y in range(rng.size):
            offset=np.array([np.cos(points[:,2]),np.sin(points[:,2])]).T
            res=computeForce([rng[x],rng[y]],points=points[:,:2]+delta*offset,
                             maxmag=0.1,circlemass=cm)
            R[x,y,:]=res#np.linalg.norm(res)
    plt.cla()
    c=plt.Circle((0,0),radius=11.75,fill=False)
    plt.gca().add_patch(c)
    plt.plot(points[:,0],points[:,1],'ks',ms=8)
    plt.plot(mouse[0],mouse[1],'go',ms=8)
    plt.xlim([-12,12])
    plt.ylim([-12,12])
    plt.gca().set_aspect(1)
    #plt.pcolormesh(rng,rng,R.T,vmax=0.1)
    #R=np.square(R)
    plt.quiver(rng,rng,R[:,:,0].T,R[:,:,1].T,np.linalg.norm(R,axis=2).T,scale=3)
    #plt.pcolormesh(rng,rng,np.linalg.norm(R,axis=2).T)
    loc,minforce,ms=findNearestLocalMin(points[:,:2],start=np.copy(mouse[:2]),
                                        circlemass=cm)
    plt.plot(loc[0],loc[1],'bd',ms=8)
    plt.grid(b=False)
    ax=plt.gca()
    ax.spines['top'].set_visible(True)
    ax.spines['right'].set_visible(True)
    ax.set_xticklabels([]);ax.set_yticklabels([])

def main():
    x0, x1 = -0.5, 0.5
    y0, y1 = -0.5, 0.5
    h = .125
    nu = 1e3
    dt = 0.01
    t0 = 0.01

    plot_rows, plot_cols = 3, 3
    plot_every = 10

    x, y = init_position.triangular(x0, x1, y0, y1, cell_size=h)

    # initial vorticity and circulation
    vort = problems.lamb_oseen.vorticity(x, y, t0, nu=nu)
    circ = h**2 * vort

    # particle colors
    from itertools import cycle, izip as zip
    colors = ('#FF0000 #00FF00 #0000FF #FFFF00 #FF00FF #00FFFF '
              '#660000 #006600 #000066 #666600 #660066 #006666 '
             ).split()
    colors = [c for (c, _) in zip(cycle(colors), x)]

    t = t0
    iteration = 0
    plot_count = 0
    while True:
        plot_now = (iteration % plot_every == 0)

        if plot_now:
            plot_count += 1
            pylab.subplot(plot_rows, plot_cols, plot_count,
                          autoscale_on=False,
                          xlim=(1.2 * x0, 1.2 * x1),
                          ylim=(1.2 * y0, 1.2 * y1))
            pylab.scatter(x, y, s=2, c=colors, edgecolors='none')
            pylab.grid(True)

        u, v = vm.eval_velocity(x, y, circ)

        #pylab.quiver(x[::10], y[::10], u[::10], v[::10])

        if plot_now:
            pylab.quiver(x, y, u, v, color=colors)
            pylab.title('$t = %.2f$' % t)

        # convect
        x += u * dt
        y += v * dt

        iteration += 1
        t += dt

        if plot_count == plot_rows * plot_cols:
            break

    pylab.suptitle(ur'Lamb Oseen vortex, $\nu = %.3f$' % nu)
    pylab.show()

    def OnCalcShiftmap(self, event):
        from PYME.Analysis import twoColour, twoColourPlot
        import pylab
        masterChan = self.chChannel.GetSelection()
        
        master = self.objFitRes[masterChan]
        x0 = master['fitResults']['x0']
        y0 = master['fitResults']['y0']
        err_x0 = master['fitError']['x0']
        err_y0 = master['fitError']['y0']
        z0 = master['fitResults']['z0']

        wxy = master['fitResults']['wxy']

        wxy_bead = float(self.tBeadWXY.GetValue())

        mask = numpy.abs(wxy - wxy_bead) < (.25*wxy_bead)
        
        self.shiftfields ={}

        pylab.figure() 
        
        nchans = self.image.data.shape[3]
        ch_i = 1
        
        for ch in range(nchans):
            if not ch == masterChan:
                res = self.objFitRes[ch]
                
                x = res['fitResults']['x0']
                y = res['fitResults']['y0']
                z = res['fitResults']['z0']
                err_x = numpy.sqrt(res['fitError']['x0']**2 + err_x0**2)
                err_y = numpy.sqrt(res['fitError']['y0']**2 + err_y0**2)
                
                dx = x - x0
                dy = y - y0
                dz = z - z0
                
                print(('dz:', numpy.median(dz[mask])))
                
                
                spx, spy = twoColour.genShiftVectorFieldLinear(x[mask], y[mask], dx[mask], dy[mask], err_x[mask], err_y[mask])
                self.shiftfields[ch] = (spx, spy, numpy.median(dz[mask]))
                #twoColourPlot.PlotShiftField2(spx, spy, self.image.data.shape[:2])
                
                pylab.subplot(1,nchans -1, ch_i)
                ch_i += 1
                twoColourPlot.PlotShiftResidualsS(x[mask], y[mask], dx[mask], dy[mask], spx, spy)
                
        pylab.figure()
        X, Y = numpy.meshgrid(numpy.linspace(0., 70.*self.image.data.shape[0], 20), numpy.linspace(0., 70.*self.image.data.shape[1], 20))
        X = X.ravel()
        Y = Y.ravel()
        for k in self.shiftfields.keys():
            spx, spy, dz = self.shiftfields[k]
            pylab.quiver(X, Y, spx.ev(X, Y), spy.ev(X, Y), color=['r', 'g', 'b'][k], scale=2e3)
            
        pylab.axis('equal')

def ContourPlot2D(u, v, p, Y, X):
    pl.figure(figsize = (11,7), dpi = 100)
    pl.contourf(X,Y,p,alpha=0.5,cmap=cm.gist_heat)# plotting the pressure field contours
    pl.colorbar()
    pl.quiver(X[::2,::2],Y[::2,::2],u[::2,::2],v[::2,::2]) # plotting velocity vectors
    pl.xlabel('X')
    pl.ylabel('Y')
    pl.title('Pressure contours and velocity vectors')

	def quiver(self):
		mins = self.aus.min(axis=0)
		maxs = self.aus.max(axis=0)
		X,Y = np.meshgrid(linspace(mins[0], maxs[0], self.quiver_res), 
								linspace(mins[1], maxs[1], self.quiver_res))
		Z = np.dstack([X,Y])
		vals = self.f(0,Z.transpose(2,0,1))
		PL.quiver(X,Y,vals[0], vals[1])

def main():
    nu = 5e-4     # ν
    gamma0 = 1.0  # Γ₀
    vortex = LambOseenVortex(total_circulation=gamma0, viscosity=nu)

    # particle initialization
    x0, x1 = -0.5, 0.5
    y0, y1 = -0.5, 0.5
    h = .125 / 4
    e = 2 * h     # ε
    x, y = init_position.triangular(x0, x1, y0, y1, cell_size=h)

    # integration parameters
    t0 = 4
    dt = 0.001      # δt
    nr_iterations = 100

    # initial circulation evaluation
    initial_vorticity = vortex(x, y, t0)
    circ = h**2 * initial_vorticity

    iteration = 0
    t = t0
    while True:

        # convection
        u, v = eval_velocity(x, y, circ, e ** 2)
        x += u * dt
        y += v * dt

        # diffusion
        dcirc = eval_circulation_change(x, y, circ, nu, e ** 2)
        circ += dcirc * dt

        if iteration % 10 == 0:
            err2 = squared_velocity_errors(x, y, u, v, t, vortex)
            print "Iteration {0}, t = {1}".format(iteration, t)
            print "Squared error: {0} ± {1}".format(err2.mean(), err2.std())
            print

            #pylab.subplot(2, 5, (iteration + 1) // 10)
            #u_a, v_a = vortex.velocity(x, y, t)
            #pylab.scatter(x, y, s=1, c='red', edgecolor='none')
            #pylab.quiver(x, y, u, v, color='red')
            #pylab.quiver(x, y, u_a, v_a, color='blue')

        if not (iteration < nr_iterations):
            break

        iteration += 1
        t += dt

    u_a, v_a = vortex.velocity(x, y, t)
    pylab.scatter(x, y, s=1, c='red', edgecolor='none')
    pylab.quiver(x, y, u, v, color='red')
    pylab.quiver(x, y, u_a, v_a, color='blue')

    pylab.show()

    def PlotVel(self, **Scatter_args):
        """Plot the current configuration."""

        X, U = self.X, self.U

        m.quiver(X[:, 0], X[:, 1], U[:, 0], U[:, 1], **Scatter_args)

        m.xlim(0, 1)
        m.ylim(0, 1)

def plot_residuals_year_2pos(year, filt, pos1, pos2, refresh=False):
    year = str(year) 
    
    dir_xym = year + '_' + filt + '/01.XYM/'

    # Try to read from a FITS table of compiled data if possible.
    if refresh == True:
        make_residuals_table_year_2pos(year, filt, pos1, pos2)

    d = read_residuals_table_year_2pos(year, filt, pos1, pos2)

    print d.xstd_p.min(), d.xstd_p.max()

    # Remember, this was all aligned to F814W in ACS. So there is a
    # relative scale change between the transformed and the raw pixels.
    # We need to put evertyhing back on the same scale.
    # The plate scale was 0.411095 (pulled from TRANS.xym2mat.ref4)
    scale = 0.411095

    dx1 = d.xmean_p[0, :] - d.stars['x']
    dx2 = d.xmean_p[1, :] - d.stars['x']
    dy1 = d.ymean_p[0, :] - d.stars['y']
    dy2 = d.ymean_p[1, :] - d.stars['y']
    lim = 0.06
    idx = np.where((np.abs(dx1) < lim) & (np.abs(dy1) < lim) & 
                   (np.abs(dx2) < lim) & (np.abs(dy2) < lim))[0]

    # Plot the offsets for each position on a common coordinate system.
    py.clf()
    q1 = py.quiver(d.xmean[idx], d.ymean[idx], dx1[idx], dy1[idx], 
                   color='red', scale=1.8)
    q1 = py.quiver(d.xmean[idx], d.ymean[idx], dx2[idx], dy2[idx], 
                   color='blue', scale=1.8)
    py.quiverkey(q1, -0.1, -0.12, 0.02/scale, '0.02 WFC3IR pixels')
    py.axis('equal')
    py.xlabel("X (pixels)")
    py.ylabel("Y (pixels)")
    py.title(year + ',' + filt + ',' + pos1 + ',' + pos2)
    py.savefig('plots/resid_final_{0}_{1}_{2}_{3}.png'.format(year, filt, pos1, pos2))

    # Plot the offsets for each position on the raw coordinate system.
    py.clf()
    xraw1 = d.xraw[0,:,:].mean(axis=0)
    yraw1 = d.yraw[0,:,:].mean(axis=0)
    xraw2 = d.xraw[1,:,:].mean(axis=0)
    yraw2 = d.yraw[1,:,:].mean(axis=0)
    
    q1 = py.quiver(xraw1[idx], yraw1[idx], dx1[idx], dy1[idx], 
                   color='red', scale=1.8)
    q1 = py.quiver(xraw2[idx], yraw2[idx], dx2[idx], dy2[idx], 
                   color='blue', scale=1.8)
    py.quiverkey(q1, -0.1, -0.12, 0.02/scale, '0.02 WFC3IR pixels')
    py.axis('equal')
    py.xlabel("X (pixels)")
    py.ylabel("Y (pixels)")
    py.title(year + ',' + filt + ',' + pos1 + ',' + pos2)
    py.savefig('plots/resid_raw_{0}_{1}_{2}_{3}.png'.format(year, filt, pos1, pos2))