def make_vector_plot(coordfile, columns=[1, 2, 3, 4], data=None,
                     figure_id=None, title=None, axes=None, every=1,
                     labelsize=8, ylimit=None, limit=None, xlower=None,
                     ylower=None, output=None, headl=4, headw=3,
                     xsh=0.0, ysh=0.0, fit=None, scale=1.0, vector=True,
                     textscale=5, append=False, linfit=False, rms=True,
                     plotname=None):
    """ Convert a XYXYMATCH file into a vector plot or set of residuals plots.

        This function provides a single interface for generating either a
        vector plot of residuals or a set of 4 plots showing residuals.
        The data being plotted can also be adjusted for a linear fit
        on-the-fly.

        Parameters
        ----------
        coordfile : string
            Name of file with matched sets of coordinates. This input file can
            be a file compatible for use with IRAF's geomap.
        columns : list [Default: [0,1,2,3]]
            Column numbers for the X,Y positions from each image
        data : list of arrays
            If specified, this can be used to input matched data directly
        title : string
            Title to be used for the generated plot
        axes : list
            List of X and Y min/max values to customize the plot axes
        every : int [Default: 1]
            Slice value for the data to be plotted
        limit : float
            Radial offset limit for selecting which sources are included in
            the plot
        labelsize : int [Default: 8] or str
            Font size to use for tick labels, either in font points or as a
            string understood by tick_params().
        ylimit : float
            Limit to use for Y range of plots.
        xlower : float
        ylower : float
            Limit in X and/or Y offset for selecting which sources are included
            in the plot
        output : string
            Filename of output file for generated plot
        headl : int [Default: 4]
            Length of arrow head to be used in vector plot
        headw : int [Default: 3]
            Width of arrow head to be used in vector plot
        xsh : float
        ysh : float
            Shift in X and Y from linear fit to be applied to source positions
            from the first image
        scale : float
            Scale from linear fit to be applied to source positions from the
            first image
        fit : array
            Array of linear coefficients for rotation (and scale?) in X and Y
            from a linear fit to be applied to source positions from the
            first image
        vector : bool [Default: True]
            Specifies whether or not to generate a vector plot. If False, task
            will generate a set of 4 residuals plots instead
        textscale : int [Default: 5]
            Scale factor for text used for labelling the generated plot
        append : bool [Default: False]
            If True, will overplot new plot on any pre-existing plot
        linfit : bool [Default: False]
            If True, a linear fit to the residuals will be generated and
            added to the generated residuals plots
        rms : bool [Default: True]
            Specifies whether or not to report the RMS of the residuals as a
            label on the generated plot(s).
        plotname : str [Default: None]
            Write out plot to a file with this name if specified.

    """
    from matplotlib import pyplot as plt

    if data is None:
        data = readcols(coordfile, cols=columns)

    xy1x = data[0]
    xy1y = data[1]
    xy2x = data[2]
    xy2y = data[3]

    numpts = xy1x.shape[0]
    if fit is not None:
        xy1x, xy1y = apply_db_fit(data, fit, xsh=xsh, ysh=ysh)
        dx = xy2x - xy1x
        dy = xy2y - xy1y
    else:
        dx = xy2x - xy1x - xsh
        dy = xy2y - xy1y - ysh
    # apply scaling factor to deltas
    dx *= scale
    dy *= scale

    print('Total # points: {:d}'.format(len(dx)))
    if limit is not None:
        indx = np.sqrt(dx**2 + dy**2) <= limit
        dx = dx[indx].copy()
        dy = dy[indx].copy()
        xy1x = xy1x[indx].copy()
        xy1y = xy1y[indx].copy()

    if xlower is not None:
        xindx = np.abs(dx) >= xlower
        dx = dx[xindx].copy()
        dy = dy[xindx].copy()
        xy1x = xy1x[xindx].copy()
        xy1y = xy1y[xindx].copy()
    print('# of points after clipping: {:d}'.format(len(dx)))

    dr = np.sqrt(dx**2 + dy**2)
    max_vector = dr.max()

    if output is not None:
        write_xy_file(output, [xy1x, xy1y, dx, dy])

    fig = plt.figure(num=figure_id)
    if not append:
        plt.clf()

    if vector:
        dxs = imagestats.ImageStats(dx.astype(np.float32))
        dys = imagestats.ImageStats(dy.astype(np.float32))
        minx = xy1x.min()
        maxx = xy1x.max()
        miny = xy1y.min()
        maxy = xy1y.max()
        plt_xrange = maxx - minx
        plt_yrange = maxy - miny

        qplot = plt.quiver(xy1x[::every], xy1y[::every], dx[::every],
                           dy[::every], units='y', headwidth=headw,
                           headlength=headl)
        key_dx = 0.01 * plt_xrange
        key_dy = 0.005 * plt_yrange * textscale
        maxvec = max_vector / 2.
        key_len = round(maxvec + 0.005, 2)

        plt.xlabel('DX: %.4f to %.4f +/- %.4f' % (dxs.min, dxs.max,
                                                  dxs.stddev))
        plt.ylabel('DY: %.4f to %.4f +/- %.4f' % (dys.min, dys.max,
                                                  dys.stddev))
        plt.title(r"$Vector\ plot\ of\ %d/%d\ residuals:\ %s$" %
                  (xy1x.shape[0], numpts, title))
        plt.quiverkey(qplot, minx + key_dx, miny - key_dy, key_len,
                      "%0.2f pixels" % (key_len),
                      coordinates='data', labelpos='E', labelcolor='Maroon',
                      color='Maroon')
    else:
        plot_defs = [[xy1x, dx, "X (pixels)", "DX (pixels)"],
                     [xy1y, dx, "Y (pixels)", "DX (pixels)"],
                     [xy1x, dy, "X (pixels)", "DY (pixels)"],
                     [xy1y, dy, "Y (pixels)", "DY (pixels)"]]
        if axes is None:
            # Compute a global set of axis limits for all plots
            minx = min(xy1x.min(), xy1y.min())
            maxx = max(xy1x.max(), xy1y.max())
            miny = min(dx.min(), dy.min())
            maxy = max(dx.max(), dy.max())
        else:
            minx = axes[0][0]
            maxx = axes[0][1]
            miny = axes[1][0]
            maxy = axes[1][1]

        if ylimit is not None:
            miny = -ylimit
            maxy = ylimit

        rms_labelled = False
        if title is None:
            fig.suptitle("Residuals [%d/%d]" % (xy1x.shape[0], numpts),
                         ha='center', fontsize=labelsize + 6)
        else:
            # This definition of the title supports math symbols in the title
            fig.suptitle(r"$" + title + "$", ha='center',
                         fontsize=labelsize + 6)

        for pnum, p in enumerate(plot_defs):
            pn = pnum + 1
            ax = fig.add_subplot(2, 2, pn)
            plt.plot(
                p[0], p[1], 'b.',
                label='RMS(X) = %.4f, RMS(Y) = %.4f' % (dx.std(), dy.std())
            )
            lx = [int((p[0].min() - 500) / 500) * 500,
                  int((p[0].max() + 500) / 500) * 500]
            plt.plot(lx, [0.0, 0.0], 'k', linewidth=3)
            plt.axis([minx, maxx, miny, maxy])
            if rms and not rms_labelled:
                leg_handles, leg_labels = ax.get_legend_handles_labels()
                fig.legend(leg_handles, leg_labels, loc='center left',
                           fontsize='small', frameon=False,
                           bbox_to_anchor=(0.33, 0.51), borderaxespad=0)
                rms_labelled = True

            ax.tick_params(labelsize=labelsize)

            # Fine-tune figure; hide x ticks for top plots and y ticks for
            # right plots
            if pn <= 2:
                plt.setp(ax.get_xticklabels(), visible=False)
            else:
                ax.set_xlabel(plot_defs[pnum][2])

            if pn % 2 == 0:
                plt.setp(ax.get_yticklabels(), visible=False)
            else:
                ax.set_ylabel(plot_defs[pnum][3])

            if linfit:
                lxr = int((lx[-1] - lx[0]) / 100)
                lyr = int((p[1].max() - p[1].min()) / 100)
                a = np.vstack([p[0], np.ones(len(p[0]))]).T
                m, c = np.linalg.lstsq(a, p[1])[0]
                yr = [m * lx[0] + c, lx[-1] * m + c]
                plt.plot([lx[0], lx[-1]], yr, 'r')
                plt.text(
                    lx[0] + lxr, p[1].max() + lyr,
                    "%0.5g*x + %0.5g [%0.5g,%0.5g]" % (m, c, yr[0], yr[1]),
                    color='r'
                )

    plt.draw()

    if plotname:
        suffix = plotname[-4:]
        if '.' not in suffix:
            output += '.png'
            format = 'png'
        else:
            if suffix[1:] in ['png', 'pdf', 'ps', 'eps', 'svg']:
                format = suffix[1:]
        plt.savefig(plotname, format=format)