Python pylab.quiverkey函数代码示例

def plot_vector_diff_F160W():
    """
    Using the xym2mat analysis on the F160W filter in the 2010 data set,
    plot up the positional offset vectors for each reference star in
    each star list relative to the average list. This should show
    us if there are systematic problems with the distortion solution
    or PSF variations.
    """
    x, y, m, xe, ye, me, cnt = load_catalog()

    dx = x - x[0]
    dy = y - y[0]
    dr = np.hypot(dx, dy)

    py.clf()
    py.subplots_adjust(left=0.05, bottom=0.05, right=0.95, top=0.95)

    for ii in range(1, x.shape[0]):
        idx = np.where((x[ii,:] != 0) & (y[ii,:] != 0) & (dr[ii,:] < 0.1) &
                       (xe < 0.05) & (ye < 0.05))[0]

        py.clf()
        q = py.quiver(x[0,idx], y[0,idx], dx[ii, idx], dy[ii, idx], scale=1.0)
        py.quiverkey(q, 0.9, 0.9, 0.03333, label='2 mas', color='red')
        py.title('{0} stars in list {1}'.format(len(idx), ii))
        py.xlim(0, 4500)
        py.ylim(0, 4500)

        foo = input('Continue?')
        if foo == 'q' or foo == 'Q':
            break

def check_vpd_ks2_astrometry():
    """
    Check the VPD and quiver plots for our KS2-extracted, re-transformed astrometry.
    """
    catFile = workDir + '20.KS2_PMA/wd1_catalog.fits'
    tab = atpy.Table(catFile)

    good = (tab.xe_160 < 0.05) & (tab.ye_160 < 0.05) & \
        (tab.xe_814 < 0.05) & (tab.ye_814 < 0.05) & \
        (tab.me_814 < 0.05) & (tab.me_160 < 0.05)

    tab2 = tab.where(good)

    dx = (tab2.x_160 - tab2.x_814) * ast.scale['WFC'] * 1e3
    dy = (tab2.y_160 - tab2.y_814) * ast.scale['WFC'] * 1e3

    py.clf()
    q = py.quiver(tab2.x_814, tab2.y_814, dx, dy, scale=5e2)
    py.quiverkey(q, 0.95, 0.85, 5, '5 mas', color='red', labelcolor='red')
    py.savefig(workDir + '20.KS2_PMA/vec_diffs_ks2_all.png')

    py.clf()
    py.plot(dy, dx, 'k.', ms=2)
    lim = 30
    py.axis([-lim, lim, -lim, lim])
    py.xlabel('Y Proper Motion (mas)')
    py.ylabel('X Proper Motion (mas)')
    py.savefig(workDir + '20.KS2_PMA/vpd_ks2_all.png')

    idx = np.where((np.abs(dx) < 10) & (np.abs(dy) < 10))[0]
    print('Cluster Members (within dx < 10 mas and dy < 10 mas)')
    print(('   dx = {dx:6.2f} +/- {dxe:6.2f} mas'.format(dx=dx[idx].mean(),
                                                        dxe=dx[idx].std())))
    print(('   dy = {dy:6.2f} +/- {dye:6.2f} mas'.format(dy=dy[idx].mean(),
                                                        dye=dy[idx].std())))

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))

    def PlotBaseFly(self):
        self.BaseFly()
        # Creates lists of x and y coordinates
        x_coords = [0, (self.displ * numpy.cos(self.head * (2 * numpy.pi) / 360))]
        y_coords = [0, (self.displ * numpy.sin(self.head * (2 * numpy.pi) / 360))]

        fig = pylab.figure()
        # Adjusts margins
        fig.subplots_adjust(left=None, bottom=0.2, right=None, top=None, wspace=None, hspace=None)

        # Plots start point
        pylab.plot(x_coords[0], y_coords[0], self.marker1, markersize=10, label='Start')
        # Plots end point
        pylab.plot(x_coords[1], y_coords[1], self.marker2, markersize=10, label='End')

        # Generates 2D wind vector field
        if self.displ < 10:
            X, Y = numpy.meshgrid(numpy.arange(-30, 30, 6), numpy.arange(-30, 30, 6))
            Q = pylab.quiver(X, Y, self.wind_U, self.wind_V, label='Wind Direction')
            QK = pylab.quiverkey(Q, 0.1, -0.08, 10, 'Wind Direction', labelpos='W',
                                 fontproperties={'weight': 'bold'})
        else:
            X, Y = numpy.meshgrid(numpy.arange(-self.displ - 20, self.displ + 20, self.displ / 5), numpy.arange(-self.displ - 20, self.displ + 20, self.displ / 5))
            Q = pylab.quiver(X, Y, self.wind_U, self.wind_V, label='Wind Direction')
            QK = pylab.quiverkey(Q, 0.1, -0.08, self.displ / 5, 'Wind Direction', labelpos='W',
                                 fontproperties={'weight': 'bold'})

        # Shows displ, head, time
        # Refer to http://docs.python.org/2/library/string.html for formatting
        # and http://stackoverflow.com/questions/19261089/how-to-use-new-style-string-formatting-in-matplotlib-figure-text
        '''
        Text = '{0} {1:10}\n'.format('Displacement: ', int(self.displ)) + \
               '{0} {1:10}\n'.format('Heading:         ', int(self.head)) + \
               '{0} {1:10}\n'.format('Flight Time:  ', int(self.time))

        fig.text(0.25, -0.01, Text, family='monospace')
        '''

        TextL = 'Displacement (m):\nHeading (Degrees):\nDescent Time (s):'
        TextR = str(int(self.displ)) + '\n' + str(int(self.head)) + '\n' + str(int(self.time))
        fig.text(0.4, 0.01, TextL, multialignment='left')
        fig.text(0.6, 0.01, TextR, multialignment='right')

        pylab.title('UAV Displacement in Y against UAV Displacement in X')
        pylab.xlabel('UAV Displacement in X')
        pylab.ylabel('UAV Displacement in Y')
        pylab.xlim(-self.displ - 20, self.displ + 20)
        pylab.ylim(-self.displ - 20, self.displ + 20)

        pylab.legend()

        pylab.show()

