def _estimate_2dhist_shift(imgxy, refxy, searchrad=3.0):
    """ Create a 2D matrix-histogram which contains the delta between each
        XY position and each UV position. Then estimate initial offset
        between catalogs.
    """
    print("Computing initial guess for X and Y shifts...")

    # create ZP matrix
    zpmat = _xy_2dhist(imgxy, refxy, r=searchrad)

    nonzeros = np.count_nonzero(zpmat)
    if nonzeros == 0:
        # no matches within search radius. Return (0, 0):
        print("WARNING: No matches found within a search radius of {:g} "
              "pixels.".format(searchrad))
        return 0.0, 0.0, 0, 0, zpmat, False

    elif nonzeros == 1:
        # only one non-zero bin:
        yp, xp = np.unravel_index(np.argmax(zpmat), zpmat.shape)
        maxval = int(np.ceil(zpmat[yp, xp]))
        xp -= searchrad
        yp -= searchrad
        print("Found initial X and Y shifts of {:.4g}, {:.4g} "
              "based on a single non-zero bin and {} matches"
              .format(xp, yp, maxval))
        return xp, yp, maxval, maxval, zpmat, True

    (xp, yp), fit_status, fit_sl = _find_peak(zpmat, peak_fit_box=5,
                                              mask=zpmat > 0)

    if fit_status.startswith('ERROR'):
        print("WARNING: No valid shift found within a search radius of {:g} "
              "pixels.".format(searchrad))
        maxval = int(np.ceil(zpmat.max()))
        return 0.0, 0.0, maxval, maxval, zpmat, False

    xp -= searchrad
    yp -= searchrad

    if fit_status == 'WARNING:EDGE':
        print(
            "WARNING: Found peak in the 2D histogram lies at the edge of "
            "the histogram. Try increasing 'searchrad' for improved results."
        )

    # Attempt to estimate "significance of detection":
    maxval = zpmat.max()
    zpmat_mask = (zpmat > 0) & (zpmat < maxval)

    if np.any(zpmat_mask):
        bkg = zpmat[zpmat_mask].mean()
        sig = maxval / np.sqrt(bkg)

    flux = int(zpmat[fit_sl].sum())
    print("Found initial X and Y shifts of {:.4g}, {:.4g} "
          "with significance of {:.4g} and {:d} matches"
          .format(xp, yp, sig, flux))

    return xp, yp, int(np.ceil(maxval)), flux, zpmat, True
Create a 2D matrix-histogram which contains the delta between each
        XY position and each UV position. Then estimate initial offset
        between catalogs.