Merge pull request #33 from GernotMaier/package-data

Add stars to package data

Author: GernotMaier

README.md

@@ -57,8 +57,7 @@ python v2dl5/scripts/generate_runlist.py \
 
 ```console
 python v2dl5/scripts/reflected_region_analysis.py \
-    --obs_table ../../../VTS/DL3/v490/dl3_pointlike_moderate2tel/obs-index.fits.gz \
-    --runlist my_output_dir/runlist.txt \
+    --run_list my_output_dir/runlist.txt \
     --config examples/reflected_region_analysis.yml \
     --output_dir my_output_dir
 ```
@@ -70,3 +69,24 @@ python v2dl5/scripts/plot_binary_light_curves.py \
     --instrument VERITAS \
     --configuration examples/binary_lightcurve_plotting.yml
 ```
+
+## Auxiliary data
+
+Auxiliary data is stored in v2dl5/data and available at run time. This includes:
+
+### Star catalogues
+
+Hippargos catalog for stars with magnitude < 9: [v2dl5/data/hip_mag9.fits.gz](v2dl5/data/hip_mag9.fits.gz).
+Star catalogs are listed in the configuration files as
+
+```yaml
+datasets:
+    exclusion_region:
+        on_radius: 0.5 deg
+        magnitude_B: 7
+        star_exclusion_radius: 0.30 deg
+        fov: 3.5 deg
+        star_file: hip_mag9.fits.gz
+```
+
+Star catalogs are expected to be in the [v2dl5/data](v2dl5/data) directory and of fits format.README.md

examples/reflected_region.yml

@@ -26,7 +26,7 @@ datasets:
         magnitude_B: 7
         star_exclusion_radius: 0.30 deg
         fov: 3.5 deg
-        star_file: ./data/hip_mag9.fits.gz
+        star_file: hip_mag9.fits.gz
     containment_correction: false
     safe_mask:
         methods: ['aeff-max']

pyproject.toml

@@ -4,6 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [tool.setuptools]
 packages=["v2dl5"]
+include-package-data = true
 
 [tool.setuptools_scm]
 write_to = "v2dl5/_version.py"
@@ -45,6 +46,9 @@ dependencies = [
 "bug tracker" = "https://github.com/GernotMaier/V2DL5/issues"
 
 [project.scripts]
+v2dl5-generate-runlist = "v2dl5.scripts.generate_runlist:main"
+v2dl5-print-simplified-light-curve = "v2dl5.scripts.print_simplified_light_curve:main"
+v2dl5-plot-binary-light-curves = "v2dl5.scripts.plot_binary_light_curves:main"
 v2dl5-reflected-region-analysis = "v2dl5.scripts.reflected_region_analysis:main"
 
 [tool.pytest.ini_options]

data/.gitignore renamed to v2dl5/data/.gitignore

@@ -62,14 +62,11 @@ def _read_fluxes_from_file(self, data_config):
         Read flux from file.
 
         """
-
         try:
-            if data_config["file_name"].endswith(".csv") or data_config["file_name"].endswith(
-                ".ecsv"
-            ):
+            if data_config["file_name"].endswith((".csv", ".ecsv")):
                 return self._read_fluxes_from_ecsv_file(data_config)
         except KeyError:
-            self._logger.error("File name not found in configuration file")
+            self._logger.error(f"File name not found in configuration {data_config}")
             raise KeyError
 
     def _add_orbital_parameters(self, data):
@@ -131,34 +128,38 @@ def _read_fluxes_from_ecsv_file(self, data_config, TimeMinMax=True, MJD_min=-1.0
         """
         table = Table.read(data_config["file_name"])
         f = {}
-        if not TimeMinMax:
-            table["time_min"] = table["time"].data
-            table["time_max"] = table["time"].data + 0.1
-
-        # MJD filter
-        condition = np.ones(len(table), dtype=bool)
-        if MJD_min > -1:
-            condition &= table["time_min"] > MJD_min
-        if MJD_max > -1:
-            condition &= table["time_max"] < MJD_max
-        table = table[condition]
-
-        f = {}
-        f["time_min"] = table["time_min"].data.tolist()
-        f["time_max"] = table["time_max"].data.tolist()
-        f["flux"] = table["flux"].data.flatten().tolist()
-        if "flux_err" in table.colnames:
-            f["flux_err"] = table["flux_err"].data.flatten().tolist()
-        else:
-            up = table["flux_up"].data.flatten().tolist()
-            down = table["flux_down"].data.flatten().tolist()
-            f["flux_err"] = [0.5 * abs(u - d) for u, d in zip(up, down)]
-        f["MJD"] = [0.5 * (a + b) for a, b in zip(f["time_min"], f["time_max"])]
-        f["MJD_err"] = [0.5 * (b - a) for a, b in zip(f["time_min"], f["time_max"])]
-        if "flux_ul" in table.colnames:
-            flux_ul = table["flux_ul"].data.flatten().tolist()
-            is_ul = table["is_ul"].data.flatten().tolist()
-            f["flux_ul"] = [flux if is_ul else -1.0 for flux, is_ul in zip(flux_ul, is_ul)]
-        else:
-            f["flux_ul"] = [-1.0 for _ in f["flux"]]
+        try:
+            if not TimeMinMax:
+                table["time_min"] = table["time"].data
+                table["time_max"] = table["time"].data + 0.1
+
+            # MJD filter
+            condition = np.ones(len(table), dtype=bool)
+            if MJD_min > -1:
+                condition &= table["time_min"] > MJD_min
+            if MJD_max > -1:
+                condition &= table["time_max"] < MJD_max
+            table = table[condition]
+
+            f = {}
+            f["time_min"] = table["time_min"].data.tolist()
+            f["time_max"] = table["time_max"].data.tolist()
+            f["flux"] = table["flux"].data.flatten().tolist()
+            if "flux_err" in table.colnames:
+                f["flux_err"] = table["flux_err"].data.flatten().tolist()
+            else:
+                up = table["flux_up"].data.flatten().tolist()
+                down = table["flux_down"].data.flatten().tolist()
+                f["flux_err"] = [0.5 * abs(u - d) for u, d in zip(up, down)]
+            f["MJD"] = [0.5 * (a + b) for a, b in zip(f["time_min"], f["time_max"])]
+            f["MJD_err"] = [0.5 * (b - a) for a, b in zip(f["time_min"], f["time_max"])]
+            if "flux_ul" in table.colnames:
+                flux_ul = table["flux_ul"].data.flatten().tolist()
+                is_ul = table["is_ul"].data.flatten().tolist()
+                f["flux_ul"] = [flux if is_ul else -1.0 for flux, is_ul in zip(flux_ul, is_ul)]
+            else:
+                f["flux_ul"] = [-1.0 for _ in f["flux"]]
+        except KeyError:
+            self._logger.error(f"Incomplete data file; key not found in {data_config['file_name']}")
+            raise KeyError
         return f

data/hip_mag9.fits.gz renamed to v2dl5/data/hip_mag9.fits.gz

@@ -66,10 +66,7 @@ def _parse():
 
 
 def main():
-    """
-    Binary light-curve plotting
-
-    """
+    """Binary light-curve plotting."""
     logging.root.setLevel(logging.INFO)
 
     args = _parse()

v2dl5/light_curves/data_reader.py

@@ -0,0 +1,52 @@
+#!/usr/bin/python
+"""
+Read gammapy-generated light curve files and print a simplified and readable version of the data.
+
+"""
+
+import argparse
+import logging
+
+import numpy as np
+from astropy.table import Table
+
+
+def extract_single_value(column):
+    return [
+        float(value[0]) if isinstance(value, (list, np.ndarray)) else float(value)
+        for value in column
+    ]
+
+
+def main():
+    """Main function to read and print light curve data."""
+    parser = argparse.ArgumentParser(description="Read and print light curve data.")
+
+    parser.add_argument(
+        "--file",
+        type=str,
+        required=True,
+        help="File (ecsv format) containing light curve data.",
+    )
+
+    args = parser.parse_args()
+
+    logging.basicConfig(level=logging.INFO)
+
+    complicated_table = Table.read(args.file, format="ascii.ecsv")
+    data_short = {
+        "e_min": extract_single_value(complicated_table["e_min"]),
+        "time_min": complicated_table["time_min"],
+        "time_max": complicated_table["time_max"],
+        "flux": extract_single_value(complicated_table["flux"]),
+        "flux_err": extract_single_value(complicated_table["flux_err"]),
+        "flux_ul": extract_single_value(complicated_table["flux_ul"]),
+        "significance": extract_single_value(complicated_table["sqrt_ts"]),
+    }
+
+    short_table = Table(data_short)
+    print(short_table)
+
+
+if __name__ == "__main__":
+    main()

v2dl5/scripts/plot_binary_light_curves.py

@@ -3,13 +3,13 @@
 """
 
 import logging
-import os
 
 from astropy import units as u
 from astropy.coordinates import Angle, SkyCoord, name_resolve
 from astropy.io import fits
 from astropy.table import Table
 from gammapy.maps import WcsGeom
+from importlib_resources import files
 from regions import CircleSkyRegion
 
 
@@ -155,6 +155,9 @@ def _read_bright_star_catalogue(self, exclusion_region_dict=None, max_wobble_dis
         """
         Read bright star catalogue from file.
 
+        Catalogue files are expected to be in the v2dl5/data directory and are distributed as part
+        of the v2dl5 package.
+
         Parameters
         ----------
         exclusion_region_dict: dict
@@ -176,7 +179,9 @@ def _read_bright_star_catalogue(self, exclusion_region_dict=None, max_wobble_dis
         self._logger.info(
             "Reading bright star catalogue from %s", exclusion_region_dict["star_file"]
         )
-        hip = fits.open(os.path.expandvars(exclusion_region_dict["star_file"]))
+        star_file = files("v2dl5.data").joinpath("data/" + exclusion_region_dict["star_file"])
+        hip = fits.open(star_file)
+
         catalogue = Table(hip[1].data)
         catalogue = catalogue[
             catalogue["Vmag"] + catalogue["B-V"] < exclusion_region_dict["magnitude_B"]