def show_shwfs_vecs(
    shiftvec, subaps, refpos=None, img=None, extent=None, title=None, scale=10, pause=False, fignum=None, keep=True
):
    """
	Show a SHWFS measurements with subapertures and potentially a background 
	image. If **refpos** is given, SHWFS shift vectors will be plotted at 
	those locations, otherwise at the center of the sub apertures.

	Additionally, a background image (wavefront, sensor image) can be 
	plotted along with the measurements in the background. If the coordinate 
	system is different than that of the sub apertures, this can be fixed by 
	using the **extent** keyword.

	**Scale** tweaks the arrow length, unity being the approximate real size.
	The default is 10 to enhance the visibility.

	@param [in] shiftvec (N,2) vector with shifts
	@param [in] subaps (N,4) vector with sub aperture positions as (low0, high0, low1, high1)
	@param [in] refpos Positions to plot shift vectors at (if None use subaps center)
	@param [in] img Background image, e.g. a wavefront
	@param [in] extent Extent for the image for imshow()
	@param [in] title Plot title
	@param [in] scale Scale used for arrow width, unity is approximately
	@param [in] pause Pause before continuing
	@param [in] fignum Figure number to use for figure()
	@param [in] keep Keep plot window open after returning
	"""

    import pylab as plt
    from matplotlib.collections import PatchCollection

    sasize = np.r_[subaps[:, :2].ptp(1).mean(), subaps[:, 2:].ptp(1).mean()]
    subaps_im = [
        plt.Rectangle((subap[1], subap[0]), sasize[0], sasize[1], fc="none", ec="k") for subap in subaps[:, ::2]
    ]

    plt.figure(fignum)
    plt.clf()
    plt.title(title or "SHWFS vectors on MLA grid")

    if img != None:
        plt.imshow(img, extent=extent)

    if refpos == None:
        refpos = subaps[:, ::2] + sasize / 2

    q = plt.quiver(
        refpos[:, 1], refpos[:, 0], shiftvec[:, 1], shiftvec[:, 0], angles="xy", scale=refpos.ptp() / 10.0, color="r"
    )
    p = plt.quiverkey(q, refpos.min(0)[1], refpos.min(0)[0], 5, "5 pix.", coordinates="data", color="r")

    thisgca = plt.gca()
    thisgca.add_collection(PatchCollection(subaps_im, match_original=True))

    if pause:
        raw_input("Press any key to continue...")
    if not keep:
        plt.close()

