From 4577655eaf3368f56183c452004a25552fea2ad4 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Tue, 4 Feb 2025 15:59:23 +0100 Subject: [PATCH 1/8] config.yaml: cloud_correction block was modified lst_m1_m2_cloud_correction.py: interpolation function and additional cleaning procedure were added --- magicctapipe/resources/config.yaml | 9 +- .../lst1_magic/lst_m1_m2_cloud_correction.py | 453 ++++++++++++++++-- 2 files changed, 414 insertions(+), 48 deletions(-) diff --git a/magicctapipe/resources/config.yaml b/magicctapipe/resources/config.yaml index 8b4e14c6..028fba95 100644 --- a/magicctapipe/resources/config.yaml +++ b/magicctapipe/resources/config.yaml @@ -279,9 +279,10 @@ dl2_to_dl3: source_dec: null # used when the source name cannot be resolved cloud_correction: - specify_cloud: yes - base_height: "5900. m" - thickness: "1320. m" - vertical_transmission: 1.0 + max_gap_lidar_shots: "900 s" + lidar_report_file: "/home/natti/CTA/lodz_cluster/lstchain_10_11_test/scripts/pic_server/lidar_reports_clouds_ST.03.16.yaml" number_of_layers: 10 + cleaning_multiplication_factor: 1.25 + use_additional_cleaning: True + diff --git a/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py b/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py index f7655a98..35de2e7a 100644 --- a/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py +++ b/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py @@ -12,22 +12,32 @@ """ import argparse import logging +import os +import time +from pathlib import Path import astropy.units as u +import ctapipe import numpy as np import pandas as pd +import tables import yaml from astropy.coordinates import AltAz, SkyCoord +from astropy.time import Time from ctapipe.coordinates import TelescopeFrame from ctapipe.image import ( concentration_parameters, hillas_parameters, leakage_parameters, + tailcuts_clean, timing_parameters, ) from ctapipe.instrument import SubarrayDescription from ctapipe.io import read_table +from scipy.interpolate import interp1d +import magicctapipe +from magicctapipe.image import MAGICClean from magicctapipe.io import save_pandas_data_in_table logger = logging.getLogger(__name__) @@ -37,7 +47,7 @@ def model0(imp, h, zd): """ - Calculated the geometrical part of the model relating the emission height with the angular distance from the arrival direction + Calculates the geometrical part of the model relating the emission height with the angular distance from the arrival direction Parameters ---------- @@ -53,7 +63,7 @@ def model0(imp, h, zd): numpy ndarray Angular distance in units of degree """ - d = h / np.cos(zd) + d = h / np.cos(np.deg2rad(zd)) return np.arctan((imp / d).to("")).to_value("deg") @@ -107,24 +117,237 @@ def trans_height(x, Hc, dHc, trans): return t -def process_telescope_data(input_file, config, tel_id, camgeom, focal_eff): +def lidar_cloud_interpolation( + mean_subrun_timestamp, max_gap_lidar_shots, lidar_report_file +): + """ + Retrieves or interpolates LIDAR cloud parameters based on the closest timestamps to an input mean timestamp of the processed subrun. + + Parameters + ----------- + mean_subrun_timestamp : int or float + The mean timestamp of the processed subrun (format: unix). + + max_gap_lidar_shots : int or float + Maximum allowed time gap for interpolation (in seconds). + + lidar_report_file : str + Path to the yaml file containing LIDAR laser reports with columns: + - "timestamp" (format: ISO 8601) + - "base_height", "top_height", "transmission", "lidar_zenith" + + Returns + -------- + tuple or None + + A tuple containing interpolated or nearest values for: + - base_height (float): The base height of the cloud layer in meters. + - top_height (float): The top height of the cloud layer in meters. + - vertical_transmission (float): Transmission factor of the cloud layer adjusted for zenith angle. + + If no nodes are found within the maximum allowed time gap, returns None. + """ + + if not os.path.isfile(lidar_report_file): + raise FileNotFoundError(f"LIDAR report file not found: {lidar_report_file}") + + with open(lidar_report_file, "r") as f: + data = yaml.safe_load(f) + + records = [] + for entry in data["data"]: + timestamp = pd.to_datetime(entry["timestamp"], errors="coerce") + lidar_zenith = entry["lidar_zenith"] + + lowest_transmission_layer = min( + entry["layers"], key=lambda layer: layer["transmission"] + ) + + records.append( + { + "timestamp": timestamp, + "lidar_zenith": lidar_zenith, + "base_height": lowest_transmission_layer["base_height"], + "top_height": lowest_transmission_layer["top_height"], + "transmission": lowest_transmission_layer["transmission"], + } + ) + + df = pd.DataFrame(records) + + df["timestamp"] = pd.to_datetime(df["timestamp"], format="ISO8601") + df["unix_timestamp"] = df["timestamp"].astype(np.int64) / 10**9 + df["time_diff"] = df["unix_timestamp"] - mean_subrun_timestamp + + df["lidar_zenith"] = (pd.to_numeric(df["lidar_zenith"], errors="coerce")).to_numpy() + vertical_transmission = df["transmission"] ** np.cos(np.deg2rad(df["lidar_zenith"])) + df["vertical_transmission"] = vertical_transmission + + closest_node_before = df[ + (df["time_diff"] < 0) & (np.abs(df["time_diff"]) <= max_gap_lidar_shots) + ].nlargest(1, "time_diff") + closest_node_after = df[ + (df["time_diff"] > 0) & (np.abs(df["time_diff"]) <= max_gap_lidar_shots) + ].nsmallest(1, "time_diff") + + # Check whether the conditions for interpolation are met or not + if not closest_node_before.empty and not closest_node_after.empty: + node_before = closest_node_before.iloc[0] + node_after = closest_node_after.iloc[0] + logger.info( + f"\nFound suitable interpolation nodes within the allowed temporal gap for timestamp {Time(mean_subrun_timestamp, format='unix').iso}" + f"\nUsing following interpolation nodes:\n" + f"\n******************** Node before ******************* \n{node_before}" + f"\n******************** Node after ******************** \n{node_after}" + f"\n\nInterpolation results:" + ) + + interp_values = {} + for param in ["base_height", "top_height", "vertical_transmission"]: + interp_func = interp1d( + [node_before["unix_timestamp"], node_after["unix_timestamp"]], + [node_before[param], node_after[param]], + kind="linear", + bounds_error=False, + ) + interp_values[param] = interp_func(mean_subrun_timestamp) + logger.info(f"\t {param}: {interp_values[param]:.4f}") + + return ( + interp_values["base_height"], + interp_values["top_height"], + interp_values["vertical_transmission"], + ) + + # Handle cases where only one node is available + closest_node = ( + closest_node_before if not closest_node_before.empty else closest_node_after + ) + + if closest_node is not None and not closest_node.empty: + closest = closest_node.iloc[0] + logger.info( + f"\nOnly one suitable LIDAR report found for timestamp {Time(mean_subrun_timestamp, format='unix').iso} " + f"within the maximum allowed temporal gap. \nSkipping interpolation. Using nearest node values instead." + f"\n\n{closest}" + ) + return ( + closest["base_height"], + closest["top_height"], + closest["vertical_transmission"], + ) + + logger.info( + f"\nNo node is within the maximum allowed temporal gap for timestamp {Time(mean_subrun_timestamp, format='unix').iso}. Exiting ..." + ) + return None + + +def clean_image_with_modified_thresholds( + event_image, event_pulse_time, unsuitable_mask, magic_clean, config_clean_magic, cmf +): + + """ + Creates a wrapper function to modify the thresholds + + Parameters + ---------- + event_image : numpy.ndarray + Input array with event image + + event_pulse_time : numpy.ndarray + Input array with event times + + unsuitable_mask : numpy.ndarray + Array of unsuitable pixels + + magic_clean : magicctapipe.image.cleaning.MAGICClean + Cleaning implementation used by MAGIC + + config_clean_magic : dict + Cleaning parameters read from the config file + + cmf : float + Multiplication factor for additional cleaning + + Returns + ------- + clean_mask : numpy.ndarray + Mask with pixels surviving the cleaning + + image : numpy.ndarray + Image with surviving pixels + + peak_time : numpy.ndarray + Times with only surviving pixels + """ + + # Multiply the thresholds by cmf + modified_picture_thresh = config_clean_magic["picture_thresh"] * cmf + modified_boundary_thresh = config_clean_magic["boundary_thresh"] * cmf + + # Directly modify magic_clean instance variables with the modified thresholds + magic_clean.picture_thresh = modified_picture_thresh + magic_clean.boundary_thresh = modified_boundary_thresh + + # Now call the clean_image method on the magic_clean instance with the modified thresholds + clean_mask, image, peak_time = magic_clean.clean_image( + event_image=event_image, + event_pulse_time=event_pulse_time, + unsuitable_mask=unsuitable_mask, + ) + + return clean_mask, image, peak_time + + +def process_telescope_data( + dl1_params, + dl1_images, + config, + tel_id, + tel_ids, + camgeom, + focal_eff, + nlayers, + mean_subrun_zenith, + Hc, + dHc, + trans, + cmf, +): """ Corrects LST-1 and MAGIC data affected by a cloud presence Parameters ---------- - input_file : str - Path to an input .h5 DL1 file + dl1_params : str + Path to an input .h5 DL1 table with parameters + dl1_images : str + Path to an input .h5 DL1 table with images config : dict Configuration for the LST-1 + MAGIC analysis tel_id : numpy.int16 LST-1 and MAGIC telescope ids + tel_ids : dict + List of LST-1 and MAGIC telescope names and ids from config file camgeom : ctapipe.instrument.camera.geometry.CameraGeometry An instance of the CameraGeometry class containing information about the camera's configuration, including pixel type, number of pixels, rotation angles, and the reference frame. focal_eff : astropy.units.quantity.Quantity Effective focal length + nlayers : astropy.units.quantity.Quantity + Array with heights of each cloud layer a.g.l. + mean_subrun_zenith : numpy.float64 + Mean value of zenith per subran + Hc : astropy.units.quantity.Quantity + Height of the base of the cloud a.g.l. + dHc : astropy.units.quantity.Quantity + Cloud thickness + trans : numpy.float64 + Transmission of the cloud + cmf : float + Multiplication factor for additional cleaning Returns ------- @@ -132,21 +355,28 @@ def process_telescope_data(input_file, config, tel_id, camgeom, focal_eff): Data frame of corrected DL1 parameters """ - assigned_tel_ids = config["mc_tel_ids"] + # AC LST + cleaning_level_lst = config["LST"]["tailcuts_clean"] + modified_boundary_thresh_lst = cleaning_level_lst["boundary_thresh"] * cmf + modified_picture_thresh_lst = cleaning_level_lst["picture_thresh"] * cmf + # AC MAGIC + # Configure the MAGIC image cleaning + config_clean_magic = config["MAGIC"]["magic_clean"] + magic_clean = MAGICClean(camgeom, config_clean_magic) - correction_params = config.get("cloud_correction", {}) - Hc = u.Quantity(correction_params.get("base_height")) - dHc = u.Quantity(correction_params.get("thickness")) - trans0 = correction_params.get("vertical_transmission") - nlayers = correction_params.get("number_of_layers") + unsuitable_mask = None all_params_list = [] - dl1_params = read_table(input_file, "/events/parameters") - dl1_images = read_table(input_file, "/events/dl1/image_" + str(tel_id)) + if dl1_images is None: + return None m2deg = np.rad2deg(1) / focal_eff * u.degree + Hcl = np.linspace(Hc, Hc + dHc, nlayers) # position of each layer + transl = trans_height(Hcl, Hc, dHc, trans) # transmissions of each layer + transl = np.append(transl, transl[-1]) + inds = np.where( np.logical_and(dl1_params["intensity"] > 0.0, dl1_params["tel_id"] == tel_id) )[0] @@ -161,29 +391,21 @@ def process_telescope_data(input_file, config, tel_id, camgeom, focal_eff): multiplicity = dl1_params["multiplicity"][index] combo_type = dl1_params["combo_type"][index] - if assigned_tel_ids["LST-1"] == tel_id: + if tel_ids["LST-1"] == tel_id: event_id_image, obs_id_image = event_id_lst, obs_id_lst else: event_id_image, obs_id_image = event_id_magic, obs_id_magic pointing_az = dl1_params["pointing_az"][index] pointing_alt = dl1_params["pointing_alt"][index] - zenith = 90.0 - np.rad2deg(pointing_alt) - psi = dl1_params["psi"][index] * u.deg time_diff = dl1_params["time_diff"][index] n_islands = dl1_params["n_islands"][index] signal_pixels = dl1_params["n_pixels"][index] - trans = trans0 ** (1 / np.cos(zenith)) - Hcl = np.linspace(Hc, Hc + dHc, nlayers) # position of each layer - transl = trans_height(Hcl, Hc, dHc, trans) # transmissions of each layer - transl = np.append(transl, transl[-1]) - alt_rad = np.deg2rad(dl1_params["alt"][index]) az_rad = np.deg2rad(dl1_params["az"][index]) impact = dl1_params["impact"][index] * u.m - psi = dl1_params["psi"][index] * u.deg cog_x = (dl1_params["x"][index] * m2deg).value * u.deg cog_y = (dl1_params["y"][index] * m2deg).value * u.deg @@ -205,6 +427,8 @@ def process_telescope_data(input_file, config, tel_id, camgeom, focal_eff): src_x, src_y = -src_y, -src_x cog_x, cog_y = -cog_y, -cog_x + psi = np.arctan2(src_x - cog_x, src_y - cog_y) + pix_x_tel = (camgeom.pix_x * m2deg).to(u.deg) pix_y_tel = (camgeom.pix_y * m2deg).to(u.deg) @@ -217,7 +441,7 @@ def process_telescope_data(input_file, config, tel_id, camgeom, focal_eff): d2_src_pix = (src_x - pix_x_tel) ** 2 + (src_y - pix_y_tel) ** 2 distance[d2_cog_pix > d2_cog_src + d2_src_pix] = 0 - dist_corr_layer = model2(impact, Hcl, zenith) * u.deg + dist_corr_layer = model2(impact, Hcl, mean_subrun_zenith) * u.deg ilayer = np.digitize(distance, dist_corr_layer) trans_pixels = transl[ilayer] @@ -232,37 +456,85 @@ def process_telescope_data(input_file, config, tel_id, camgeom, focal_eff): raise ValueError("Error: 'inds_img' list is empty!") index_img = inds_img[0] image = dl1_images["image"][index_img] - cleanmask = dl1_images["image_mask"][index_img] + clean_mask = dl1_images["image_mask"][index_img] peak_time = dl1_images["peak_time"][index_img] image /= trans_pixels - corr_image = image.copy() - clean_image = corr_image[cleanmask] - clean_camgeom = camgeom[cleanmask] + # additional cleaning + if config["cloud_correction"]["use_additional_cleaning"]: + clean_mask = np.ones_like( + image, dtype=bool + ) # Assuming full mask if not defined + if tel_ids["LST-1"] == tel_id: + clean = tailcuts_clean( + camgeom, + image, + boundary_thresh=modified_boundary_thresh_lst, + picture_thresh=modified_picture_thresh_lst, + min_number_picture_neighbors=cleaning_level_lst[ + "min_number_picture_neighbors" + ], + ) + # Create a mask indicating which pixels survived the cleaning + clean_mask = np.zeros_like(image, dtype=bool) + clean_mask[clean] = True + + if np.sum(clean_mask) == 0: + continue + + # Apply the mask to relevant data arrays + clean_peak_time = peak_time[clean_mask] + clean_image = image[clean_mask] + clean_camgeom = camgeom[clean_mask] + + elif tel_ids["MAGIC-I"] == tel_id or tel_ids["MAGIC-II"] == tel_id: + # Use the wrapper function to clean the image + clean_mask, image, peak_time = clean_image_with_modified_thresholds( + event_image=image, + event_pulse_time=peak_time, + unsuitable_mask=unsuitable_mask, + magic_clean=magic_clean, # Pass the magic_clean instance here + config_clean_magic=config_clean_magic, + cmf=cmf, + ) + + clean_camgeom = camgeom[clean_mask] + clean_image = image[clean_mask] + clean_peak_time = peak_time[clean_mask] + + # Check if clean_image is empty or all zeros + if clean_image.size == 0 or not np.any(clean_image): + continue # Skip this iteration if clean_image is empty or all zeros + + else: + clean_image = image[clean_mask] + clean_camgeom = camgeom[clean_mask] + clean_peak_time = peak_time[clean_mask] + # re-calculation of dl1 parameters hillas_params = hillas_parameters(clean_camgeom, clean_image) timing_params = timing_parameters( - clean_camgeom, clean_image, peak_time[cleanmask], hillas_params + clean_camgeom, clean_image, clean_peak_time, hillas_params ) - leakage_params = leakage_parameters(camgeom, corr_image, cleanmask) - if assigned_tel_ids["LST-1"] == tel_id: + leakage_params = leakage_parameters(camgeom, image, clean_mask) + if tel_ids["LST-1"] == tel_id: conc_params = concentration_parameters( clean_camgeom, clean_image, hillas_params - ) # "For LST-1 we compute concentration from the cleaned image and for MAGIC from the full image to reproduce the current behaviour in the standard code + ) # For LST-1 we compute concentration from the cleaned image and for MAGIC from the full image to reproduce the current behaviour in the standard code else: - conc_params = concentration_parameters(camgeom, corr_image, hillas_params) + conc_params = concentration_parameters(camgeom, image, hillas_params) + + prefixed_conc_params = { + f"concentration_{key}": value for key, value in conc_params.items() + } event_params = { **hillas_params, **timing_params, **leakage_params, + **prefixed_conc_params, } - prefixed_conc_params = { - f"concentration_{key}": value for key, value in conc_params.items() - } - event_params.update(prefixed_conc_params) - event_info_dict = { "obs_id": obs_id, "event_id": event_id, @@ -290,6 +562,7 @@ def process_telescope_data(input_file, config, tel_id, camgeom, focal_eff): def main(): """Main function.""" + start_time = time.time() parser = argparse.ArgumentParser() parser.add_argument( @@ -302,9 +575,9 @@ def main(): ) parser.add_argument( - "--output_file", + "--output_dir", "-o", - dest="output_file", + dest="output_dir", type=str, default="./data", help="Path to a directory where to save an output corrected DL1 file", @@ -345,14 +618,58 @@ def main(): tel_ids = config["mc_tel_ids"] - dfs = [] # Initialize an empty list to store DataFrames + correction_params = config.get("cloud_correction", {}) + max_gap_lidar_shots = u.Quantity(correction_params.get("max_gap_lidar_shots")) + lidar_report_file = correction_params.get( + "lidar_report_file" + ) # path to the lidar report + nlayers = correction_params.get("number_of_layers") + cmf = correction_params.get("cleaning_multiplication_factor") + + dl1_params = read_table(args.input_file, "/events/parameters") + + mean_subrun_timestamp = np.mean(dl1_params["timestamp"]) + mean_subrun_zenith = np.mean(90.0 - np.rad2deg(dl1_params["pointing_alt"])) + + cloud_params = lidar_cloud_interpolation( + mean_subrun_timestamp, max_gap_lidar_shots, lidar_report_file + ) + + Hc = u.Quantity(cloud_params[0], u.m) + dHc = u.Quantity(cloud_params[1] - cloud_params[0], u.m) + vertical_trans = u.Quantity(cloud_params[2]) + trans = vertical_trans ** (1 / np.cos(np.deg2rad(mean_subrun_zenith))) + + dfs = [] for tel_name, tel_id in tel_ids.items(): if tel_id != 0: # Only process telescopes that have a non-zero ID + # Read images for each telescope + image_node_path = "/events/dl1/image_" + str(tel_id) + try: + dl1_images = read_table(args.input_file, image_node_path) + except tables.NoSuchNodeError: + raise RuntimeError( + f"Fatal error: No image found for telescope with ID {tel_id}." + ) + df = process_telescope_data( - args.input_file, config, tel_id, camgeom[tel_id], focal_eff[tel_id] + dl1_params, + dl1_images, + config, + tel_id, + tel_ids, + camgeom[tel_id], + focal_eff[tel_id], + nlayers, + mean_subrun_zenith, + Hc, + dHc, + trans, + cmf, ) - dfs.append(df) + if df is not None: + dfs.append(df) df_all = pd.concat(dfs, ignore_index=True) @@ -374,16 +691,64 @@ def main(): lambda x: x.value if isinstance(x, u.Quantity) else x ) + df_all["psi"] = np.degrees(df_all["psi"]) + df_all["phi"] = np.degrees(df_all["phi"]) + for col in columns_to_convert: df_all[col] = pd.to_numeric(df_all[col], errors="coerce") df_all = df_all.drop(columns=["deviation"], errors="ignore") + Path(args.output_dir).mkdir(exist_ok=True, parents=True) + input_file_name = Path(args.input_file).name + output_file_name = input_file_name.replace("dl1_stereo", "dl1_corr") + output_file = f"{args.output_dir}/{output_file_name}" + save_pandas_data_in_table( - df_all, args.output_file, group_name="/events", table_name="parameters" + df_all, output_file, group_name="/events", table_name="parameters" ) - subarray_info.to_hdf(args.output_file) + with tables.open_file(output_file, mode="a") as f: + cloud_metadata_group = f.create_group("/", "weather", "Cloud parameters") + + cloud_base_height_data = Hc.to_value(u.m) + cloud_base_height_array = f.create_array( + cloud_metadata_group, + "cloud_base_height", + np.array([cloud_base_height_data]), + ) + cloud_base_height_array.attrs["unit"] = str(u.m) + + cloud_thickness_data = dHc.to_value(u.m) + cloud_thickness_array = f.create_array( + cloud_metadata_group, "cloud_thickness", np.array([cloud_thickness_data]) + ) + cloud_thickness_array.attrs["unit"] = str(u.m) + + cloud_vertical_trans_data = vertical_trans.to_value(u.dimensionless_unscaled) + cloud_vertical_trans_array = f.create_array( + cloud_metadata_group, + "cloud_vertical_transmission", + np.array([cloud_vertical_trans_data]), + ) + cloud_vertical_trans_array.attrs["unit"] = str(u.dimensionless_unscaled) + + subarray_info.to_hdf(output_file) + + logger.info(f"Correction parameters: {correction_params}") + logger.info(f"ctapipe version: {ctapipe.__version__}") + logger.info(f"magicctapipe version: {magicctapipe.__version__}") + + process_time = time.time() - start_time + logger.info(f"\nProcess time: {process_time:.0f} [sec]\n") + logger.info(f"\nOutput file: {output_file}") + + logger.info( + f"\n******************** Cloud parameters ******************* \n" + f"base_height: {cloud_params[0]:.2f} \n" + f"top_height: {cloud_params[1]:.2f}\n" + f"vertical_transmission: {cloud_params[2]:.6f}" + ) logger.info("\nDone.") From 17c019e93fb6ac4fdc12217a512d32efbdc6e46e Mon Sep 17 00:00:00 2001 From: Mario Pecimotika <45829235+mpecimotika@users.noreply.github.com> Date: Wed, 5 Feb 2025 09:56:39 +0100 Subject: [PATCH 2/8] Updated error_codes.py with new code for cloud correction error --- magicctapipe/utils/error_codes.py | 1 + 1 file changed, 1 insertion(+) diff --git a/magicctapipe/utils/error_codes.py b/magicctapipe/utils/error_codes.py index b8a15bf2..e8cceab8 100644 --- a/magicctapipe/utils/error_codes.py +++ b/magicctapipe/utils/error_codes.py @@ -9,3 +9,4 @@ GENERIC_ERROR_CODE = 1 NO_EVENTS_WITHIN_MAXIMUM_DISTANCE = 2 NO_COINCIDENT_EVENTS = 3 +NO_PROPER_LIDAR_REPORTS = 4 From 09e4589b2de6f13f9df6ce6e758edd26f339bc18 Mon Sep 17 00:00:00 2001 From: Mario Pecimotika <45829235+mpecimotika@users.noreply.github.com> Date: Wed, 5 Feb 2025 15:51:55 +0100 Subject: [PATCH 3/8] interpolation function update; Include NO_PROPER_LIDAR_REPORTS in __init__.py --- .gitignore | 1 + magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py | 4 +++- magicctapipe/utils/__init__.py | 2 ++ 3 files changed, 6 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 08120e67..f0ed6323 100644 --- a/.gitignore +++ b/.gitignore @@ -52,6 +52,7 @@ distribute-*.tar.gz # Other generated files htmlcov .coverage +.DS_Store # Sphinx documentation docs/api diff --git a/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py b/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py index 35de2e7a..6e645929 100644 --- a/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py +++ b/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py @@ -13,6 +13,7 @@ import argparse import logging import os +import sys import time from pathlib import Path @@ -39,6 +40,7 @@ import magicctapipe from magicctapipe.image import MAGICClean from magicctapipe.io import save_pandas_data_in_table +from magicctapipe.utils import NO_PROPER_LIDAR_REPORTS logger = logging.getLogger(__name__) logger.addHandler(logging.StreamHandler()) @@ -240,7 +242,7 @@ def lidar_cloud_interpolation( logger.info( f"\nNo node is within the maximum allowed temporal gap for timestamp {Time(mean_subrun_timestamp, format='unix').iso}. Exiting ..." ) - return None + sys.exit(NO_PROPER_LIDAR_REPORTS) def clean_image_with_modified_thresholds( diff --git a/magicctapipe/utils/__init__.py b/magicctapipe/utils/__init__.py index fc2e4a86..0cc66113 100644 --- a/magicctapipe/utils/__init__.py +++ b/magicctapipe/utils/__init__.py @@ -4,6 +4,7 @@ GENERIC_ERROR_CODE, NO_COINCIDENT_EVENTS, NO_EVENTS_WITHIN_MAXIMUM_DISTANCE, + NO_PROPER_LIDAR_REPORTS, ) from .functions import ( HEIGHT_ORM, @@ -28,6 +29,7 @@ "GENERIC_ERROR_CODE", "NO_COINCIDENT_EVENTS", "NO_EVENTS_WITHIN_MAXIMUM_DISTANCE", + "NO_PROPER_LIDAR_REPORTS", "identify_time_edges", "intersect_time_intervals", "GTIGenerator", From 465d6a93b24248281f89a3329f94e6f07dd2cadf Mon Sep 17 00:00:00 2001 From: Mario Pecimotika <45829235+mpecimotika@users.noreply.github.com> Date: Sun, 9 Feb 2025 15:09:12 +0100 Subject: [PATCH 4/8] script to transform MAGIC LIDAR reports in yaml format needed for cloud correction --- .../lst1_magic_lidar_reports_to_yaml.py | 235 ++++++++++++++++++ 1 file changed, 235 insertions(+) create mode 100644 magicctapipe/scripts/lst1_magic/lst1_magic_lidar_reports_to_yaml.py diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_lidar_reports_to_yaml.py b/magicctapipe/scripts/lst1_magic/lst1_magic_lidar_reports_to_yaml.py new file mode 100644 index 00000000..4a90add9 --- /dev/null +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_lidar_reports_to_yaml.py @@ -0,0 +1,235 @@ +#!/usr/bin/env python +# coding: utf-8 + +""" +This script processes LIDAR report files obtained with MARS-based software quate to extract cloud-reveant parameters present during the observations with MAGIC telescopes (if any), which can then be used to correct event images affected by the clouds using lst1_magic_cloud_correction.py script. + +Usage: +$ python lst1_magic_lidar_report_to_yaml.py +--input-file "magic_lidar.*.results" +--output-file "magic_lidar_report_clouds_CrabNebula.yaml" +""" + +import os +import re +import glob +import argparse +from datetime import datetime +from ruamel.yaml import YAML +from ruamel.yaml.scalarstring import SingleQuotedScalarString as SQS +import logging + +logger = logging.getLogger(__name__) +logger.addHandler(logging.StreamHandler()) +logger.setLevel(logging.INFO) + +def extract_data_from_file(input_filename): + """ + Extracts cloud-relevant data from a LIDAR report file (if any). + + Parameters + ---------- + input_filename : str + The name of an input LIDAR report file to process. + + Returns + ------- + dict or None + A dictionary with extracted data if valid cloud data is found, or None otherwise. + """ + + try: + with open(input_filename, "r") as file: + content = file.read() + logger.info(f"\n----------------------------------------------------" + f"\nProcessing file: {input_filename}") + + if not re.search(r"ERROR_CODE:\s*4", content): + logger.info("No valid LIDAR report found (ERROR_CODE != 4). \nSkipping file.") + return None + + zenith_match = re.search(r"ZENITH_deg:\s*([\d.]+)", content) + clouds_match = re.search(r"CLOUDS:\s*(\d+)", content) + + if not zenith_match or not clouds_match: + logger.warning("Required data fields missing (ZENITH and/or CLOUDS). \nSkipping file.") + return None + + zenith = float(zenith_match.group(1)) + num_clouds = int(clouds_match.group(1)) + + if num_clouds == 0: + logger.info("No cloud layers found in the report. \nSkipping file.") + return None + + layers = extract_cloud_layers(content, num_clouds) + timestamp = extract_timestamp_from_filename(input_filename) + + if timestamp is None: + return None + + extracted_data = { + "timestamp": timestamp.strftime("%Y-%m-%d %H:%M:%S"), + "lidar_zenith": zenith, + "layers": layers + } + + return extracted_data + + except Exception as e: + logger.error(f"Error processing file '{input_filename}': {e}") + return None + +def extract_cloud_layers(content, num_clouds): + """ + Extracts cloud layer data from the input file containig LIDAR quate reports. + + Parameters + ---------- + content : str + The content of the LIDAR report file. + num_clouds : int + The number of cloud layers. + + Returns + ------- + list + A list of dictionaries containing cloud layer parameters. + """ + + layers = [] + for i in range(num_clouds): + cloud_pattern = rf"CLOUD{i}\s*:\s*MEAN_HEIGHT:\s*([\d.]+)\s*\+\-\s*[\d.]+\s*BASE:\s*([\d.]+)\s*TOP:\s*([\d.]+)\s*FWHM:\s*[\d.]+\s*TRANS:\s*([\d.]+)" + cloud_match = re.search(cloud_pattern, content) + if cloud_match: + mean_height, base, top, trans = cloud_match.groups() + layers.append({ + "base_height": float(base), + "top_height": float(top), + "transmission": float(trans) + }) + else: + logger.warning(f"Cloud data not found.") + return layers + +def extract_timestamp_from_filename(input_filename): + """ + Extracts the timestamp from the input filename. + + Parameters + ---------- + input_filename : str + The name of the input file. + + Returns + ------- + datetime or None + The extracted timestamp as a datetime object, or None if the format is incorrect. + """ + + # Extract only the base name for filename validation + basename = os.path.basename(input_filename) + pattern = r"^(magic_lidar\.)?(\d{8})\.(\d{6})\.results$" + match = re.match(pattern, basename) + + if not match: + logger.error(f"Unexpected filename format for timestamp in file: '{input_filename}'") + logger.error("Expected format: 'magic_lidar.YYYYMMDD.HHMMSS.results'.") + return None + + date_str = match.group(2) # YYYYMMDD + time_str = match.group(3) # HHMMSS + + timestamp_str = f"{date_str}.{time_str}" + return datetime.strptime(timestamp_str, "%Y%m%d.%H%M%S") + + +def write_to_yaml(data, output_filename): + """ + Writes extracted cloud parameters to a YAML file, including units metadata. + + Parameters + ---------- + data : list + A list of dictionaries, each containing cloud data for a specific timestamp. + output_filename : str + The name of the output YAML file. + + Returns + ------- + None + """ + if not data: + logger.info("No valid data found to write.") + return + + output_data = { + "units": { + "timestamp": SQS("iso"), + "lidar_zenith": SQS("deg"), + "base_height": SQS("m"), + "top_height": SQS("m"), + "transmission": SQS("dimensionless"), + }, + "data": data + } + + yaml = YAML() + yaml.default_flow_style = False + yaml.sort_keys = False + + try: + with open(output_filename, 'w') as yaml_file: + yaml.dump(output_data, yaml_file) + logger.info(f"\nData written to output file '{output_filename}'.") + except IOError as e: + logger.error(f"Failed to write data to '{output_filename}': {e}") + +def main(): + """Main function.""" + parser = argparse.ArgumentParser() + + parser.add_argument( + "-i", + "--input-file", + type=str, + required=True, + help="Input file name or full path to input file. Wildcards allowed." + ) + + parser.add_argument( + "-o", + "--output-file", + type=str, + required=True, + help="Output YAML file name. One may provide full path where the output file should be created." + ) + + args = parser.parse_args() + results = [] + + files = glob.glob(args.input_file) + if not files: + logger.error(f"No file matching pattern '{args.input_file}'. Exiting.") + return + + for filepath in files: + basename = os.path.basename(filepath) + if not re.match(r"^(magic_lidar\.)?(\d{8})\.(\d{6})\.results$", basename): + logger.error(f"Invalid filename format: '{basename}'. Expected format:" + f"'magic_lidar.YYYYMMDD.HHMMSS.results'.") + continue + + data = extract_data_from_file(filepath) + if data: + results.append(data) + + if results: + results = sorted(results, key=lambda x: x["timestamp"]) + write_to_yaml(results, args.output_file) + logger.info("All data processed and written successfully.") + else: + logger.info("No valid LIDAR data extracted from input files.") + +if __name__ == "__main__": + main() From 9c1483956cde93372c744524fa29a1d4df44a76c Mon Sep 17 00:00:00 2001 From: Mario Pecimotika <45829235+mpecimotika@users.noreply.github.com> Date: Sun, 9 Feb 2025 16:29:50 +0100 Subject: [PATCH 5/8] fixing issues --- .../lst1_magic_lidar_reports_to_yaml.py | 99 +++++++++++-------- 1 file changed, 59 insertions(+), 40 deletions(-) diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_lidar_reports_to_yaml.py b/magicctapipe/scripts/lst1_magic/lst1_magic_lidar_reports_to_yaml.py index 4a90add9..84de9d2c 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_lidar_reports_to_yaml.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_lidar_reports_to_yaml.py @@ -2,7 +2,7 @@ # coding: utf-8 """ -This script processes LIDAR report files obtained with MARS-based software quate to extract cloud-reveant parameters present during the observations with MAGIC telescopes (if any), which can then be used to correct event images affected by the clouds using lst1_magic_cloud_correction.py script. +This script processes LIDAR report files obtained with MARS-based software quate to extract cloud-reveant parameters present during the observations with MAGIC telescopes (if any), which can then be used to correct event images affected by the clouds using lst1_magic_cloud_correction.py script. Usage: $ python lst1_magic_lidar_report_to_yaml.py @@ -10,49 +10,57 @@ --output-file "magic_lidar_report_clouds_CrabNebula.yaml" """ +import argparse +import glob +import logging import os import re -import glob -import argparse from datetime import datetime + from ruamel.yaml import YAML from ruamel.yaml.scalarstring import SingleQuotedScalarString as SQS -import logging logger = logging.getLogger(__name__) logger.addHandler(logging.StreamHandler()) logger.setLevel(logging.INFO) + def extract_data_from_file(input_filename): """ - Extracts cloud-relevant data from a LIDAR report file (if any). - + Extracts cloud-relevant data from a LIDAR report file (if any). + Parameters - ---------- + ----------- input_filename : str The name of an input LIDAR report file to process. Returns - ------- + -------- dict or None A dictionary with extracted data if valid cloud data is found, or None otherwise. """ - + try: with open(input_filename, "r") as file: content = file.read() - logger.info(f"\n----------------------------------------------------" - f"\nProcessing file: {input_filename}") + logger.info( + f"\n----------------------------------------------------" + f"\nProcessing file: {input_filename}" + ) if not re.search(r"ERROR_CODE:\s*4", content): - logger.info("No valid LIDAR report found (ERROR_CODE != 4). \nSkipping file.") + logger.info( + "No valid LIDAR report found (ERROR_CODE != 4). \nSkipping file." + ) return None zenith_match = re.search(r"ZENITH_deg:\s*([\d.]+)", content) clouds_match = re.search(r"CLOUDS:\s*(\d+)", content) if not zenith_match or not clouds_match: - logger.warning("Required data fields missing (ZENITH and/or CLOUDS). \nSkipping file.") + logger.warning( + "Required data fields missing (ZENITH and/or CLOUDS). \nSkipping file." + ) return None zenith = float(zenith_match.group(1)) @@ -71,58 +79,62 @@ def extract_data_from_file(input_filename): extracted_data = { "timestamp": timestamp.strftime("%Y-%m-%d %H:%M:%S"), "lidar_zenith": zenith, - "layers": layers + "layers": layers, } - + return extracted_data - + except Exception as e: logger.error(f"Error processing file '{input_filename}': {e}") return None + def extract_cloud_layers(content, num_clouds): """ Extracts cloud layer data from the input file containig LIDAR quate reports. Parameters - ---------- + ----------- content : str The content of the LIDAR report file. num_clouds : int The number of cloud layers. Returns - ------- + -------- list A list of dictionaries containing cloud layer parameters. """ - + layers = [] for i in range(num_clouds): cloud_pattern = rf"CLOUD{i}\s*:\s*MEAN_HEIGHT:\s*([\d.]+)\s*\+\-\s*[\d.]+\s*BASE:\s*([\d.]+)\s*TOP:\s*([\d.]+)\s*FWHM:\s*[\d.]+\s*TRANS:\s*([\d.]+)" cloud_match = re.search(cloud_pattern, content) if cloud_match: mean_height, base, top, trans = cloud_match.groups() - layers.append({ - "base_height": float(base), - "top_height": float(top), - "transmission": float(trans) - }) + layers.append( + { + "base_height": float(base), + "top_height": float(top), + "transmission": float(trans), + } + ) else: - logger.warning(f"Cloud data not found.") + logger.warning("Cloud data not found.") return layers + def extract_timestamp_from_filename(input_filename): """ Extracts the timestamp from the input filename. Parameters - ---------- + ----------- input_filename : str The name of the input file. Returns - ------- + -------- datetime or None The extracted timestamp as a datetime object, or None if the format is incorrect. """ @@ -133,7 +145,9 @@ def extract_timestamp_from_filename(input_filename): match = re.match(pattern, basename) if not match: - logger.error(f"Unexpected filename format for timestamp in file: '{input_filename}'") + logger.error( + f"Unexpected filename format for timestamp in file: '{input_filename}'" + ) logger.error("Expected format: 'magic_lidar.YYYYMMDD.HHMMSS.results'.") return None @@ -142,22 +156,23 @@ def extract_timestamp_from_filename(input_filename): timestamp_str = f"{date_str}.{time_str}" return datetime.strptime(timestamp_str, "%Y%m%d.%H%M%S") - + def write_to_yaml(data, output_filename): """ Writes extracted cloud parameters to a YAML file, including units metadata. Parameters - ---------- + ----------- data : list A list of dictionaries, each containing cloud data for a specific timestamp. output_filename : str The name of the output YAML file. - Returns + Return ------- None + This function returs nothing. """ if not data: logger.info("No valid data found to write.") @@ -171,7 +186,7 @@ def write_to_yaml(data, output_filename): "top_height": SQS("m"), "transmission": SQS("dimensionless"), }, - "data": data + "data": data, } yaml = YAML() @@ -179,22 +194,23 @@ def write_to_yaml(data, output_filename): yaml.sort_keys = False try: - with open(output_filename, 'w') as yaml_file: + with open(output_filename, "w") as yaml_file: yaml.dump(output_data, yaml_file) logger.info(f"\nData written to output file '{output_filename}'.") except IOError as e: logger.error(f"Failed to write data to '{output_filename}': {e}") + def main(): """Main function.""" parser = argparse.ArgumentParser() - + parser.add_argument( "-i", "--input-file", type=str, required=True, - help="Input file name or full path to input file. Wildcards allowed." + help="Input file name or full path to input file. Wildcards allowed.", ) parser.add_argument( @@ -202,9 +218,9 @@ def main(): "--output-file", type=str, required=True, - help="Output YAML file name. One may provide full path where the output file should be created." + help="Output YAML file name. One may provide full path where the output file should be created.", ) - + args = parser.parse_args() results = [] @@ -216,10 +232,12 @@ def main(): for filepath in files: basename = os.path.basename(filepath) if not re.match(r"^(magic_lidar\.)?(\d{8})\.(\d{6})\.results$", basename): - logger.error(f"Invalid filename format: '{basename}'. Expected format:" - f"'magic_lidar.YYYYMMDD.HHMMSS.results'.") + logger.error( + f"Invalid filename format: '{basename}'. Expected format:" + f"'magic_lidar.YYYYMMDD.HHMMSS.results'." + ) continue - + data = extract_data_from_file(filepath) if data: results.append(data) @@ -231,5 +249,6 @@ def main(): else: logger.info("No valid LIDAR data extracted from input files.") + if __name__ == "__main__": main() From 15fffd60e6b0e2cd0060676824ad41cff73842eb Mon Sep 17 00:00:00 2001 From: nzywucka Date: Mon, 17 Feb 2025 14:12:00 +0100 Subject: [PATCH 6/8] The additional cleaning part is updated. --- .../lst1_magic/lst_m1_m2_cloud_correction.py | 40 ++++++++++++++----- 1 file changed, 29 insertions(+), 11 deletions(-) diff --git a/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py b/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py index 35de2e7a..ff8ba738 100644 --- a/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py +++ b/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py @@ -26,6 +26,7 @@ from astropy.time import Time from ctapipe.coordinates import TelescopeFrame from ctapipe.image import ( + apply_time_delta_cleaning, concentration_parameters, hillas_parameters, leakage_parameters, @@ -34,6 +35,7 @@ ) from ctapipe.instrument import SubarrayDescription from ctapipe.io import read_table +from lstchain.image.cleaning import apply_dynamic_cleaning from scipy.interpolate import interp1d import magicctapipe @@ -121,18 +123,18 @@ def lidar_cloud_interpolation( mean_subrun_timestamp, max_gap_lidar_shots, lidar_report_file ): """ - Retrieves or interpolates LIDAR cloud parameters based on the closest timestamps to an input mean timestamp of the processed subrun. + Retrieves or interpolates LIDAR cloud parameters based on the closest timestamps to an input mean timestamp of the processed subrun. For the moment in the case of multiple clouds only the one with the lowest transmission is taken into account. Parameters ----------- - mean_subrun_timestamp : int or float + mean_subrun_timestamp : int The mean timestamp of the processed subrun (format: unix). - max_gap_lidar_shots : int or float + max_gap_lidar_shots : float Maximum allowed time gap for interpolation (in seconds). lidar_report_file : str - Path to the yaml file containing LIDAR laser reports with columns: + Path to the yaml file created by lst1_magic_lidar_reports_to_yaml.py containing LIDAR laser reports with columns: - "timestamp" (format: ISO 8601) - "base_height", "top_height", "transmission", "lidar_zenith" @@ -356,9 +358,20 @@ def process_telescope_data( """ # AC LST + cleaning_level_lst = config["LST"]["tailcuts_clean"] modified_boundary_thresh_lst = cleaning_level_lst["boundary_thresh"] * cmf modified_picture_thresh_lst = cleaning_level_lst["picture_thresh"] * cmf + keep_isolated_pixels = cleaning_level_lst["keep_isolated_pixels"] + + use_time_delta_cleaning = config["LST"]["time_delta_cleaning"] + min_number_neighbors = use_time_delta_cleaning["min_number_neighbors"] + time_limit = use_time_delta_cleaning["time_limit"] + + use_dynamic_cleaning = config["LST"]["dynamic_cleaning"] + threshold = use_dynamic_cleaning["threshold"] + fraction = use_dynamic_cleaning["fraction"] + # AC MAGIC # Configure the MAGIC image cleaning config_clean_magic = config["MAGIC"]["magic_clean"] @@ -466,23 +479,28 @@ def process_telescope_data( image, dtype=bool ) # Assuming full mask if not defined if tel_ids["LST-1"] == tel_id: - clean = tailcuts_clean( + clean_mask = tailcuts_clean( camgeom, image, boundary_thresh=modified_boundary_thresh_lst, picture_thresh=modified_picture_thresh_lst, + keep_isolated_pixels=keep_isolated_pixels, min_number_picture_neighbors=cleaning_level_lst[ "min_number_picture_neighbors" ], ) - # Create a mask indicating which pixels survived the cleaning - clean_mask = np.zeros_like(image, dtype=bool) - clean_mask[clean] = True - if np.sum(clean_mask) == 0: - continue + clean_mask = apply_time_delta_cleaning( + camgeom, + clean_mask, + peak_time, + min_number_neighbors, + time_limit, + ) - # Apply the mask to relevant data arrays + clean_mask = apply_dynamic_cleaning( + image, clean_mask, threshold, fraction + ) clean_peak_time = peak_time[clean_mask] clean_image = image[clean_mask] clean_camgeom = camgeom[clean_mask] From a7d6e4aa3daa9f8cb9df1dbcc8d75ec49ac54fa2 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Mon, 17 Feb 2025 14:23:16 +0100 Subject: [PATCH 7/8] The additional cleaning of LST is updated. --- magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py | 1 - 1 file changed, 1 deletion(-) diff --git a/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py b/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py index 63071696..78b105fb 100644 --- a/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py +++ b/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py @@ -360,7 +360,6 @@ def process_telescope_data( """ # AC LST - cleaning_level_lst = config["LST"]["tailcuts_clean"] modified_boundary_thresh_lst = cleaning_level_lst["boundary_thresh"] * cmf modified_picture_thresh_lst = cleaning_level_lst["picture_thresh"] * cmf From 0ac1c8de862ce027e83a4786b261ca2bad82acc3 Mon Sep 17 00:00:00 2001 From: nzywucka Date: Mon, 17 Feb 2025 20:15:14 +0100 Subject: [PATCH 8/8] Another update of additional cleaning of LST-1 part. --- .../scripts/lst1_magic/lst_m1_m2_cloud_correction.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py b/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py index 78b105fb..4cc6b72c 100644 --- a/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py +++ b/magicctapipe/scripts/lst1_magic/lst_m1_m2_cloud_correction.py @@ -470,7 +470,7 @@ def process_telescope_data( raise ValueError("Error: 'inds_img' list is empty!") index_img = inds_img[0] image = dl1_images["image"][index_img] - clean_mask = dl1_images["image_mask"][index_img] + clean_mask_file = dl1_images["image_mask"][index_img] peak_time = dl1_images["peak_time"][index_img] image /= trans_pixels @@ -502,6 +502,12 @@ def process_telescope_data( clean_mask = apply_dynamic_cleaning( image, clean_mask, threshold, fraction ) + + clean_mask = clean_mask * clean_mask_file + + if np.sum(clean_mask) == 0: + continue + clean_peak_time = peak_time[clean_mask] clean_image = image[clean_mask] clean_camgeom = camgeom[clean_mask]