Added hipparcos star catalog

.gitignore

@@ -12,6 +12,7 @@ __pycache__
 ruststartracker/libruststartracker.*.so
 ruststartracker/libruststartracker.so
 ruststartracker/gaia*.csv
+ruststartracker/hipparcos*.csv
 
 .coverage
 coverage.xml

Cargo.toml

@@ -28,3 +28,4 @@ rand_distr = "0.5.1"
 default = []
 improc = ["opencv"]
 gaia = []
+hipparcos = []

README.md

@@ -91,3 +91,14 @@ Gaia DR3 data is © European Space Agency and is released under the [**Creative
 
 > Gaia Collaboration, Vallenari et al. (2022), *A\&A* **674**, A1.
 > [DOI: 10.1051/0004-6361/202243940](https://doi.org/10.1051/0004-6361/202243940)
+
+### Hipparcos and Tycho Data
+
+This project includes data from the European Space Agency (ESA) mission **Hipparcos**.
+
+The Hipparcos and Tycho Catalogues were processed by the Hipparcos and Tycho Data Analysis Consortium.
+
+The Hipparcos and Tycho Catalogues are © European Space Agency and are released under the [**Creative Commons Attribution 3.0 IGO (CC BY 3.0 IGO)**](https://creativecommons.org/licenses/by/3.0/igo/) license.
+
+> Perryman, M. A. C., et al. (1997), *Astronomy & Astrophysics* **323**, L49-L52.
+1997A&A...323L..49P

build_script.py

@@ -3,6 +3,8 @@
 This script is automatically run when the pyproject is installed.
 """
 
+import csv
+import io
 import pathlib
 import shutil
 import subprocess
@@ -16,7 +18,7 @@ def download_gaia_data(output_file: pathlib.Path) -> None:
     """
     from astroquery.gaia import Gaia
 
-    Gaia.ROW_LIMIT = 10000
+    Gaia.ROW_LIMIT = -1
 
     query = """
     SELECT source_id, ra, dec, parallax, pmra, pmdec, phot_g_mean_mag
@@ -42,6 +44,61 @@ def download_gaia_data(output_file: pathlib.Path) -> None:
     print(f"Saved to: {output_file}")
 
 
+def download_hipparcos_data(output_file: pathlib.Path) -> None:
+    """Download Hipparcos data and save it to the specified output file.
+
+    Args:
+        output_file: Path to save the downloaded Hipparcos data (csv).
+    """
+    from astroquery.vizier import Vizier
+
+    Vizier.ROW_LIMIT = -1
+
+    print("Executing query on Hipparcos database...")
+    query = Vizier.query_constraints(
+        catalog=["I/311/hip2"],
+        Hpmag="<7",
+    )
+
+    if not query:
+        print("Job did not finish successfully.")
+        return
+
+    # Initialize a string buffer to write the CSV data to
+    csv_output = io.StringIO()
+    writer = csv.writer(csv_output)
+
+    # Write the header row
+    # Use column names that are common in Gaia and Hipparcos
+    header = ["source_id", "ra", "dec", "parallax", "pmra", "pmdec", "phot_g_mean_mag"]
+    writer.writerow(header)
+
+    # Iterate over the rows of the Astropy table
+    for row in query[0]:
+        # Extract the required data points
+        source_id = row["HIP"]
+        ra = row["RArad"]
+        dec = row["DErad"]
+        parallax = row["Plx"]
+        pmra = row["pmRA"]
+        pmdec = row["pmDE"]
+        # Use Vmag as the proxy for phot_g_mean_mag
+        phot_g_mean_mag = row["Hpmag"]
+
+        # Create a list for the new row and write it to the CSV writer
+        new_row = [source_id, ra, dec, parallax, pmra, pmdec, phot_g_mean_mag]
+        writer.writerow(new_row)
+
+    # Get the complete CSV string from the buffer
+    csv_string = csv_output.getvalue()
+
+    # You can save this string to a file
+    with output_file.open("w") as f:
+        f.write(csv_string)
+
+    print(f"Saved to: {output_file}")
+
+
 def build_script() -> None:
     """Build rust backend and move shared library to correct folder."""
     cwd = pathlib.Path(__file__).parent.expanduser().absolute()
@@ -50,8 +107,12 @@ def build_script() -> None:
     if not gaia_file.exists():
         download_gaia_data(gaia_file)
 
+    hipparcos_file = cwd / "ruststartracker/hipparcos_data_j1991.25.csv"
+    if not hipparcos_file.exists():
+        download_hipparcos_data(hipparcos_file)
+
     subprocess.check_call(  # noqa: S603
-        ["cargo", "build", "--release", "--features", "improc,gaia"],  # noqa: S607
+        ["cargo", "build", "--release", "--features", "improc,gaia,hipparcos"],  # noqa: S607
         cwd=cwd,
         stdout=None,
         stderr=None,

pyproject.toml

@@ -7,6 +7,7 @@ readme = "README.md"
 license = "MIT"
 include = [
     { path = "ruststartracker/gaia_data_j2016.csv", format = ["sdist", "wheel"] },
+    { path = "ruststartracker/hipparcos_data_j1991.25.csv", format = ["sdist", "wheel"] },
     { path = "ruststartracker/libruststartracker.so", format = ["wheel"] },
 ]
 
@@ -33,6 +34,7 @@ scipy-stubs = "^1.14.1.5"
 ruff = "^0.9.3"
 bump-my-version = "^1.2.0"
 astroquery = "^0.4.10"
+typing-extensions = "^4.12.2"
 
 [tool.poetry.build]
 script = "build_script.py"

ruststartracker/catalog.py

@@ -7,14 +7,17 @@
 
 import numpy as np
 import numpy.typing as npt
+from typing_extensions import Self
 
 AU: float = 149597870.693
 """Astronomical unit."""
 
-INTERNAL_CATALOG_FILE = (
-    pathlib.Path(__file__).parent.expanduser().absolute() / "gaia_data_j2016.csv"
+GAIA_CATALOG_FILE = pathlib.Path(__file__).parent.expanduser().absolute() / "gaia_data_j2016.csv"
+"""Location of the internal Gaia star catalog file."""
+HIPPARCOS_CATALOG_FILE = (
+    pathlib.Path(__file__).parent.expanduser().absolute() / "hipparcos_data_j1991.25.csv"
 )
-"""Location of the internal star catalog file."""
+"""Location of the internal Hipparcos star catalog file."""
 
 
 def time_to_epoch(t: datetime.datetime) -> float:
@@ -44,24 +47,41 @@ class StarCatalog:
     """Proper motion of the declination in mad."""
     magnitude: npt.NDArray[np.float32]
     """Magnitude values."""
-    epoch: float = 2016.0
+    epoch: float
     """Epoch of the catalog in years."""
 
+    @classmethod
+    def from_gaia(cls, *, max_magnitude: float = 6.0) -> Self:
+        """Read internally provided Gaia catalog file.
+
+        Args:
+            max_magnitude: Maximum magnitude to include.
+        """
+        return cls(GAIA_CATALOG_FILE, epoch=2016.0, max_magnitude=max_magnitude)
+
+    @classmethod
+    def from_hipparcos(cls, *, max_magnitude: float = 6.0) -> Self:
+        """Read internally provided Hipparcos catalog file.
+
+        Args:
+            max_magnitude: Maximum magnitude to include.
+        """
+        return cls(HIPPARCOS_CATALOG_FILE, epoch=1991.25, max_magnitude=max_magnitude)
+
     def __init__(
-        self, filename: pathlib.Path | str | None = None, max_magnitude: float = 6.0
+        self, filename: pathlib.Path | str, *, epoch: float, max_magnitude: float = 6.0
     ) -> None:
         """Read catalog from file.
 
         Args:
             filename: Star catalog filename.
+            epoch: Epoch of the catalog in years.
             max_magnitude: Maximum magnitude to include.
         """
-        # Use locally included star catalog is no file is given
-        filename = INTERNAL_CATALOG_FILE if filename is None else pathlib.Path(filename)
-
+        self.epoch = epoch
         keep_columns = ("ra", "dec", "parallax", "pmra", "pmdec", "phot_g_mean_mag")
 
-        with filename.open("r") as f:
+        with pathlib.Path(filename).open("r") as f:
             it = csv.reader(f, delimiter=",", strict=True)
 
             # Get columns

ruststartracker/libruststartracker.pyi

@@ -84,6 +84,8 @@ def extract_observations(
 class StarCatalog:
     @classmethod
     def from_gaia(cls, *, max_magnitude: float | None) -> Self: ...
+    @classmethod
+    def from_hipparcos(cls, *, max_magnitude: float | None) -> Self: ...
     def normalized_positions(
         self, *, epoch: float | None, observer_position: np.ndarray | None
     ) -> npt.NDArray[np.float32]: ...

ruststartracker/test_catalog.py

@@ -38,16 +38,16 @@ def test_time_to_epoch_astropy(iso_date: str):
 
 
 def test_extract_observations():
-    positions = ruststartracker.catalog.StarCatalog().normalized_positions()
+    positions = ruststartracker.catalog.StarCatalog.from_gaia().normalized_positions()
 
     assert positions.ndim == 2
     assert positions.shape[1] == 3
     np.testing.assert_allclose(np.linalg.norm(positions, axis=-1), 1.0, rtol=1e-5)
 
 
-def test_python_rust():
+def test_gaia_python_rust():
     t0 = time.monotonic()
-    positions = ruststartracker.catalog.StarCatalog().normalized_positions(epoch=2025.0)
+    positions = ruststartracker.catalog.StarCatalog.from_gaia().normalized_positions(epoch=2025.0)
     print(f"Python catalog took {time.monotonic() - t0:.3f} seconds")
     t0 = time.monotonic()
     positions2 = ruststartracker.libruststartracker.StarCatalog.from_gaia(
@@ -57,5 +57,45 @@ def test_python_rust():
     np.testing.assert_allclose(positions, positions2, rtol=1e-5, atol=1e-5)
 
 
+def test_hipparcos_python_rust():
+    t0 = time.monotonic()
+    positions = ruststartracker.catalog.StarCatalog.from_hipparcos().normalized_positions(
+        epoch=2025.0
+    )
+    print(f"Python catalog took {time.monotonic() - t0:.3f} seconds")
+    t0 = time.monotonic()
+    positions2 = ruststartracker.libruststartracker.StarCatalog.from_hipparcos(
+        max_magnitude=6.0
+    ).normalized_positions(epoch=2025.0, observer_position=None)
+    print(f"Rust catalog took {time.monotonic() - t0:.3f} seconds")
+    np.testing.assert_allclose(positions, positions2, rtol=1e-5, atol=1e-5)
+
+
+def test_compare_hipparcos_gaia():
+    def filt(cat: ruststartracker.catalog.StarCatalog, m: float = 0):
+        positions = cat.normalized_positions(epoch=2025.0)
+        mags = cat.magnitude
+        mask = (positions[..., 2] > 0.9) * (mags > 5 - m) * (mags < 8 - m)
+        return positions[mask]
+
+    positions_hipparcos = filt(ruststartracker.catalog.StarCatalog.from_hipparcos(max_magnitude=8))
+    positions_gaia = filt(ruststartracker.catalog.StarCatalog.from_gaia(max_magnitude=8), m=0.9)
+
+    dists = np.linalg.norm(
+        positions_hipparcos[np.newaxis, :, :] - positions_gaia[:, np.newaxis, :], axis=-1
+    )
+    matches = np.min(dists, axis=0) < 0.00005
+
+    assert np.mean(matches > 0.9)
+
+    if False:
+        import matplotlib.pyplot as plt
+
+        plt.plot(positions_gaia[..., 0], positions_gaia[..., 1], "+", label="gaia")
+        plt.plot(positions_hipparcos[..., 0], positions_hipparcos[..., 1], "x", label="hipparcos")
+        plt.legend()
+        plt.show()
+
+
 if __name__ == "__main__":
     pytest.main([__file__])

ruststartracker/test_integration.py

@@ -24,7 +24,7 @@ def prepare() -> tuple[ruststartracker.StarTracker, np.ndarray]:
 
     dist_coefs = np.array([-0.44120807, -0.15954202, 0.00767012, -0.00213292, -1.64788247])
 
-    catalog = ruststartracker.StarCatalog()
+    catalog = ruststartracker.StarCatalog.from_gaia()
     star_catalog_vecs = catalog.normalized_positions(epoch=2024)
     star_catalog_magnitudes = catalog.magnitude
 

ruststartracker/test_star.py

@@ -20,12 +20,10 @@ def test_extract_observations(impl: str):
     t0 = time.monotonic()
     if impl == "python":
         centers, intensities = ruststartracker.star._extract_observations(img, threshold=30)
-    elif impl == "rust":
+    else:
         centers, intensities = ruststartracker.libruststartracker.extract_observations(
             img, 30, 3, 300
         )
-    else:
-        raise AssertionError
     print(f"Extracting observations took {time.monotonic() - t0:.5f} seconds")
 
     assert isinstance(centers, np.ndarray)

src/lib.rs

@@ -254,6 +254,15 @@ impl StarCatalog {
         })
     }
 
+    #[cfg(feature = "hipparcos")]
+    #[classmethod]
+    fn from_hipparcos(_cls: &Bound<'_, PyType>, max_magnitude: Option<f64>) -> PyResult<Self> {
+        Ok(StarCatalog {
+            inner: starcat::StarCatalog::new_from_hipparcos(max_magnitude)
+                .map_err(PyRuntimeError::new_err)?,
+        })
+    }
+
     pub fn normalized_positions(
         &self,
         epoch: Option<f64>,

src/starcat.rs

@@ -8,11 +8,17 @@ use nalgebra::Vector3;
 use serde::Deserialize;
 
 const AU: f64 = 149_597_870.693;
+
 #[cfg(feature = "gaia")]
 const GAIA_EPOCH: f64 = 2016.0;
 #[cfg(feature = "gaia")]
 const GAIA_2016_CSV: &str = include_str!("../ruststartracker/gaia_data_j2016.csv");
 
+#[cfg(feature = "hipparcos")]
+const HIPPARCOS_EPOCH: f64 = 1991.25;
+#[cfg(feature = "hipparcos")]
+const HIPPARCOS_1991_25_CSV: &str = include_str!("../ruststartracker/hipparcos_data_j1991.25.csv");
+
 pub trait Cast<T> {
     fn cast(self) -> T;
 }
@@ -110,6 +116,11 @@ impl StarCatalog {
         StarCatalog::new_from_string(GAIA_2016_CSV, GAIA_EPOCH, max_magnitude)
     }
 
+    #[cfg(feature = "hipparcos")]
+    pub fn new_from_hipparcos(max_magnitude: Option<f64>) -> Result<Self, String> {
+        StarCatalog::new_from_string(HIPPARCOS_1991_25_CSV, HIPPARCOS_EPOCH, max_magnitude)
+    }
+
     fn new_from_buffer<T: std::io::Read>(
         reader: BufReader<T>,
         epoch: f64,
@@ -259,4 +270,17 @@ mod tests {
         let cat2 = StarCatalog::new_from_gaia(Some(5.5)).unwrap();
         assert!(cat1.stars.len() == cat2.stars.len());
     }
+
+    #[cfg(feature = "hipparcos")]
+    #[test]
+    fn test_hipparcos() {
+        let cat1 = StarCatalog::new_from_file(
+            "ruststartracker/hipparcos_data_j1991.25.csv",
+            GAIA_EPOCH,
+            Some(5.5),
+        )
+        .unwrap();
+        let cat2 = StarCatalog::new_from_hipparcos(Some(5.5)).unwrap();
+        assert!(cat1.stars.len() == cat2.stars.len());
+    }
 }