def test_darPlusDistortion():
    data_dir = module_dir + '/distortion/'
    file_geox_darunfix = data_dir + 'nirc2dist_xgeoim.fits'
    file_geoy_darunfix = data_dir + 'nirc2dist_ygeoim.fits'

    data_dir = '/u/ghezgroup/data/m92_test/08jul_new_on/'
    file_geox_darfix = data_dir + 'reduce/kp/gc_f1/ce0249geo_x.fits'
    file_geoy_darfix = data_dir + 'reduce/kp/gc_f1/ce0249geo_y.fits'

    xon = pyfits.getdata(file_geox_darfix)
    yon = pyfits.getdata(file_geoy_darfix)
    xoff = pyfits.getdata(file_geox_darunfix)
    yoff = pyfits.getdata(file_geoy_darunfix)

    # Make arrays with the coordinates for each 
    imgsize = 1024
    axisX = np.arange(imgsize, dtype=float)
    axisY = np.arange(imgsize, dtype=float)
    xcoo2d, ycoo2d = np.meshgrid(axisX, axisY)

    # Lets trim so that we only keep every 20th pixel
    idx = np.arange(25, imgsize, 50)
    xon = xon.take(idx, axis=0).take(idx, axis=1)
    yon = yon.take(idx, axis=0).take(idx, axis=1)
    xoff = xoff.take(idx, axis=0).take(idx, axis=1)
    yoff = yoff.take(idx, axis=0).take(idx, axis=1)
    xcoo2d = xcoo2d.take(idx, axis=0).take(idx, axis=1)
    ycoo2d = ycoo2d.take(idx, axis=0).take(idx, axis=1)

    # Calculate differences
    xdiff = xon - xoff
    ydiff = yon - yoff

    # Make vector plots
    py.clf()
    qvr = py.quiver2([xcoo2d], [ycoo2d], [xdiff], [ydiff],
                     units='width', scale=5, 
                     width=0.005, headwidth=3, headlength=3, 
                     headaxislength=3)
    py.quiverkey(qvr, 100, 1120, 1.0, '1 pixel', coordinates='data', color='r')
    py.xlabel('NIRC2 X (pixel)')
    py.ylabel('NIRC2 Y (pixel)')
    py.title('Arrows point to DAR Fix')
    #py.savefig('plots/vector_daroffon.png')
    py.show()

def plot_arches_pm():
    """
    Plot the proper motion vectors for the Arches (Clarkson+ 2012).
    """
    dataFile = "/u/jlu/data/arches/clarkson2012_table5.txt"

    d = atpy.Table(dataFile, type="mrt")

    # Put into a "field" reference frame.
    idx = np.where(d.FMem > 0.9)[0]
    f_pmX = d.pmX[idx].mean()
    f_pmY = d.pmY[idx].mean()

    d.pmX -= f_pmX
    d.pmY -= f_pmY

    py.clf()
    Q = py.quiver(d.DelX, d.DelY, d.pmX, d.pmY)
    py.quiverkey(Q, 95, 95, 1, "40 km/s", coordinates="figure")

