Merge pull request #32 from GernotMaier/light-curve-plotting

Improved light curve plotting

Author: GernotMaier

README.md

@@ -36,7 +36,7 @@ Activate the environment to start the analysis:
 
 ```bash
 conda activate v2dl5
-export PYTHONPATH=$PYTHONPATH:"${PWD}"
+pip install -e .
 ```
 
 ## Run list generator
@@ -46,10 +46,27 @@ It generates a lot of printout which should be used to fine tune the run selecti
 
 Example:
 
-```bash
-
+```console
 python v2dl5/scripts/generate_runlist.py \
     --obs_table ../../../VTS/DL3/v490/dl3_pointlike_moderate2tel/obs-index.fits.gz \
     --config examples/run_selection.yml \
    --output_dir my_output_dir
 ```
+
+## Reflected region analysis
+
+```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 \
+    --config examples/reflected_region_analysis.yml \
+    --output_dir my_output_dir
+```
+
+## Binary light curve plotting
+
+```console
+python v2dl5/scripts/plot_binary_light_curves.py \
+    --instrument VERITAS \
+    --configuration examples/binary_lightcurve_plotting.yml
+```

examples/binary_lightcurve_plotting.yml

@@ -0,0 +1,8 @@
+data:
+    - instrument: VERITAS
+      file_name: results/light_curve_nightly.ecsv
+      plot_instrument: true
+      plot_variable: flux
+      flux_axis_label: "Flux E>350 GeV (1/(cm$^2$s))"
+      marker_type: 'o'
+      marker_color: 'r'

pyproject.toml

@@ -30,7 +30,7 @@ dependencies = [
     "astroquery",
     "gammapy",
     "Ipython",
-    "pydantic<2",
+    "pydantic",
 ]
 
 [project.optional-dependencies]

v2dl5/light_curves/binary_plotting.py

@@ -50,11 +50,22 @@ def plot_flux_vs_time(self, time_axis="MJD", mjd_min=None, mjd_max=None, file_ty
         colors = plotting_utilities.get_color_list(len(self.data))
 
         for idx, (instrument, data) in enumerate(self.data.items()):
-            x, y, e = self._get_light_curve_in_MJD_limits(data, time_axis, mjd_min, mjd_max)
+            x, y, e, x_ul, y_ul = self._get_light_curve_in_MJD_limits(
+                data, time_axis, mjd_min, mjd_max
+            )
 
             if instrument in self.config and self.config["plot_instrument"] is False:
                 continue
-
+            color = (
+                colors[idx]
+                if self.config[idx].get("marker_color") is None
+                else self.config[idx]["marker_color"]
+            )
+            marker = (
+                plotting_utilities.get_marker_list()[idx]
+                if self.config[idx].get("marker_type") is None
+                else self.config[idx]["marker_type"]
+            )
             plt.errorbar(
                 x,
                 y,
@@ -65,21 +76,26 @@ def plot_flux_vs_time(self, time_axis="MJD", mjd_min=None, mjd_max=None, file_ty
                     if self.config[idx].get("plot_label") is None
                     else self.config[idx]["plot_label"]
                 ),
-                color=(
-                    colors[idx]
-                    if self.config[idx].get("marker_color") is None
-                    else self.config[idx]["marker_color"]
-                ),
-                marker=(
-                    plotting_utilities.get_marker_list()[idx]
-                    if self.config[idx].get("marker_type") is None
-                    else self.config[idx]["marker_type"]
-                ),
+                color=color,
+                marker=marker,
                 linestyle="none",
                 fillstyle="none",
                 linewidth=plotting_utilities.get_line_width(),
                 markersize=plotting_utilities.get_marker_size(),
             )
+            if len(y_ul) > 0:
+                plt.errorbar(
+                    x_ul,
+                    y_ul,
+                    yerr=0.1 * max(y_ul),
+                    color=color,
+                    fmt="_",
+                    linestyle="none",
+                    fillstyle="none",
+                    uplims=True,
+                    linewidth=plotting_utilities.get_line_width(),
+                    markersize=plotting_utilities.get_marker_size(),
+                )
 
         plt.xlabel(self._get_time_axis_label(time_axis))
         if mjd_min is not None and mjd_max is not None:
@@ -127,21 +143,28 @@ def _get_light_curve_in_MJD_limits(self, data, time_axis, mjd_min, mjd_max):
         x = []
         y = []
         e = []
+        x_ul = []
+        y_ul = []
         for i in range(len(data["MJD"])):
             if mjd_min is not None and data["MJD"][i] < mjd_min:
                 continue
             if mjd_max is not None and data["MJD"][i] > mjd_max:
                 continue
             if time_axis == "MJD":
-                x.append(data["MJD"][i])
+                w_x = data["MJD"][i]
             elif time_axis == "orbital phase":
-                x.append(data["phase"][i])
+                w_x = data["phase"][i]
             else:
                 self._logger.error(f"Unknown time axis: {time_axis}")
                 raise ValueError
-            y.append(data["flux"][i])
-            e.append(data["flux_err"][i])
-        return x, y, e
+            if "flux_ul" in data and data["flux_ul"][i] > 0:
+                x_ul.append(w_x)
+                y_ul.append(data["flux_ul"][i])
+            elif "flux_ul" not in data or data["flux_ul"][i] < 0:
+                x.append(w_x)
+                y.append(data["flux"][i])
+                e.append(data["flux_err"][i])
+        return x, y, e, x_ul, y_ul
 
 
 #

v2dl5/light_curves/data_reader.py

@@ -1,8 +1,9 @@
-import csv
 import logging
 
 import astropy.units as u
+import numpy as np
 import yaml
+from astropy.table import Table
 
 import v2dl5.light_curves.orbital_period as orbit
 
@@ -66,65 +67,11 @@ def _read_fluxes_from_file(self, data_config):
             if data_config["file_name"].endswith(".csv") or data_config["file_name"].endswith(
                 ".ecsv"
             ):
-                return self._read_fluxes_from_csv_file(data_config)
+                return self._read_fluxes_from_ecsv_file(data_config)
         except KeyError:
             self._logger.error("File name not found in configuration file")
             raise KeyError
 
-    def _read_fluxes_from_csv_file(self, data_config, TimeMinMax=True, MJD_min=-1.0, MJD_max=-1.0):
-        """
-        Read gamma-ray fluxes from csv file (open gamma-ray format)
-        - ignore all lines with '#'
-        - accept a random number of spaces as delimiter
-
-        Parameters:
-        -----------
-        data_config: dict
-            configuration dictionary
-        TimeMinMax: bool
-            flag to read time_min and time_max
-        MJD_min: float
-            MJD min value for MJD cut
-        MJD_max: float
-            MJD max value for MJD cut
-
-        """
-        fp = open(data_config["file_name"])
-        rdr = csv.DictReader(
-            filter(lambda row: row[0] != "#", fp), delimiter=" ", skipinitialspace=True
-        )
-        f = {}
-        f["time_min"] = []
-        f["time_max"] = []
-        f["flux"] = []
-        f["flux_err"] = []
-        for row in rdr:
-            if TimeMinMax:
-                t_min = float(row["time_min"])
-                t_max = float(row["time_max"])
-            else:
-                t_min = float(row["time"])
-                t_max = float(row["time"]) + 0.1
-            if MJD_min > 0 and t_min < MJD_min:
-                continue
-            if MJD_max > 0 and t_max > MJD_max:
-                continue
-
-            f["time_min"].append(t_min)
-            f["time_max"].append(t_max)
-
-            f["flux"].append(float(row["flux"]))
-            if "flux_err" in row:
-                f["flux_err"].append(float(row["flux_err"]))
-            elif "flux_up" in row:
-                f["flux_err"].append(0.5 * abs(float(row["flux_up"]) - float(row["flux_down"])))
-        fp.close()
-
-        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"])]
-
-        return f
-
     def _add_orbital_parameters(self, data):
         """
         Add orbital phase to light-curve data data.
@@ -165,3 +112,53 @@ def convert_photon_to_energy_flux(C_v, C_e, E_0, gamma):
         # conversion to erg
         f = f * (E_0.to(u.erg)).value
         return [v * f for v in C_v], [e * f for e in C_e]
+
+    def _read_fluxes_from_ecsv_file(self, data_config, TimeMinMax=True, MJD_min=-1.0, MJD_max=-1.0):
+        """
+        Read gamma-ray fluxes from ecsv file (open gamma-ray format)
+
+        Parameters:
+        -----------
+        data_config: dict
+            configuration dictionary
+        TimeMinMax: bool
+            flag to read time_min and time_max
+        MJD_min: float
+            MJD min value for MJD cut
+        MJD_max: float
+            MJD max value for MJD cut
+
+        """
+        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"]]
+        return f

v2dl5/plot.py

@@ -207,7 +207,7 @@ def plot_light_curve(self, light_curve, plot_name):
         )
 
         try:
-            light_curve.plot(ax=ax, marker="o", label=plot_name, sed_type="flux")
+            light_curve.plot(ax=ax, marker="o", label=plot_name, sed_type="flux", time_format="mjd")
             ax.set_yscale("linear")
             self._plot(
                 plot_name="light_curve_" + plot_name.replace(" ", "_"), output_dir=self.output_dir