def interpolate_to_points(da, ixds):
    """Interpolate data array to points

    Parameters
    ----------
    da : DataArray
        data to interpolate
    ixds : Dataset
        dataset with lon, lat, i, j coordinates
        to interpolate at

    Returns
    -------
    DataArray
        interpolated data array
        with points dimension
        and passed-through i, j indices
    """
    ida = da.interp(coords=ixds[['lon', 'lat']], method='linear')
    # preserve i,j indices
    for name in ['j', 'i']:
        ida.coords[name] = ixds[name]
    return ida