def plot_vpd_align(alignRoot, epoch1, epoch2, outRoot='plot_align'):
    """
    Read in an aligned data set. Plot a vector-point-diagram for the 
    positional differences between two different epochs within the data set.
    """
    s = starset.StarSet(alignRoot)

    x1 = s.getArrayFromEpoch(epoch1, 'xorig')
    y1 = s.getArrayFromEpoch(epoch1, 'yorig')
    x2 = s.getArrayFromEpoch(epoch2, 'xorig')
    y2 = s.getArrayFromEpoch(epoch2, 'yorig')

    x1e = s.getArrayFromEpoch(epoch1, 'xpixerr_p')
    y1e = s.getArrayFromEpoch(epoch1, 'ypixerr_p')
    x2e = s.getArrayFromEpoch(epoch2, 'xpixerr_p')
    y2e = s.getArrayFromEpoch(epoch2, 'ypixerr_p')

    dx = x2 - x1
    dy = y2 - y1
    dxe = np.hypot(x1e, x2e)
    dye = np.hypot(y1e, y2e)
    dr = np.hypot(dx, dy)

    idx = np.where((dxe < 0.02) & (dye < 0.02) & (dr < 3.0))[0]
    print len(x1), len(idx)

    drmax = 1.5

    py.figure(1)
    py.clf()
    py.plot(dx[idx], dy[idx], 'k.', ms=2)
    py.xlabel('Delta X (pixels)')
    py.ylabel('Delta Y (pixels)')
    py.axis([-drmax, drmax, -drmax, drmax])
    py.savefig(outRoot + '_vpd.png')

    py.figure(2)
    py.clf()
    q = py.quiver(x1[idx], y1[idx], dx[idx], dy[idx], scale=10)
    py.quiverkey(q, 0.5, 0.95, 0.05, '0.05 pix', color='red')
    py.savefig(outRoot + '_dr_vec.png')

def whisker4QReduce(X2WIN_IMAGE=None,Y2WIN_IMAGE=None,XYWIN_IMAGE=None):
    """
    This function make the whisker plot based on the sextractor export from QuickReduce
    J.Hao, 12/4/2012
    """
    xpos = np.genfromtxt('xpos_ypos_fp.txt').T[0]
    ypos = np.genfromtxt('xpos_ypos_fp.txt').T[1]
    temp = np.zeros(62).astype('complex')
    temp.real = X2WIN_IMAGE - Y2WIN_IMAGE
    temp.imag = 2*XYWIN_IMAGE
    data=np.array([xpos, ypos,X2WIN_IMAGE + Y2WIN_IMAGE,temp]).T
    data = averageN30new(data)
    data = subMeanAll(data)
    pl.figure(figsize=(11,5.5))
    pl.subplot(1,2,1)
    phi22 = 0.5*np.arctan2(data[:,3].imag,data[:,3].real)
    x = data[:,0].real
    y = data[:,1].real
    u = np.abs(data[:,3])*np.cos(phi22)
    v = np.abs(data[:,3])*np.sin(phi22)
    qvr = pl.quiver(x,y,u,v,width = 0.004, color='r',pivot='middle',headwidth=0.,headlength=0.,headaxislength=0.,scale_units='width')
    qk = pl.quiverkey(qvr, -150,-240,0.3,str(0.3)+' pix^2',coordinates='data',color='blue')
    pl.plot(x,y,'b,')
    pl.xlim(-250,250)
    pl.ylim(-250,250)
    pl.grid(color='g')
    pl.xlabel('Camera WEST [mm]')
    pl.ylabel('Camera NORTH [mm]')
    pl.title('M22')
    pl.subplot(1,2,2)
    m20sqr = np.sqrt(data[:,2].real)
    x = data[:,0].real
    y = data[:,1].real
    m20sqr_med = np.median(m20sqr) 
    m20sqr_diff = m20sqr - m20sqr_med
    m20sqr_diff_absmed = np.median(np.abs(m20sqr_diff))
    plotScale = 1./m20sqr_diff_absmed*100
    pos = m20sqr_diff >=0
    neg = m20sqr_diff < 0
    pl.scatter(x[pos],y[pos],s=m20sqr_diff[pos]*plotScale,c='r',alpha=0.5)
    pl.scatter(x[neg],y[neg],s=-m20sqr_diff[neg]*plotScale,c='b',alpha=0.5)
    pl.scatter(-230,-210,s=0.01*plotScale,c='b',alpha=0.5)
    pl.text(-200,-215,'-'+str(0.01)+' pix')
    pl.scatter(-230,-230,s=0.01*plotScale,c='r',alpha=0.5)
    pl.text(-200,-235,str(0.01)+' pix')
    pl.plot(x,y,'y,')
    pl.grid(color='g')
    pl.xlim(-250,250)
    pl.ylim(-250,250)
    pl.xlabel('Camera WEST [mm]')
    pl.ylabel('Camera NORTH [mm]')
    pl.title('median '+r'$\sqrt{M20}$: '+str(round(0.27*m20sqr_med,3))+' [arcsec]')
    return '---done!--'

