From 89f4e376bde180ab09eae4c7865638317467331b Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Tue, 22 Oct 2024 14:36:47 +0200 Subject: [PATCH] Update based on latest comments --- docs/changes/2625.features.rst | 3 +- src/ctapipe/utils/astro.py | 188 +++++++++++++-------------------- 2 files changed, 74 insertions(+), 117 deletions(-) diff --git a/docs/changes/2625.features.rst b/docs/changes/2625.features.rst index 4ea9e5b020a..70644eed742 100644 --- a/docs/changes/2625.features.rst +++ b/docs/changes/2625.features.rst @@ -1,4 +1,5 @@ -Add ''get_bright_stars_with_motion'' to ctapipe.utils.astro +Modifies ''get_bright_stars'' in ctapipe.utils.astro to include proper motion +Add ''get_hipparcos_stars'' to use the Hippoarcos catalog to get star positions with proper motion this function allows to look up the position of stars around a specified location including proper motion. different catalogues can be used. diff --git a/src/ctapipe/utils/astro.py b/src/ctapipe/utils/astro.py index a3e94a313d4..8c259d74b9a 100644 --- a/src/ctapipe/utils/astro.py +++ b/src/ctapipe/utils/astro.py @@ -7,24 +7,19 @@ import logging from enum import Enum -from pathlib import Path import astropy.units as u from astropy.coordinates import Angle, SkyCoord -from astropy.table import Table from astropy.time import Time log = logging.getLogger("main") -__all__ = ["get_bright_stars"] - - -CACHE_FILE = Path("~/.psf_stars.ecsv").expanduser() +__all__ = ["get_bright_stars", "get_hipparcos_stars"] class StarCatalogues(Enum): - Yale = "V/50/catalog" # Yale bright star catalogue - Hippoarcos = "I/239/hip_main" # hipparcos catalogue + Yale = "V/50/catalog" #: Yale bright star catalogue + Hippoarcos = "I/239/hip_main" #: hipparcos catalogue def select_stars(stars, pointing=None, radius=None, Bmag_cut=None, Vmag_cut=None): @@ -72,7 +67,7 @@ def select_stars(stars, pointing=None, radius=None, Bmag_cut=None, Vmag_cut=None return stars -def get_star_catalog(catalog, filename): +def get_star_catalog(catalog): """ Utility function to download a star catalog for the get_bright_stars function @@ -80,15 +75,20 @@ def get_star_catalog(catalog, filename): ---------- catalog: string Name of the catalog to be used. Usable names are found in the Enum StarCatalogues. Default: Yale - filename: string - Name and location of the catalog file to be produced + + Returns + ---------- + Astropy table: + List of all stars after cuts with catalog numbers, magnitudes, + and coordinates as SkyCoord objects including proper motion """ from astroquery.vizier import Vizier vizier = Vizier( catalog=StarCatalogues[catalog].value, columns=[ - "recno", + "HIP", #: HIP is the Hippoarcos ID available for that catalog + "HR", #: HR is the Harvard Revised Number available for the Yale catalog "RAJ2000", "DEJ2000", "RAICRS", @@ -97,30 +97,61 @@ def get_star_catalog(catalog, filename): "pmDE", "Vmag", "Bmag", - "BTmag", - "VTmag", ], row_limit=1000000, ) stars = vizier.query_constraints(Vmag="0.0..10.0")[0] - if "Bmag" not in stars.keys(): - log.exception("The chosen catalog does not have Bmag data") - if "Vmag" not in stars.keys(): - log.exception("The chosen catalog does not have Vmag data") stars.meta["Catalog"] = StarCatalogues[catalog].value - stars.write(CACHE_FILE, overwrite=True) + return stars + + +def get_bright_stars(pointing=None, radius=None, magnitude_cut=None): + """ + Get an astropy table of bright stars from the Yale bright star catalog + Parameters + ---------- + pointing: astropy Skycoord + pointing direction in the sky (if none is given, full sky is returned) + radius: astropy angular units + Radius of the sky region around pointing position. Default: full sky + magnitude_cut: float + Return only stars above a given absolute magnitude. Default: None (all entries) -def get_bright_stars( - pointing=None, radius=None, Bmag_cut=None, Vmag_cut=None, catalog="Yale" -): # max_magnitude): + Returns + ------- + Astropy table: + List of all stars after cuts with catalog numbers, magnitudes, + and coordinates as SkyCoord objects including proper motion """ - Get an astropy table of bright stars from a VizieR catalog - You can browse the catalogs at https://vizier.cds.unistra.fr/viz-bin/VizieR + from ctapipe.utils import get_table_dataset + + stars = get_table_dataset("yale_bright_star_catalog5", role="bright star catalog") + + stars["ra_dec"] = SkyCoord( + ra=Angle(stars["RAJ2000"], unit=u.deg), + dec=Angle(stars["DEJ2000"], unit=u.deg), + pm_ra_cosdec=stars["pmRA"].quantity, # yes, pmRA is already pm_ra_cosdec + pm_dec=stars["pmDE"].quantity, + frame="icrs", + obstime=Time("J2000.0"), + ) + stars.remove_columns(["RAJ2000", "DEJ2000"]) + + stars = select_stars( + stars, pointing=pointing, radius=radius, Vmag_cut=magnitude_cut + ) + + return stars + + +def get_hipparcos_stars(pointing=None, radius=None, magnitude_cut=None): + """ + Get an astropy table of bright stars from the Hippoarcos star catalog Parameters ---------- @@ -128,107 +159,32 @@ def get_bright_stars( pointing direction in the sky (if none is given, full sky is returned) radius: astropy angular units Radius of the sky region around pointing position. Default: full sky - Vmag_cut: float - Return only stars above a given apparent magnitude. Default: None (all entries) - Bmag_cut: float + magnitude_cut: float Return only stars above a given absolute magnitude. Default: None (all entries) - catalog: string - Name of the catalog to be used. Usable names are found in the Enum StarCatalogues. Default: Yale Returns ------- Astropy table: - List of all stars after cuts with names, catalog numbers, magnitudes, - and coordinates + List of all stars after cuts with catalog numbers, magnitudes, + and coordinates as SkyCoord objects including proper motion """ - stars = None - - if CACHE_FILE.exists(): - log.info("Loading stars from cached table") - try: - stars = Table.read(CACHE_FILE) - if stars.meta["Catalog"] == StarCatalogues[catalog].value: - if Bmag_cut is not None and "Bmag_cut" in stars.meta.keys(): - if stars.meta["Bmag_cut"] >= Bmag_cut: - log.debug(f"Loaded table is valid for { Bmag_cut= }") - else: - log.debug( - "Loaded cache table has smaller magnitude_cut, reloading" - ) - stars = None - if Vmag_cut is not None and "Vmag_cut" in stars.meta.keys(): - if stars.meta["Vmag_cut"] >= Vmag_cut: - log.debug(f"Loaded table is valid for {Vmag_cut= }") - else: - log.debug( - "Loaded cache table has smaller magnitude_cut, reloading" - ) - stars = None - else: - stars = None - except Exception: - log.exception("Cache file exists but reading failed. Recreating") - - if stars is None: - from astroquery.vizier import Vizier - - log.info("Querying Vizier for bright stars catalog") - # query vizier for stars with 0 <= Vmag <= max_magnitude - vizier = Vizier( - catalog=StarCatalogues[catalog].value, - columns=[ - "recno", - "RAJ2000", - "DEJ2000", - "RAICRS", - "DEICRS", - "pmRA", - "pmDE", - "Vmag", - "Bmag", - "BTmag", - "VTmag", - ], - row_limit=1000000, - ) - - stars = vizier.query_constraints(Vmag="0.0..10.0")[0] - if "Bmag" in stars.keys(): - if Bmag_cut is not None: - stars.meta["Bmag_cut"] = Bmag_cut - elif Bmag_cut is not None: - log.exception("The chosen catalog does not have Bmag data") - if "Vmag" in stars.keys(): - if Vmag_cut is not None: - stars.meta["Vmag_cut"] = Vmag_cut - elif Vmag_cut is not None: - log.exception("The chosen catalog does not have Vmag data") - stars.meta["Catalog"] = StarCatalogues[catalog].value - - stars.write(CACHE_FILE, overwrite=True) - - if "RAJ2000" in stars.keys(): - stars["ra_dec"] = SkyCoord( - ra=Angle(stars["RAJ2000"], unit=u.deg), - dec=Angle(stars["DEJ2000"], unit=u.deg), - pm_ra_cosdec=stars["pmRA"].quantity, # yes, pmRA is already pm_ra_cosdec - pm_dec=stars["pmDE"].quantity, - frame="icrs", - obstime=Time("J2000.0"), - ) - elif "RAICRS" in stars.keys(): - stars["ra_dec"] = SkyCoord( - ra=Angle(stars["RAICRS"], unit=u.deg), - dec=Angle(stars["DEICRS"], unit=u.deg), - pm_ra_cosdec=stars["pmRA"].quantity, # yes, pmRA is already pm_ra_cosdec - pm_dec=stars["pmDE"].quantity, - frame="icrs", - obstime=Time("J1991.25"), - ) + from ctapipe.utils import get_table_dataset + + stars = get_table_dataset("hippoarcos_star_catalog5", role="bright star catalog") + + stars["ra_dec"] = SkyCoord( + ra=Angle(stars["RAICRS"], unit=u.deg), + dec=Angle(stars["DEICRS"], unit=u.deg), + pm_ra_cosdec=stars["pmRA"].quantity, # yes, pmRA is already pm_ra_cosdec + pm_dec=stars["pmDE"].quantity, + frame="icrs", + obstime=Time("J1991.25"), + ) + stars.remove_columns(["RAICRS", "DEICRS"]) stars = select_stars( - stars, pointing=pointing, radius=radius, Bmag_cut=Bmag_cut, Vmag_cut=Vmag_cut + stars, pointing=pointing, radius=radius, Vmag_cut=magnitude_cut ) return stars