panoptes.pipeline.utils.sources¶

def read_catalog(catalog_filename) -> pandas.DataFrame:
    """Read a catalog file, choosing the reader from its suffix.

    The parsed catalog is cached (:py:data:`CATALOG_CACHE_SIZE`), keyed on the
    file's path, modification time and size, so re-reading the same unchanged
    file is free and a rewritten one is picked up. The cache lives in the
    calling process: a pool of frame workers reads the catalog once each, not
    once per frame. Call ``read_catalog.cache_clear()`` to drop it.

    What comes back is a shallow copy, so a caller that adds or replaces a
    column cannot corrupt the cached catalog for everything after it. That is
    copy-on-write doing the work, and it costs nothing until something is
    actually written -- which is why `pyproject.toml` floors pandas at 3, where
    copy-on-write is unconditional. On 2.x it can be off, and a caller writing
    through `.loc` would reach the cached frame through the shared blocks.

    `picid` comes back as a `category`. It identifies a star rather than
    measuring anything, and nothing downstream should be doing arithmetic on it.
    One thing to know when consuming it: `groupby` over a categorical iterates
    every category by default, so pass `observed=True` when grouping on it.

    Args:
        catalog_filename (str|Path): Path to a `.parquet`/`.pq`, `.ecsv` or
            `.csv`/`.tsv` file, optionally compressed (`.ecsv.gz`).

    Returns:
        `pandas.DataFrame`: The catalog, with `picid` as a categorical.

    Raises:
        FileNotFoundError: If the path does not exist.
        ValueError: If the suffix is not a recognized format, a required column
            is missing, or `picid` cannot be read as an integer identifier.

    """
    path = Path(catalog_filename)
    if not path.exists():
        raise FileNotFoundError(f"Catalog file does not exist: {path}")

    # `resolve` so two spellings of one file -- relative and absolute, or
    # through a symlink -- are one cache entry rather than two copies.
    path = path.resolve()
    stat = path.stat()

    return _parse_catalog(path, stat.st_mtime_ns, stat.st_size).copy(deep=False)

def get_stars_from_wcs(
    wcs0: WCS, round_to: int = 0, pad: float = 1.0, pad_size=(20, 10), **kwargs: Any
) -> pandas.DataFrame:
    """Lookup star information from WCS footprint.

    Generates the correct layout for an SQL `POLYGON` that can be passed to
    :py:func:`get_stars`.

    Args:
        wcs0 (astropy.wcs.WCS): A valid (i.e. `wcs.is_celestial`) World Coordinate System object.
        round_to (int): Round the limits to this decimal place, default 0. Keeps the
            requested bounds stable between frames of the same sequence.
        pad (float): The amount of padding in degrees to add to each of the RA and Dec
            limits, default 0.5 [degrees].
        **kwargs (Any): Optional keywords to pass to :py:func:`get_stars`.

    """
    wcs_footprint = wcs0.calc_footprint()
    print(f"Looking up catalog stars for WCS: {wcs_footprint}")

    ra_max, dec_max = (wcs0.wcs.crval + np.array(pad_size)).round(round_to)
    ra_min, dec_min = (wcs0.wcs.crval - np.array(pad_size)).round(round_to)

    limits = dict(ra_max=ra_max % 360, ra_min=ra_min % 360, dec_max=dec_max, dec_min=dec_min)

    print(f"Searching square shape with {round_to=} and {pad=}: {limits!r}")
    catalog_stars = get_stars(shape=limits, **kwargs)

    return catalog_stars