def plot_pos_diff_xym1mat(xym1mat_cfile, scale=8, errMax=0.05):
    """
    Send in a match_c.txt file from xym1mat (3rd file) that has been trimmed
    of all bogus entries. Then plot of the differences in the positions
    from the 1st and 2nd starlists that went into the list matching process.
    This allows me to see any large scale systematics due to residual distortion
    or alignment errors.

    This assumes that the input starlists were in the XYMEEE format. We will
    use the frame #2 astrometric errors to trim out bad points.
    """
    t = atpy.Table(xym1mat_cfile, type='ascii')

    if errMax != None:
        x2e = t.col17
        y2e = t.col18

        err = np.hypot(x2e, y2e)

        origLen = len(t)
        t = t.where(err < errMax)
        newLen = len(t)

        print 'Trimmed %d of %d sources with errors > %.2f pix' % \
            (origLen-newLen, origLen, errMax)
        
    x1 = t.col8
    y1 = t.col9
    x2 = t.col11
    y2 = t.col12
    dx = x1 - x2
    dy = y1 - y2


    py.clf()
    q = py.quiver(x1, y1, dx, dy, scale=scale)
    py.quiverkey(q, 0.9, 0.95, 0.1, '0.1 pix', coordinates='axes', color='red')
    py.xlabel('X (pix)')
    py.ylabel('Y (pix)')
    py.savefig('pos_diff_vecs_' + xym1mat_cfile.replace('txt', 'png'))

  def sliceuv(self,ind,time=0,plot=False,**opts):
    savename=False
    if savename in opts.keys(): savename=opts['savename']

    if ind=='bar':
      isK,uname,vname,ind = True, 'ubar','vbar', 9999
    elif ind in ('s','surf','surface'):
      isK,uname,vname,ind = True,  'u','v', -1
    elif ind>=0:
      isK,uname,vname,ind = True,  'u','v', ind
    elif ind <0:
      isK=False
      isK,uname,vname,ind = False, 'u','v', ind


    if isK:
      xu,yu,zu,u=self.slicek(uname,ind,time)
      xv,yv,zv,v=self.slicek(vname,ind,time)
    else:
      xu,yu,zu,u=self.slicez(uname,ind,time)
      xv,yv,zv,v=self.slicez(vname,ind,time)

    z=(zu[1:,:]+zu[:-1,:])/2.
    u=( u[1:,:]+ u[:-1,:])/2.
    v=( v[:,1:]+ v[:,:-1])/2.

    x=self.grid.lonp
    y=self.grid.latp

    # rotate uv:
    ang=rt.rho2uvp(self.grid.angle,'p')
    u,v=calc.rot2d(u,v,-ang)


    if plot:
      p=self.grid.plot(bathy=None)
      xm, ym = p(x,y)
      mm=0*m
      mm[::3,::3]=1
      mm=mm*m==1
      s=np.sqrt(u**2+v**2)
      q=pylab.quiver(xm[mm],ym[mm],u[mm],v[mm],s[mm])
      pylab.colorbar(shrink=.7)
      pylab.axis([p.xmin, p.xmax, p.ymin, p.ymax])
      qk = pylab.quiverkey(q, .05, .01, 0.2, '0.2 m/s',labelpos='E',
                                              coordinates='axes')
      if savename:
        pylab.savefig(savename,dpi=pylab.gcf().dpi)
        pylab.close(pylab.gcf())

    return x,y,z,u,v

