Proper motion

docs/changes/2625.features.rst

@@ -0,0 +1,4 @@
+Add ''get_bright_stars_with_motion'' to ctapipe.utils.astro
+this function allows to look up the position of stars
+around a specified location including proper motion.
+different catalogues can be used.

pyproject.toml

@@ -60,6 +60,7 @@ tests = [
     "pytest-xdist",
     "pytest_astropy_header",
     "tomli",
+    "astroquery",
 ]
 
 docs = [

src/ctapipe/utils/astro.py

@@ -4,63 +4,231 @@
 not provided by external packages and/or to satisfy particular needs of
 usage within ctapipe.
 """
-from astropy import units as u
+
+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"]
 
 
-def get_bright_stars(pointing=None, radius=None, magnitude_cut=None):
+CACHE_FILE = Path("~/.psf_stars.ecsv").expanduser()
+
+
+class StarCatalogues(Enum):
+    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):
+    """
+    utility function to cut stars based on magnitude and/or location
+
+    Parameters
+    ----------
+    stars: astropy table
+        Table of stars, including magnitude and coordinates
+    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
+    Vmag_cut: float
+        Return only stars above a given apparent magnitude. Default: None (all entries)
+    Bmag_cut: float
+        Return only stars above a given absolute magnitude. Default: None (all entries)
+
+    Returns
+    -------
+    Astropy table:
+       List of all stars after cuts with same keys as the input table stars
     """
-    Get an astropy table of bright stars.
 
-    This function is using the Yale bright star catalog, available through ctapipe
-    data downloads.
+    if Bmag_cut is not None and "Bmag" in stars.keys():
+        stars = stars[stars["Bmag"] < Bmag_cut]
+    elif Bmag_cut is not None and "BTmag" in stars.keys():
+        stars = stars[stars["BTmag"] < Bmag_cut]
+    if Vmag_cut is not None and "Vmag" in stars.keys():
+        stars = stars[stars["Vmag"] < Vmag_cut]
+    elif Vmag_cut is not None and "VTmag" in stars.keys():
+        stars = stars[stars["VTmag"] < Vmag_cut]
 
-    The included Yale bright star catalog contains all 9096 stars, excluding the
-    Nova objects present in the original catalog and is complete down to magnitude
-    ~6.5, while the faintest included star has mag=7.96. :cite:p:`bright-star-catalog`
+    if radius is not None:
+        if pointing is None:
+            raise ValueError(
+                "Sky pointing, pointing=SkyCoord(), must be "
+                "provided if radius is given."
+            )
+        separations = stars["ra_dec"].separation(pointing)
+        stars["separation"] = separations
+        stars = stars[separations < radius]
+
+    return stars
+
+
+def get_star_catalog(catalog, filename):
+    """
+    Utility function to download a star catalog for the get_bright_stars function
+
+    Parameters
+    ----------
+    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
+    """
+    from astroquery.vizier import Vizier
+
+    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" 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)
+
+
+def get_bright_stars(
+    pointing=None, radius=None, Bmag_cut=None, Vmag_cut=None, catalog="Yale"
+):  # max_magnitude):
+    """
+    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
 
     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 magnitude. Default: None (all entries)
+    Vmag_cut: float
+        Return only stars above a given apparent magnitude. Default: None (all entries)
+    Bmag_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
     """
-    from ctapipe.utils import get_table_dataset
 
-    catalog = get_table_dataset("yale_bright_star_catalog5", role="bright star catalog")
+    stars = None
 
-    starpositions = SkyCoord(
-        ra=Angle(catalog["RAJ2000"], unit=u.deg),
-        dec=Angle(catalog["DEJ2000"], unit=u.deg),
-        frame="icrs",
-        copy=False,
-    )
-    catalog["ra_dec"] = starpositions
+    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 magnitude_cut is not None:
-        catalog = catalog[catalog["Vmag"] < magnitude_cut]
+    if stars is None:
+        from astroquery.vizier import Vizier
 
-    if radius is not None:
-        if pointing is None:
-            raise ValueError(
-                "Sky pointing, pointing=SkyCoord(), must be "
-                "provided if radius is given."
-            )
-        separations = catalog["ra_dec"].separation(pointing)
-        catalog["separation"] = separations
-        catalog = catalog[separations < radius]
+        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
 
-    catalog.remove_columns(["RAJ2000", "DEJ2000"])
+        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"),
+        )
+
+    stars = select_stars(
+        stars, pointing=pointing, radius=radius, Bmag_cut=Bmag_cut, Vmag_cut=Vmag_cut
+    )
 
-    return catalog
+    return stars

src/ctapipe/utils/tests/test_astro.py

@@ -3,19 +3,46 @@
 This module contains the utils.astro unit tests
 """
 from astropy import units as u
-from astropy.coordinates import SkyCoord
+from astropy.coordinates import Angle, SkyCoord
+from astropy.time import Time
 
 from ..astro import get_bright_stars
 
 
 def test_get_bright_stars():
     """
-    unit test for utils.astro.get_bright_stars(). Tests that only Zeta Tau is
-    returned close to the Crab Nebula as object brighter than mag=3.5.
+    unit test for utils.astro.get_bright_stars_with_motion().
     """
-    pointing = SkyCoord(ra=83.275 * u.deg, dec=21.791 * u.deg, frame="icrs")
+    # TODO add tests for all catalogues, specifically by trying to find some particular bright star, also test that motion is properly included
+    eta_tau = SkyCoord(
+        ra=Angle("03 47 29.1", unit=u.deg),
+        dec=Angle("+24 06 18", unit=u.deg),
+        frame="icrs",
+    )
+    obstime = Time("2020-01-01")
 
-    table = get_bright_stars(pointing, radius=2.0 * u.deg, magnitude_cut=3.5)
+    # lets find 25 Eta Tau
+
+    table = get_bright_stars(pointing=eta_tau, radius=1.0 * u.deg, Vmag_cut=3.5)
+
+    assert len(table) == 1  # this looks if 25 Eita Tau was found
+    # check that the star moves
+    assert table[0]["ra_dec"].apply_space_motion(obstime).ra != eta_tau.ra
+    assert table[0]["ra_dec"].apply_space_motion(obstime).dec != eta_tau.dec
+
+    # now test the other catalog. First get object 766 from the catalog
+
+    HIP_star = SkyCoord(
+        ra=Angle("00 08 23.2585712", unit=u.deg),
+        dec=Angle("+29 05 25.555166", unit=u.deg),
+        frame="icrs",
+    )
+
+    table = get_bright_stars(
+        pointing=HIP_star, radius=1.0 * u.deg, Vmag_cut=3.5, catalog="Hippoarcos"
+    )
 
     assert len(table) == 1
-    assert table[0]["Name"] == "123Zet Tau"
+    # now check that stars move
+    assert table[0]["ra_dec"].apply_space_motion(obstime).ra != eta_tau.ra
+    assert table[0]["ra_dec"].apply_space_motion(obstime).dec != eta_tau.dec