def get_stars(shape=None, vmag_min=7, vmag_max=14, catalog_filename=None, **kwargs: Any) -> pandas.DataFrame:
    """Look up star information from a local copy of the PANOPTES Input Catalog.

    `picid` is exactly the Gaia DR3 `source_id`, so the catalog is a Gaia cone
    search with its columns renamed and there is no crossmatch step -- see
    `scripts/fetch_catalog.py`, which builds one. It is read from a local file --
    parquet, ECSV or CSV, see :py:func:`read_catalog` -- and getting that file
    onto disk is a separate fetch step, not something this function does. There
    is no network lookup; see improvement plan 4.3.

    The file is expected to carry the mapped column names
    (:py:data:`REQUIRED_CATALOG_COLUMNS`) rather than the raw upstream ones.

    Note:

        In the upstream catalog the GAIA `bp` and `rp` magnitude and error
        columns are switched. Whatever produces the catalog is responsible for
        correcting that; this function does not second-guess the file it is
        given.

    Args:
        shape (dict|None): A dictionary containing the keys `ra_min`, `ra_max`,
            `dec_min`, `dec_max`, in degrees. If None, no positional filtering.
        vmag_min (float, optional): Minimum Vmag to include, inclusive.
        vmag_max (float, optional): Maximum Vmag to include, exclusive.
        catalog_filename (str|Path): Path to the catalog file: parquet, ECSV
            or CSV. Required; there is no default.
        **kwargs (Any): Ignored, for call-site compatibility.

    Returns:
        `pandas.DataFrame`: The catalog entries inside the requested bounds.

    Raises:
        ValueError: If no catalog path is given, or the file is missing columns.
        FileNotFoundError: If the catalog path does not exist.

    """
    if catalog_filename is None:
        raise ValueError(
            "A local catalog is required. Pass catalog_filename, or set "
            "`params.catalog.catalog_filename` in the pipeline settings. "
            "There is no network catalog lookup."
        )

    catalog_stars = read_catalog(catalog_filename)

    # Vmag range is [vmag_min, vmag_max), as documented.
    selected = catalog_stars.catalog_vmag.between(vmag_min, vmag_max, inclusive="left")

    if shape is not None:
        selected &= catalog_stars.catalog_dec.between(shape["dec_min"], shape["dec_max"])

        # Right Ascension wraps from 360 to 0, so the box can straddle the origin.
        ra = catalog_stars.catalog_ra
        if shape["ra_max"] < shape["ra_min"]:
            selected &= (ra >= shape["ra_min"]) | (ra <= shape["ra_max"])
        else:
            selected &= (ra >= shape["ra_min"]) & (ra <= shape["ra_max"])

    results = catalog_stars[selected].reset_index(drop=True)

    # Filtering a categorical keeps every category, so a field cut from an
    # all-sky catalog would otherwise carry millions of ids that are not in it.
    results["picid"] = results["picid"].cat.remove_unused_categories()

    print(f"Found {len(results)} in Vmag=[{vmag_min}, {vmag_max}) and bounds=[{shape}]")

    return results