def plot_residuals_year_pos(year, filt, pos, refresh=False):
    year = str(year) 
    
    dir_xym = year + '_' + filt + '_' + pos + '/01.XYM/'

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

    d = read_residuals_table_year_pos(year, filt, pos)

    print d.xstd.min(), d.xstd.max()
    
    # Look at the variance in the raw X pixels.
    py.clf()
    q = py.quiver(d.xmean, d.ymean, d.xstd, d.ystd, scale=0.7)
    py.quiverkey(q, -0.1, -0.12, 0.02, '0.02 pixels', color='red')
    py.xlim(0, 1050)
    py.xlabel("X (pixels)")
    py.ylabel("Y (pixels)")
    py.title(year + ',' + filt + ',' + pos)
    
    py.savefig('plots/resid_errors_{0}_{1}_{2}.png'.format(year, filt, pos))

    def sliceuv(self, ind, time=0, plot=False, **opts):
        savename = False
        if savename in opts.keys():
            savename = opts["savename"]

        if ind == "bar":
            isK, uname, vname, ind = True, "ubar", "vbar", 9999
        elif ind in ("s", "surf", "surface"):
            isK, uname, vname, ind = True, "u", "v", -1
        elif ind >= 0:
            isK, uname, vname, ind = True, "u", "v", ind
        elif ind < 0:
            isK = False
            isK, uname, vname, ind = False, "u", "v", ind

        if isK:
            xu, yu, zu, u = self.slicek(uname, ind, time)
            xv, yv, zv, v = self.slicek(vname, ind, time)
        else:
            xu, yu, zu, u = self.slicez(uname, ind, time)
            xv, yv, zv, v = self.slicez(vname, ind, time)

        z = (zu[1:, :] + zu[:-1, :]) / 2.0
        u = (u[1:, :] + u[:-1, :]) / 2.0
        v = (v[:, 1:] + v[:, :-1]) / 2.0

        x = self.grid.lonp
        y = self.grid.latp

        # rotate uv:
        ang = rt.rho2uvp(self.grid.angle, "p")
        u, v = calc.rot2d(u, v, -ang)

        if plot:
            p = self.grid.plot(bathy=None)
            xm, ym = p(x, y)
            mm = 0 * m
            mm[::3, ::3] = 1
            mm = mm * m == 1
            s = np.sqrt(u ** 2 + v ** 2)
            q = pylab.quiver(xm[mm], ym[mm], u[mm], v[mm], s[mm])
            pylab.colorbar(shrink=0.7)
            pylab.axis([p.xmin, p.xmax, p.ymin, p.ymax])
            qk = pylab.quiverkey(q, 0.05, 0.01, 0.2, "0.2 m/s", labelpos="E", coordinates="axes")
            if savename:
                pylab.savefig(savename, dpi=pylab.gcf().dpi)
                pylab.close(pylab.gcf())

        return x, y, z, u, v

