inputs:
    ------
    az0, el0: N-D array of azimuth, elevation. May be masked arrays
    az, el: 1-D vectors of azimuth, elevation points from other camera to find closest angle for joint FOV.

    output:
    row, col:  index of camera 0 closest to camera 1 FOV for each unmasked pixel

    I think with some minor tweaks this could be numba.jit if too slow.

def _findindex(az0, el0, az, el):
    """
    inputs:
    ------
    az0, el0: N-D array of azimuth, elevation. May be masked arrays
    az, el: 1-D vectors of azimuth, elevation points from other camera to find closest angle for joint FOV.

    output:
    row, col:  index of camera 0 closest to camera 1 FOV for each unmasked pixel

    I think with some minor tweaks this could be numba.jit if too slow.
    """

    assert az0.size == el0.size  # just for clarity
    assert az.ndim == el.ndim == 1, 'expect vector of test points'
    ic = np.empty(az.size, dtype=int)

    for i, (a, e) in enumerate(zip(az, el)):
        # we do this point by point because we need to know the closest pixel for each point
        # errang = haver.anglesep(az,el, apt,ept, deg=False)
        ic[i] = haver.anglesep_meeus(az0, el0, a, e, deg=False).argmin()

    """
    THIS UNRAVEL_INDEX MUST BE ORDER = 'C'
    """
    r, c = np.unravel_index(ic, az0.shape, order='C')

    mask = (c == 0) | (c == az0.shape[1] - 1) | (r == 0) | (r == az0.shape[0] - 1)

    r = np.ma.masked_where(mask, r)
    c = np.ma.masked_where(mask, c)

    return r, c