def get_catalog_match(
    point_sources,
    wcs=None,
    catalog_stars=None,
    max_separation_arcsec=None,
    ra_column="measured_ra",
    dec_column="measured_dec",
    **kwargs: Any,
) -> pandas.DataFrame:
    """Match the point source positions to the catalog.

    `picid` in the catalog is exactly the Gaia DR3 `source_id`, so the catalog is
    a Gaia cone search with its columns renamed and there is no crossmatch step
    -- see `scripts/fetch_catalog.py`, which builds one.

    The catalog is read from a local file. This function matches the `ra_column`
    and `dec_column` positions to the `catalog_ra` and `catalog_dec` columns of
    the catalog. When `catalog_stars` is not supplied the lookup is done via
    :py:func:`get_stars_from_wcs`.

    The matched catalog row is joined onto each source and the result returned.
    The columns added are whatever the catalog file carries
    (:py:data:`REQUIRED_CATALOG_COLUMNS`, plus :py:data:`GAIA_CATALOG_COLUMNS`
    when the file has them), together with:

        * catalog_sep -- separation from the matched catalog position, in arcsec
        * catalog_wcs_x, catalog_wcs_y -- the catalog position in pixels
        * catalog_wcs_x_int, catalog_wcs_y_int -- the same, truncated to int

    Note:

        Every source is matched to its *nearest* catalog star, so a detection
        with no real counterpart still gets a row -- `max_separation_arcsec` is
        what removes those. Catalog entries with no detection are not returned,
        so the result has at most one row per source in `point_sources`.

    If a `max_separation_arcsec` is given then results will be filtered if their
    separation from the catalog was larger than the number given. Typical values
    would be in the range of 20-30 arcsecs, which corresponds to 2-3 pixels.

    Args:
        point_sources (`pandas.DataFrame`): The DataFrame containing point sources
            to be matched. This usually comes from the output of
            :py:func:`.images.detect_sources` but could be done manually.
        wcs (`astropy.wcs.WCS`): Required, despite the `None` default. Supplying
            `catalog_stars` skips the lookup but not the pixel positions, which
            :py:func:`get_xy_positions` computes from this WCS for every call,
            so omitting it raises `AttributeError` after the matching work is done.
        catalog_stars (`pandas.DataFrame`, optional): If provided, the catalog match
            will be performed against this set of stars rather than performing a lookup.
        max_separation_arcsec (float|None, optional): If not None, sources more
            than this many arcsecs from catalog will be filtered.
        ra_column (str): The column name to use for the RA coordinates, default `measured_ra`.
        dec_column (str): The column name to use for the Dec coordinates, default
            `measured_dec`. Both defaults are stale -- :py:func:`.images.detect_sources`
            names its columns `photutils_sky_centroid_ra` and `..._dec`, which is
            what :py:func:`.images.match_sources` passes.
        **kwargs (Any): Extra options are passed to `get_stars_from_wcs`, which
            passes them to `get_stars`. Only used when `catalog_stars` is None.

    Returns:
        `pandas.DataFrame`: A dataframe with the catalog information added to
            the sources.

    """
    assert point_sources is not None

    if catalog_stars is None:
        print(f"Looking up stars for wcs={wcs.wcs.crval}")
        # Lookup stars in catalog
        catalog_stars = get_stars_from_wcs(wcs, **kwargs)

    if catalog_stars is None:
        print("No catalog matches, returning table without ids")
        return point_sources

    # Get coords for catalog stars
    catalog_coords = SkyCoord(
        ra=catalog_stars["catalog_ra"].values * u.deg,
        dec=catalog_stars["catalog_dec"].values * u.deg,
        frame="icrs",
    )

    # Get coords from detected point sources
    stars_coords = SkyCoord(
        ra=point_sources[ra_column].values * u.deg,
        dec=point_sources[dec_column].values * u.deg,
        frame="icrs",
    )

    # Do catalog matching
    print(f"Matching {len(catalog_coords)} catalog stars to {len(stars_coords)} detected stars")
    idx, d2d, d3d = stars_coords.match_to_catalog_sky(catalog_coords)
    print(f"Got {len(idx)} matched sources (includes duplicates) for wcs={wcs.wcs.crval}")

    catalog_matches = catalog_stars.iloc[idx].copy()
    catalog_matches["catalog_sep"] = d2d.to_value(u.arcsec)

    # Get the XY positions
    catalog_matches = get_xy_positions(wcs, catalog_matches)

    # Add the matches and their separation.
    matched_sources = point_sources.reset_index(drop=True).join(catalog_matches.reset_index(drop=True))

    # All point sources so far are matched.
    # matched_sources['status'] = 'matched'
    # matched_sources.status = matched_sources.status.astype('category')

    # Sources that didn't match.
    #     if return_unmatched:
    #         print(f'Adding unmatched sources to table for wcs={wcs.wcs.crval!r}')
    #         unmatched = catalog_stars.iloc[catalog_stars.index.difference(idx)].copy()

    #         unmatched['status'] = 'unmatched'
    #         point_sources = point_sources.append(unmatched)

    # Reorder columns so id cols are first then alpha.
    # new_column_order = sorted(list(matched_sources.columns))
    # id_cols = ['picid', 'unit_id', 'camera_id', 'time', 'gaia', 'twomass', 'status']
    # for i, col in enumerate(id_cols):
    #     new_column_order.remove(col)
    #     new_column_order.insert(i, col)
    # matched_sources = matched_sources.reindex(columns=new_column_order)

    print(f"Point sources: {len(matched_sources)} for wcs={wcs.wcs.crval!r}")

    # Remove catalog matches that are too far away.
    if max_separation_arcsec is not None:
        print(f"Removing matches > {max_separation_arcsec} arcsec from catalog.")
        matched_sources = matched_sources.query("catalog_sep <= @max_separation_arcsec")

    print(f"Returning matched sources: {len(matched_sources)} for wcs={wcs.wcs.crval!r}")
    return matched_sources