def get_external_catalog(phot, filter_file=
    '/usr/local/share/eazy-photoz/filters/FILTER.RES.latest', ZP=23.9,
    sys_err=0.3, verbose=True, external_limits=3, timeout=300):
    """
    Fetch photometry from vizier
    """
    import numpy as np
    import astropy.units as u
    try:
        from .. import utils
    except:
        from grizli import utils
    from eazy.filters import FilterFile
    res = FilterFile(filter_file)
    vizier_catalog = list(utils.VIZIER_VEGA.keys())
    ra = np.median(phot['ra'])
    dec = np.median(phot['dec'])
    dr = np.sqrt((phot['ra'] - ra) ** 2 * np.cos(phot['dec'] / 180 * np.pi) **
        2 + (phot['dec'] - dec) ** 2) * 60
    radius = 1.5 * dr.max()
    tabs = utils.get_Vizier_photometry(ra, dec, templates=None, radius=
        radius * 60, vizier_catalog=vizier_catalog, filter_file=filter_file,
        MW_EBV=0, convert_vega=False, raw_query=True, verbose=True, timeout
        =timeout)
    extern_phot = utils.GTable()
    N = len(phot)
    for t_i in tabs:
        if 'RAJ2000' in t_i.colnames:
            other_radec = 'RAJ2000', 'DEJ2000'
        elif 'RA_ICRS' in t_i.colnames:
            other_radec = 'RA_ICRS', 'DE_ICRS'
        else:
            other_radec = 'ra', 'dec'
        idx, dr = phot.match_to_catalog_sky(t_i, other_radec=other_radec)
        if (dr < 2 * u.arcsec).sum() == 0:
            continue
        tab = t_i[dr < 2 * u.arcsec]
        idx = idx[dr < 2 * u.arcsec]
        err_scale = 1
        if tab.meta['name'] == utils.UKIDSS_LAS_VIZIER:
            flux_scale = 1.33
        else:
            flux_scale = 1.0
        convert_vega = utils.VIZIER_VEGA[tab.meta['name']]
        bands = utils.VIZIER_BANDS[tab.meta['name']]
        to_flam = 10 ** (-0.4 * 48.6) * 3e+18
        for ib, b in enumerate(bands):
            f_number = res.search(b, verbose=False)[0] + 1
            filt = res.filters[f_number - 1]
            if convert_vega:
                to_ab = filt.ABVega()
            else:
                to_ab = 0.0
            fcol, ecol = bands[b]
            fnu = 10 ** (-0.4 * (tab[fcol] + to_ab - ZP))
            efnu = tab[ecol] * np.log(10) / 2.5 * fnu * err_scale
            efnu = np.sqrt(efnu ** 2 + (sys_err * fnu) ** 2)
            fnu.fill_value = -99
            efnu.fill_value = -99
            comment = 'Filter {0} from {1} (N={2})'.format(bands[b][0], t_i
                .meta['name'], len(idx))
            if verbose:
                print(comment)
            if ((~efnu.mask).sum() > 4) & (external_limits > 0):
                fill = np.percentile(efnu.data[~efnu.mask], [90])
                efill = external_limits * fill
            else:
                fill = -99
                efill = -99
            extern_phot.meta['F{0}'.format(f_number)] = b, comment
            extern_phot['F{0}'.format(f_number)] = fill * np.ones(N)
            extern_phot['F{0}'.format(f_number)][idx] = fnu.filled()
            extern_phot['E{0}'.format(f_number)] = efill * np.ones(N)
            extern_phot['E{0}'.format(f_number)][idx] = efnu.filled()
    return extern_phot