def display_2nd_moments(data=None):
    # remove the mean for all moments
    datamean = data.mean(axis = 0)
    data = subMeanAll(data)
    pl.figure(figsize=(11,5.5))
    pl.subplot(1,2,1)
    phi22 = 0.5*np.arctan2(data[:,3].imag,data[:,3].real)
    x = data[:,0].real
    y = data[:,1].real
    #phi22[x<0] = phi22+np.deg2rad(180)
    u = np.abs(data[:,3])*np.cos(phi22)
    v = np.abs(data[:,3])*np.sin(phi22)
    qvr = pl.quiver(x,y,u,v,width = 0.004, color='r',pivot='middle',headwidth=0.,headlength=0.,headaxislength=0.,scale_units='width')
    #qk = pl.quiverkey(qvr, -150,-240,np.max(np.sqrt(u**2+v**2)),str(round(np.max(np.sqrt(u**2+v**2)),3))+' pix^2',coordinates='data',color='blue')
    qk = pl.quiverkey(qvr, -150,-240,0.3,str(0.3)+' pix^2',coordinates='data',color='blue')
    pl.plot(x,y,'b,')
    pl.xlim(-250,250)
    pl.ylim(-250,250)
    pl.grid(color='g')
    pl.xlabel('Camera WEST [mm]')
    pl.ylabel('Camera NORTH [mm]')
    pl.title('mean |M22|= '+str(round(abs(datamean[:,3]),5))+' pix^2')
    pl.subplot(1,2,2)
    m20sqr = np.sqrt(data[:,2].real)
    x = data[:,0].real
    y = data[:,1].real
    m20sqr_med = np.median(m20sqr) 
    m20sqr_diff = m20sqr - m20sqr_med
    m20sqr_diff_absmed = np.median(np.abs(m20sqr_diff))
    plotScale = 1./m20sqr_diff_absmed*100
    pos = m20sqr_diff >=0
    neg = m20sqr_diff < 0
    pl.scatter(x[pos],y[pos],s=m20sqr_diff[pos]*plotScale,c='r',alpha=0.5)
    pl.scatter(x[neg],y[neg],s=-m20sqr_diff[neg]*plotScale,c='b',alpha=0.5)
    pl.scatter(-230,-210,s=0.01*plotScale,c='b',alpha=0.5)
    pl.text(-200,-215,'-'+str(0.01)+' pix')
    pl.scatter(-230,-230,s=0.01*plotScale,c='r',alpha=0.5)
    pl.text(-200,-235,str(0.01)+' pix')
    pl.plot(x,y,'y,')
    pl.grid(color='g')
    pl.xlim(-250,250)
    pl.ylim(-250,250)
    pl.xlabel('Camera WEST [mm]')
    pl.ylabel('Camera NORTH [mm]')
    pl.title('median '+r'$\sqrt{M20}$: '+str(round(scale*m20sqr_med,3))+' [arcsec]')
    return '---done!--'

def add_arrows(fieldx,fieldy):
	norm=abs(np.concatenate((fieldx,fieldy)).max())/1. #1d norm
	freq=4
	M = py.sqrt(pow(fieldx, 2) + pow(fieldy, 2))
	X,Y = py.meshgrid( py.arange(0,nx,1),py.arange(0,ny,1) )

        if(1==1):
            Q = py.quiver( 
		X[::freq, ::freq], Y[::freq, ::freq], 
		fieldx[::freq, ::freq], fieldy[::freq, ::freq], 
                #		M[::freq, ::freq],cmap=py.cm.jet,
                #		,color='yellow'
		#M[::freq, ::freq],cmap=py.cm.gray, alpha=0.2,
#		facecolors=mycmap(M[::freq, ::freq]),
                #		width=0.4,linewidths=(.5), edgecolors=('k'), 
#                color=colors,
		width=(0.2*py.sqrt(freq)),linewidths=(.2),
#		pivot='mid', 
#		pivot='tail', 
		units='x',
		headaxislength=5 ,
		scale=norm/freq)
            colors=mycmap_rgba(M[::freq, ::freq])
            Q.set_facecolor(colors)
            Q.set_edgecolor('k')
            Q.set_edgecolor(colors)
        else:
            Q = py.quiver( 
		X[::freq, ::freq], Y[::freq, ::freq], 
		fieldx[::freq, ::freq], fieldy[::freq, ::freq], 
                #		M[::freq, ::freq],cmap=py.cm.gray, alpha=1,
		width=0.2,
		units='x',
		headaxislength=5 ,
		scale=norm/2)

        legend='{:.2g}km/s'.format(norm)
        qk = py.quiverkey(Q, 0.9, 1.05, norm, legend,
                          labelpos='E',
                          fontproperties={'weight': 'bold'})