ij_1d = np.linalg.norm(delta, axis=2).argmin()
i, j = np.unravel_index(ij_1d, lon.shape)
return i, j

def loopy_app(lon, lat, pts):
    return [find_indices(lon, lat, pts[i, 0], pts[i, 1]) for i in range(pts
        .shape[0])]

def vectorized_app(lon, lat, points):
    lonlat = np.column_stack((lon.ravel(), lat.ravel()))
idx = scipy.spatial.distance.cdist(lonlat, points).argmin(0)
return np.column_stack(np.unravel_index(idx, lon.shape)).tolist()
Efficiently find indices of nearest points on non-rectangular 2D grid