From cf8db3379087234dd9819a80073af2dd11975cda Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Thu, 28 May 2026 09:52:53 +0200 Subject: [PATCH 01/10] grid density and plotting --- .../applications/plot_production_grid.py | 44 +-- .../applications/production_generate_grid.py | 14 + .../production_configuration/job_grid_io.py | 32 +- .../observation_grid.py | 108 ++++++- .../simulation_jobs.py | 78 +++++ .../visualization/plot_production_grid.py | 268 +++++----------- .../test_production_generate_grid.py | 9 + .../test_job_grid_io.py | 4 + .../test_observation_grid.py | 59 +++- .../test_simulation_jobs.py | 153 +++++++++ .../test_plot_production_grid.py | 291 +++++------------- 11 files changed, 608 insertions(+), 452 deletions(-) diff --git a/src/simtools/applications/plot_production_grid.py b/src/simtools/applications/plot_production_grid.py index 73165b4f65..763d395c41 100644 --- a/src/simtools/applications/plot_production_grid.py +++ b/src/simtools/applications/plot_production_grid.py @@ -11,14 +11,6 @@ ---------------------- grid_points_file (str, required) Path to the ECSV file containing grid points. -site (str, required) - Observatory site name used to read the reference coordinates from the site model. -model_version (str, required) - Model version used to read the site reference coordinates. -observation_time (str, optional) - Observation time in UTC ISO format used for Alt/Az <-> RA/Dec transformations. - If omitted, the application uses ``metadata.time_of_observation_utc`` from the - grid file when available. plot_ra_dec_tracks (flag, optional) If provided, plot RA/Dec guide tracks on top of the sky projection. When native RA/Dec grid points are present (grid file contains explicit ``ra`` and ``dec`` @@ -37,25 +29,17 @@ simtools-plot-production-grid \ --grid_points_file path/to/grid_points_production.ecsv \ - --site North \ - --model_version 6.0.2 \ - --observation_time "2025-06-01 00:00:00" \ --plot_ra_dec_tracks Output ------ -The output figure shows both local Alt/Az (polar projection) and equatorial -RA/Dec panels. +The output figure shows local Alt/Az (polar projection). The equatorial +RA/Dec panel is added when RA/Dec columns are available in the grid file. """ -import logging - from simtools.application_control import build_application -from simtools.model.site_model import SiteModel from simtools.visualization.plot_production_grid import ProductionGridPlotter -logger = logging.getLogger(__name__) - def _add_arguments(parser): """Register application-specific command line arguments.""" @@ -65,15 +49,6 @@ def _add_arguments(parser): required=True, help="Path to the ECSV file containing grid points.", ) - parser.add_argument( - "--observation_time", - type=str, - default=None, - help=( - "Observation time in UTC ISO format for coordinate transforms. " - "If not provided, uses observing time stored in the grid file metadata when present." - ), - ) parser.add_argument( "--plot_ra_dec_tracks", action="store_true", @@ -91,23 +66,10 @@ def _add_arguments(parser): def main(): """Run the ProductionGridPlotter.""" - app_context = build_application( - initialization_kwargs={ - "db_config": True, - "simulation_model": ["version", "site", "model_version"], - } - ) - site_model = SiteModel( - model_version=app_context.args["model_version"], - site=app_context.args["site"], - ) + app_context = build_application(initialization_kwargs={"db_config": False, "output": True}) plotter = ProductionGridPlotter( grid_points_file=app_context.args["grid_points_file"], - site_location_lat=site_model.get_parameter_value_with_unit("reference_point_latitude"), - site_location_lon=site_model.get_parameter_value_with_unit("reference_point_longitude"), - site_location_height=site_model.get_parameter_value_with_unit("reference_point_altitude"), - observation_time=app_context.args["observation_time"], output_path=app_context.io_handler.get_output_directory(), ) diff --git a/src/simtools/applications/production_generate_grid.py b/src/simtools/applications/production_generate_grid.py index 388362ce22..8bc6b50d74 100644 --- a/src/simtools/applications/production_generate_grid.py +++ b/src/simtools/applications/production_generate_grid.py @@ -22,6 +22,10 @@ ``--axis [scaling]``. Example: ``--axis azimuth 310 deg 20 deg 3 linear``. Options for scaling are: linear, log, 1/cos. +direction_grid_density (float, optional) + Direction-grid density in points per deg^2. + If provided, direction-axis binning (azimuth/zenith or ra/dec) is derived from + range and density, while min/max values are kept from ``--axis`` definitions. time_of_observation (str, optional) Time of the observation in UTC (format: 'YYYY-MM-DD HH:MM:SS'). Used only if RA/Dec axes are provided (for coordinate transforms and sidereal-time @@ -97,6 +101,16 @@ def _add_arguments(parser): "Observation time in UTC (format: 'YYYY-MM-DD HH:MM:SS'). Used only in 'ra_dec' mode." ), ) + parser.add_argument( + "--direction_grid_density", + type=float, + required=False, + default=None, + help=( + "Direction-grid density in points per deg^2. If set, direction-axis binning is " + "derived from axis ranges and this density." + ), + ) parser.add_argument( "--output_file", type=str, diff --git a/src/simtools/production_configuration/job_grid_io.py b/src/simtools/production_configuration/job_grid_io.py index be1bf5c7ea..f7978ff868 100644 --- a/src/simtools/production_configuration/job_grid_io.py +++ b/src/simtools/production_configuration/job_grid_io.py @@ -3,6 +3,7 @@ import logging from pathlib import Path +import numpy as np from astropy import units as u from astropy.table import Table @@ -43,6 +44,8 @@ "view_cone_max": ("view_cone_max_value", "view_cone_max_unit"), } +_OPTIONAL_ANGLE_FIELDS = ("ra", "dec") + def _serialize_quantity(value): """Serialize a Quantity to value/unit columns.""" @@ -72,6 +75,15 @@ def _serialize_job_row(job_row): job_row[quantity_name] ) + for angle_name in _OPTIONAL_ANGLE_FIELDS: + angle_value = job_row.get(angle_name) + if angle_value is None: + continue + if isinstance(angle_value, u.Quantity): + serialized_row[angle_name] = float(angle_value.to_value(u.deg)) + else: + serialized_row[angle_name] = float(angle_value) + return serialized_row @@ -94,6 +106,14 @@ def _deserialize_job_row(serialized_row): serialized_row[unit_key], ) + for angle_name in _OPTIONAL_ANGLE_FIELDS: + if angle_name not in serialized_row: + continue + angle_value = serialized_row[angle_name] + if np.ma.is_masked(angle_value) or angle_value is None: + continue + job_row[angle_name] = float(angle_value) * u.deg + return job_row @@ -117,7 +137,17 @@ def serialize_job_grid(job_rows, output_file, metadata=None): if output_path.suffix.lower() != ".ecsv": raise ValueError("Job grid output file must use the '.ecsv' extension.") - output_table = Table(rows=serialized_rows, names=JOB_GRID_COLUMNS) + optional_columns = [ + angle_name + for angle_name in _OPTIONAL_ANGLE_FIELDS + if any(angle_name in row for row in serialized_rows) + ] + output_columns = [*JOB_GRID_COLUMNS, *optional_columns] + output_rows = [ + {column: row.get(column) for column in output_columns} for row in serialized_rows + ] + + output_table = Table(rows=output_rows, names=output_columns) output_table.meta = metadata logger.info(f"Writing job grid with {len(job_rows)} rows to '{output_path}'.") output_table.write(output_path, format="ascii.ecsv", overwrite=True) diff --git a/src/simtools/production_configuration/observation_grid.py b/src/simtools/production_configuration/observation_grid.py index 6c343a0359..dadd25166c 100644 --- a/src/simtools/production_configuration/observation_grid.py +++ b/src/simtools/production_configuration/observation_grid.py @@ -299,27 +299,111 @@ def create_circular_binning(self, azimuth_range, num_bins): np.ndarray Array of bin centers. """ - azimuth_min, azimuth_max = azimuth_range - azimuth_min %= 360 - azimuth_max %= 360 + raw_span = abs(float(azimuth_range[1]) - float(azimuth_range[0])) + if raw_span > 0.0 and np.isclose(raw_span % 360.0, 0.0): + azimuth_start = float(azimuth_range[0]) % 360.0 + return (azimuth_start + np.linspace(0.0, 360.0, num_bins, endpoint=False)) % 360.0 - clockwise_distance = (azimuth_max - azimuth_min) % 360 - counterclockwise_distance = (azimuth_min - azimuth_max) % 360 + azimuth_min, azimuth_max = azimuth_range + azimuth_start = float(azimuth_min) % 360.0 + directed_distance = float((azimuth_max - azimuth_min) % 360.0) - if clockwise_distance <= counterclockwise_distance: - return ( - np.linspace(azimuth_min, azimuth_min + clockwise_distance, num_bins, endpoint=True) - % 360 - ) return ( np.linspace( - azimuth_min, azimuth_min - counterclockwise_distance, num_bins, endpoint=True + azimuth_start, + azimuth_start + directed_distance, + num_bins, + endpoint=True, ) - % 360 + % 360.0 + ) + + @staticmethod + def _ceil_with_tolerance(value): + """Ceil a float while avoiding near-integer floating-point artifacts.""" + nearest_integer = round(value) + if np.isclose(value, nearest_integer): + return int(nearest_integer) + return int(np.ceil(value)) + + @staticmethod + def _directed_circular_span_degrees(azimuth_range): + """Return directed circular span (degrees) from start to end.""" + raw_span = abs(float(azimuth_range[1]) - float(azimuth_range[0])) + if raw_span > 0.0 and np.isclose(raw_span % 360.0, 0.0): + return 360.0 + return float((azimuth_range[1] - azimuth_range[0]) % 360.0) + + def _is_adaptive_horizontal_density_enabled(self): + """Return whether horizontal grid generation should use row-wise adaptive azimuth bins.""" + return ( + self.coordinate_system == "horizontal" + and "azimuth" in self.axes + and "zenith_angle" in self.axes + and self.axes["azimuth"].get("direction_grid_density") is not None ) + def _generate_adaptive_horizontal_grid(self): + """Generate horizontal grid with azimuth binning adapted per zenith row.""" + density = float(self.axes["azimuth"]["direction_grid_density"]) + azimuth_range = self.axes["azimuth"]["range"] + azimuth_span = self._directed_circular_span_degrees(azimuth_range) + zenith_values = self.target_values["zenith_angle"] + + zenith_step = 1.0 / np.sqrt(density) + if len(zenith_values) > 1: + zenith_step = float(np.mean(np.abs(np.diff(zenith_values.to_value(u.deg))))) + + extra_axis_keys = [ + key for key in self.target_values if key not in ("azimuth", "zenith_angle") + ] + extra_arrays = [self.target_values[key].value for key in extra_axis_keys] + extra_units = [self.target_values[key].unit for key in extra_axis_keys] + + if extra_arrays: + extra_grid = np.meshgrid(*extra_arrays, indexing="ij") + extra_combinations = np.vstack(list(map(np.ravel, extra_grid))).T + else: + extra_combinations = [np.array([])] + + grid_points = [] + for zenith in zenith_values: + altitude_cosine = np.cos(np.deg2rad(90.0 - zenith.to_value(u.deg))) + azimuth_bins = 1 + if azimuth_span > 0.0 and altitude_cosine > 0.0: + azimuth_bins = max( + 1, + self._ceil_with_tolerance( + azimuth_span * density * zenith_step * altitude_cosine + ), + ) + + azimuth_values = self.create_circular_binning(azimuth_range, azimuth_bins) * u.deg + for extra_combination in extra_combinations: + point_base = { + key: Quantity(extra_combination[i], extra_units[i]) + for i, key in enumerate(extra_axis_keys) + } + for azimuth in azimuth_values: + point = { + **point_base, + "azimuth": azimuth, + "zenith_angle": zenith, + } + self._add_lookup_limits_to_point( + point, + zenith=zenith.to_value(u.deg), + azimuth=azimuth.to_value(u.deg), + ) + grid_points.append(point) + + return grid_points + def _generate_horizontal_grid(self): """Generate grid points for horizontal (zenith/azimuth) mode.""" + if self._is_adaptive_horizontal_density_enabled(): + return self._generate_adaptive_horizontal_grid() + value_arrays = [value.value for value in self.target_values.values()] units = [value.unit for value in self.target_values.values()] grid = np.meshgrid(*value_arrays, indexing="ij") diff --git a/src/simtools/production_configuration/simulation_jobs.py b/src/simtools/production_configuration/simulation_jobs.py index d770729bec..1a3d5f84d1 100644 --- a/src/simtools/production_configuration/simulation_jobs.py +++ b/src/simtools/production_configuration/simulation_jobs.py @@ -174,6 +174,71 @@ def _resolve_axis_configs(args_dict): return axis_configs +def _circular_span_degrees(axis_range): + """Return directed circular span (degrees) from start to end.""" + start, end = axis_range + raw_span = abs(float(end) - float(start)) + if raw_span > 0.0 and np.isclose(raw_span % 360.0, 0.0): + return 360.0 + return float((end - start) % 360.0) + + +def _ceil_with_tolerance(value): + """Ceil a float while avoiding near-integer floating-point artifacts.""" + nearest_integer = round(value) + if np.isclose(value, nearest_integer): + return int(nearest_integer) + return int(np.ceil(value)) + + +def _mean_cosine_over_dec_span(dec_range): + """Return mean cos(dec) over a declination interval (degrees).""" + dec_min, dec_max = sorted((float(dec_range[0]), float(dec_range[1]))) + dec_span_rad = np.deg2rad(abs(dec_max - dec_min)) + if np.isclose(dec_span_rad, 0.0): + return abs(np.cos(np.deg2rad(dec_min))) + return abs((np.sin(np.deg2rad(dec_max)) - np.sin(np.deg2rad(dec_min))) / dec_span_rad) + + +def _mean_sine_over_zenith_span(zenith_range): + """Return mean sin(zenith) over a zenith interval (degrees).""" + zenith_min, zenith_max = sorted((float(zenith_range[0]), float(zenith_range[1]))) + zenith_span_rad = np.deg2rad(abs(zenith_max - zenith_min)) + if np.isclose(zenith_span_rad, 0.0): + return abs(np.sin(np.deg2rad(zenith_min))) + return abs((np.cos(np.deg2rad(zenith_min)) - np.cos(np.deg2rad(zenith_max))) / zenith_span_rad) + + +def _apply_direction_grid_density(axis_configs, direction_axes, density): + """Derive direction-axis binning from density (points per deg^2).""" + if density is None: + return + if density <= 0: + raise ValueError("direction_grid_density must be strictly positive.") + + density_sqrt = np.sqrt(density) + longitudinal_axis_scale = 1.0 + if tuple(direction_axes) == _RADEC_AXES: + longitudinal_axis_scale = _mean_cosine_over_dec_span(axis_configs["dec"]["range"]) + elif tuple(direction_axes) == _HORIZONTAL_AXES: + longitudinal_axis_scale = _mean_sine_over_zenith_span(axis_configs["zenith"]["range"]) + + for axis_name in direction_axes: + axis_range = axis_configs[axis_name]["range"] + if axis_name == "azimuth": + span_degrees = _circular_span_degrees(axis_range) + else: + span_degrees = abs(axis_range[1] - axis_range[0]) + + if axis_name in ("azimuth", "ra"): + span_degrees *= longitudinal_axis_scale + + axis_configs[axis_name]["binning"] = max( + 1, + _ceil_with_tolerance(span_degrees * density_sqrt), + ) + + def _resolve_coordinate_system(axis_configs): """Resolve the coordinate system from axis definitions.""" has_horizontal_axes = all(axis_name in axis_configs for axis_name in _HORIZONTAL_AXES) @@ -236,6 +301,15 @@ def build_axes_dict_from_cli_args(args_dict): raise ValueError(f"Missing required shared axis definition(s): {missing_axes}.") direction_axes = _RADEC_AXES if coordinate_system == "ra_dec" else _HORIZONTAL_AXES + if coordinate_system == "horizontal" and args_dict.get("direction_grid_density") is not None: + axis_configs["azimuth"]["direction_grid_density"] = float( + args_dict["direction_grid_density"] + ) + _apply_direction_grid_density( + axis_configs, + direction_axes, + args_dict.get("direction_grid_density"), + ) return { GRID_AXIS_ARGUMENTS[axis_name]["engine_axis"]: axis_configs[axis_name] for axis_name in (*_REQUIRED_AXES, *direction_axes) @@ -280,6 +354,8 @@ def build_job_grid_metadata(args_dict): args_dict, ) coordinate_system = _resolve_coordinate_system_from_args(args_dict) + if coordinate_system == "ra_dec" and not args_dict.get("site"): + raise ValueError("site is required when using RA/Dec axes.") return { "site": args_dict.get("site"), "simulation_software": args_dict.get("simulation_software"), @@ -692,6 +768,8 @@ def build_simulation_jobs(args_dict): "primary": primary, "azimuth_angle": point["azimuth"], "zenith_angle": point["zenith_angle"], + "ra": point.get("ra"), + "dec": point.get("dec"), "model_version": model_version, "array_layout_name": resolved_layout_name, "corsika_le_interaction": corsika_le, diff --git a/src/simtools/visualization/plot_production_grid.py b/src/simtools/visualization/plot_production_grid.py index a138ca77b0..8075d37785 100644 --- a/src/simtools/visualization/plot_production_grid.py +++ b/src/simtools/visualization/plot_production_grid.py @@ -3,12 +3,9 @@ import logging from pathlib import Path -import astropy.units as u import matplotlib.pyplot as plt import numpy as np -from astropy.coordinates import AltAz, EarthLocation, SkyCoord from astropy.table import Table -from astropy.time import Time logger = logging.getLogger(__name__) DEFAULT_OUTPUT_FILE_STEM = "production_grid_sky_projection" @@ -26,17 +23,6 @@ class ProductionGridPlotter: ---------- grid_points_file : str or Path Path to the ECSV file containing grid points. - site_location_lat : float or astropy.units.Quantity - Site latitude in degrees. - site_location_lon : float or astropy.units.Quantity - Site longitude in degrees. - site_location_height : float or astropy.units.Quantity - Site height in meters. - observation_time : str, optional - Observation time in UTC ISO format used only for coordinate transformations - during plotting (e.g. RA/Dec <-> Alt/Az and horizon visibility). It does - not modify or re-generate the underlying grid points from the input file. - If None, metadata from the grid file is used. output_path : str or Path Path to save output plots. """ @@ -44,10 +30,6 @@ class ProductionGridPlotter: def __init__( self, grid_points_file, - site_location_lat, - site_location_lon, - site_location_height, - observation_time, output_path, ): """Initialize the ProductionGridPlotter.""" @@ -55,29 +37,12 @@ def __init__( self.output_path = Path(output_path) self.output_path.mkdir(parents=True, exist_ok=True) - latitude = u.Quantity(site_location_lat, unit=u.deg) - longitude = u.Quantity(site_location_lon, unit=u.deg) - height = u.Quantity(site_location_height, unit=u.m) - - self.location = EarthLocation( - lat=latitude, - lon=longitude, - height=height, - ) self.grid_metadata = {} + self.grid_columns = [] + self.has_radec_columns = False self.grid_points = self._load_grid_points() - observation_time_utc = observation_time or self.grid_metadata.get("time_of_observation_utc") - if observation_time_utc is None: - observation_time_utc = "2025-01-01 00:00:00" - self.time = Time(observation_time_utc, scale="utc") - logger.info(f"Loaded {len(self.grid_points)} grid points from {self.grid_points_file}") - logger.info( - f"Site location: lat={latitude.to_value(u.deg)} deg, " - f"lon={longitude.to_value(u.deg)} deg" - ) - logger.info(f"Observation time (UTC): {self.time.isot}") def _load_grid_points(self): """ @@ -106,6 +71,11 @@ def _load_grid_points(self): grid_table = Table.read(self.grid_points_file, format="ascii.ecsv") self.grid_metadata = dict(grid_table.meta) + self.grid_columns = list(grid_table.colnames) + self.has_radec_columns = {"ra", "dec"}.issubset(self.grid_columns) or { + "ra_value", + "dec_value", + }.issubset(self.grid_columns) return self._convert_ecsv_table_to_grid_points(grid_table) @staticmethod @@ -152,9 +122,18 @@ def _extract_quantity_value(point, key): return None return float(value) + @classmethod + def _extract_first_available_quantity_value(cls, point, keys): + """Extract first available quantity value from a list of candidate keys.""" + for key in keys: + value = cls._extract_quantity_value(point, key) + if value is not None: + return value + return None + def _normalize_grid_point(self, point): """ - Normalize a grid point to both Alt/Az and RA/Dec when possible. + Normalize a grid point from available columns. Parameters ---------- @@ -166,41 +145,37 @@ def _normalize_grid_point(self, point): dict or None Normalized point with native frame metadata. """ - azimuth = self._extract_quantity_value(point, "azimuth") - zenith = self._extract_quantity_value(point, "zenith_angle") - ra = self._extract_quantity_value(point, "ra") - dec = self._extract_quantity_value(point, "dec") + azimuth = self._extract_first_available_quantity_value( + point, + ("azimuth", "azimuth_angle", "azimuth_angle_value"), + ) + zenith = self._extract_first_available_quantity_value( + point, + ("zenith_angle", "zenith_angle_value"), + ) + ra = self._extract_first_available_quantity_value(point, ("ra", "ra_value")) + dec = self._extract_first_available_quantity_value(point, ("dec", "dec_value")) + point_ra = float(ra % 360.0) if ra is not None else None + point_dec = float(dec) if dec is not None else None if azimuth is not None and zenith is not None: - alt = (90.0 - zenith) * u.deg - skycoord = SkyCoord( - AltAz( - alt=alt, - az=(azimuth % 360.0) * u.deg, - obstime=self.time, - location=self.location, - ) - ) return { "native_frame": "altaz", "azimuth": azimuth % 360.0, "zenith": zenith, - "ra": float(skycoord.icrs.ra.deg % 360.0), - "dec": float(skycoord.icrs.dec.deg), + "ra": point_ra, + "dec": point_dec, "visible_in_altaz": True, } if ra is not None and dec is not None: - skycoord = SkyCoord(ra=ra * u.deg, dec=dec * u.deg, frame="icrs") - altaz = skycoord.transform_to(AltAz(obstime=self.time, location=self.location)) - visible_in_altaz = bool(altaz.alt.deg > 0.0) return { "native_frame": "radec", - "azimuth": float(altaz.az.deg % 360.0) if visible_in_altaz else None, - "zenith": float(90.0 - altaz.alt.deg) if visible_in_altaz else None, + "azimuth": None, + "zenith": None, "ra": float(ra % 360.0), "dec": float(dec), - "visible_in_altaz": visible_in_altaz, + "visible_in_altaz": False, } logger.warning(f"Skipping point without supported coordinates: {point}") @@ -285,32 +260,61 @@ def plot_sky_projection(self, plot_ra_dec_tracks=False, dec_values=None): dec_values : list of float, optional List of declination values to plot as tracks. """ - figure = plt.figure(figsize=(15, 7)) - altaz_axis = figure.add_subplot(1, 2, 1, projection="polar") - radec_axis = figure.add_subplot(1, 2, 2) - plot_points = self.normalize_grid_points() - altaz_count = self._plot_altaz_points(altaz_axis, plot_points) - radec_count = self._plot_radec_points(radec_axis, plot_points) - inferred_grid_count = 0 + has_radec_points = any( + point["ra"] is not None and point["dec"] is not None for point in plot_points + ) + show_radec_panel = self.has_radec_columns and has_radec_points - if plot_ra_dec_tracks and dec_values: - self._plot_ra_dec_tracks(altaz_axis, dec_values) - elif plot_ra_dec_tracks: - inferred_grid_count = self._plot_inferred_radec_grid(altaz_axis, plot_points) + figure = plt.figure(figsize=(15, 7) if show_radec_panel else (8, 7)) + if show_radec_panel: + altaz_axis = figure.add_subplot(1, 2, 1, projection="polar") + radec_axis = figure.add_subplot(1, 2, 2) + else: + altaz_axis = figure.add_subplot(1, 1, 1, projection="polar") + radec_axis = None - if altaz_count > 0 or (plot_ra_dec_tracks and (dec_values or inferred_grid_count > 0)): - altaz_axis.legend(loc="upper left", bbox_to_anchor=(1.0, 1.15)) + altaz_count = self._plot_altaz_points(altaz_axis, plot_points) + radec_count = self._plot_radec_points(radec_axis, plot_points) if radec_axis else 0 + if plot_ra_dec_tracks: + if dec_values: + logger.info( + "RA/Dec tracks are disabled in file-driven plotting mode " + "(ignoring manual dec_values)" + ) + else: + logger.info("RA/Dec tracks are disabled in file-driven plotting mode") + + if altaz_count > 0: + _, altaz_labels = altaz_axis.get_legend_handles_labels() + if any(label and not label.startswith("_") for label in altaz_labels): + altaz_axis.legend(loc="upper left", bbox_to_anchor=(1.0, 1.15)) + + if radec_axis and radec_count > 0: + _, radec_labels = radec_axis.get_legend_handles_labels() + if any(label and not label.startswith("_") for label in radec_labels): + radec_axis.legend(loc="upper right") + + figure.suptitle("Production Grid Points", fontsize=14, y=0.98) + subtitle_lines = [] + if self.grid_metadata.get("site"): + subtitle_lines.append(f"Site: {self.grid_metadata['site']}") + if self.grid_metadata.get("time_of_observation_utc"): + subtitle_lines.append( + f"Observation time (UTC): {self.grid_metadata['time_of_observation_utc']}" + ) - if radec_count > 0: - radec_axis.legend(loc="upper right") + for index, subtitle_line in enumerate(subtitle_lines): + figure.text( + 0.5, + 0.95 - index * 0.025, + subtitle_line, + ha="center", + va="top", + fontsize=10, + ) - figure.suptitle( - f"Production Grid Points\nObservation time: {self.time.iso}", - fontsize=14, - y=0.98, - ) - figure.tight_layout(rect=(0.0, 0.0, 1.0, 0.95)) + figure.tight_layout(rect=(0.0, 0.0, 1.0, 0.92)) output_file = self.output_path / ( f"{DEFAULT_OUTPUT_FILE_STEM}.{DEFAULT_OUTPUT_FILE_EXTENSION}" @@ -361,7 +365,7 @@ def _configure_radec_axis(self, axis, plot_points): axis.set_ylim(dec_min, dec_max) @staticmethod - def _scatter_group(axis, plot_points, label, color, x_key, y_key, x_transform=None): + def _scatter_group(axis, plot_points, color, x_key, y_key, x_transform=None): """Scatter a group of points with configurable coordinates.""" if not plot_points: return 0 @@ -378,7 +382,6 @@ def _scatter_group(axis, plot_points, label, color, x_key, y_key, x_transform=No linewidths=1.0, edgecolors=color, facecolors="none", - label=label, zorder=10, ) return len(plot_points) @@ -389,8 +392,6 @@ def _plot_frame_points( plot_points, primary_frame, secondary_frame, - primary_label, - secondary_label, primary_color, secondary_color, x_key, @@ -415,7 +416,6 @@ def _plot_frame_points( plotted_points += self._scatter_group( axis, primary_points, - label=primary_label, color=primary_color, x_key=x_key, y_key=y_key, @@ -424,7 +424,6 @@ def _plot_frame_points( plotted_points += self._scatter_group( axis, secondary_points, - label=secondary_label, color=secondary_color, x_key=x_key, y_key=y_key, @@ -455,8 +454,6 @@ def _plot_altaz_points(self, axis, plot_points): plot_points=plot_points, primary_frame="altaz", secondary_frame="radec", - primary_label="Native Alt/Az points", - secondary_label="RA/Dec transformed to Alt/Az", primary_color="tab:blue", secondary_color="tab:orange", x_key="azimuth", @@ -488,104 +485,9 @@ def _plot_radec_points(self, axis, plot_points): plot_points=plot_points, primary_frame="radec", secondary_frame="altaz", - primary_label="Native RA/Dec points", - secondary_label="Alt/Az transformed to RA/Dec", primary_color="tab:orange", secondary_color="tab:blue", x_key="ra", y_key="dec", panel_name="RA/Dec", ) - - def _plot_radec_track_group(self, axis, group, color, label=None, linestyle="-"): - """Plot a single inferred RA/Dec track group on the Alt/Az panel.""" - sky_coordinates = SkyCoord( - ra=np.array([point["ra"] for point in group]) * u.deg, - dec=np.array([point["dec"] for point in group]) * u.deg, - frame="icrs", - ) - altaz = sky_coordinates.transform_to(AltAz(obstime=self.time, location=self.location)) - visible_segments = self._split_visible_segments(altaz.alt.deg > 0.0) - - plotted_segments = 0 - show_label = label - for segment in visible_segments: - axis.plot( - altaz.az.rad[segment], - 90.0 - altaz.alt.deg[segment], - color=color, - lw=DEFAULT_GRID_LINE_WIDTH, - linestyle=linestyle, - alpha=0.8, - label=show_label, - zorder=3, - ) - plotted_segments += 1 - show_label = None - return plotted_segments - - def _plot_inferred_radec_grid(self, axis, plot_points): - """Plot thin inferred RA/Dec grid lines on the Alt/Az panel.""" - track_groups = self.infer_radec_grid_tracks(plot_points) - - plotted_tracks = 0 - for index, group in enumerate(track_groups["declination_tracks"]): - plotted_tracks += self._plot_radec_track_group( - axis, - group, - color="0.45", - linestyle="-", - label="Inferred Dec grid" if index == 0 else None, - ) - - for index, group in enumerate(track_groups["right_ascension_tracks"]): - plotted_tracks += self._plot_radec_track_group( - axis, - group, - color="0.6", - linestyle=":", - label="Inferred RA grid" if index == 0 else None, - ) - - if plotted_tracks == 0: - logger.info("No inferred RA/Dec grid tracks available for plotting") - else: - logger.info(f"Plotted {plotted_tracks} inferred RA/Dec grid track segments") - - return plotted_tracks - - def _plot_ra_dec_tracks(self, axis, dec_values): - """ - Plot RA/Dec coordinate tracks on the polar projection. - - Parameters - ---------- - axis : matplotlib.axes.Axes - Polar projection axes. - dec_values : list of float - List of declination values in degrees. - """ - colors = plt.get_cmap("viridis")(np.linspace(0, 1, len(dec_values))) - - for dec_val, color in zip(dec_values, colors): - right_ascension = np.linspace(0, 360, 400) * u.deg - declination = dec_val * u.deg - - coords = SkyCoord(ra=right_ascension, dec=declination, frame="icrs") - altaz = coords.transform_to(AltAz(obstime=self.time, location=self.location)) - - az_track = altaz.az.rad - zenith_track = 90 - altaz.alt.deg - - mask = altaz.alt.deg > 0 - if np.any(mask): - axis.plot( - az_track[mask], - zenith_track[mask], - color=color, - lw=2, - label=f"Dec = {dec_val:.1f} deg", - alpha=0.7, - ) - - logger.info(f"Plotted RA/Dec tracks for {len(dec_values)} declination values") diff --git a/tests/unit_tests/applications/test_production_generate_grid.py b/tests/unit_tests/applications/test_production_generate_grid.py index ebc85cc529..4cf70cda50 100644 --- a/tests/unit_tests/applications/test_production_generate_grid.py +++ b/tests/unit_tests/applications/test_production_generate_grid.py @@ -75,3 +75,12 @@ def test_add_arguments_accepts_zenith_angle_scaling_factor(): args = parser.parse_args(["--zenith_angle_scaling_factor", "2.5"]) assert args.zenith_angle_scaling_factor == pytest.approx(2.5) + + +def test_add_arguments_accepts_direction_grid_density(): + parser = CommandLineParser() + app._add_arguments(parser) + + args = parser.parse_args(["--direction_grid_density", "1.5"]) + + assert args.direction_grid_density == pytest.approx(1.5) diff --git a/tests/unit_tests/production_configuration/test_job_grid_io.py b/tests/unit_tests/production_configuration/test_job_grid_io.py index 0e739f7249..e667fce4c6 100644 --- a/tests/unit_tests/production_configuration/test_job_grid_io.py +++ b/tests/unit_tests/production_configuration/test_job_grid_io.py @@ -12,6 +12,8 @@ def _job_rows(): "primary": "gamma", "azimuth_angle": 45 * u.deg, "zenith_angle": 20 * u.deg, + "ra": 123 * u.deg, + "dec": -45 * u.deg, "energy_min": 30 * u.GeV, "energy_max": 10 * u.TeV, "core_scatter_number": 10, @@ -46,6 +48,8 @@ def test_serialize_and_read_job_grid_ecsv(tmp_test_directory): assert rows[0]["energy_min"] == 30 * u.GeV assert rows[0]["core_scatter_number"] == 10 assert rows[0]["array_layout_name"] == "CTAO-North-Alpha" + assert rows[0]["ra"] == 123 * u.deg + assert rows[0]["dec"] == -45 * u.deg def test_serialize_job_grid_rejects_non_ecsv_output(tmp_test_directory): diff --git a/tests/unit_tests/production_configuration/test_observation_grid.py b/tests/unit_tests/production_configuration/test_observation_grid.py index 97d60a56b1..c2b0dd191d 100644 --- a/tests/unit_tests/production_configuration/test_observation_grid.py +++ b/tests/unit_tests/production_configuration/test_observation_grid.py @@ -146,12 +146,12 @@ def test_generate_extra_axis_combinations_returns_mesh_for_remaining_axes(): assert combinations.shape == (4, 2) -def test_create_circular_binning_uses_shorter_counterclockwise_path(): +def test_create_circular_binning_uses_directed_span_from_start_to_end(): engine = ProductionGridEngine(axes={"axes": {}}) binning = engine.create_circular_binning((10, 350), 3) - assert np.allclose(binning, [10, 0, 350]) + assert np.allclose(binning, [10, 180, 350]) def test_create_circular_binning_uses_clockwise_path_when_shorter(): @@ -162,6 +162,55 @@ def test_create_circular_binning_uses_clockwise_path_when_shorter(): assert np.allclose(binning, [350, 0, 10]) +def test_create_circular_binning_covers_requested_range_0_to_240(): + engine = ProductionGridEngine(axes={"axes": {}}) + + binning = engine.create_circular_binning((0, 240), 3) + + assert np.allclose(binning, [0, 120, 240]) + + +def test_generate_horizontal_grid_uses_adaptive_azimuth_bins_per_zenith_row(): + engine = ProductionGridEngine( + axes={ + "axes": { + "zenith_angle": {"range": [0, 60], "binning": 61, "units": "deg"}, + "azimuth": { + "range": [0, 240], + "binning": 10, + "units": "deg", + "direction_grid_density": 0.05, + }, + "nsb_level": {"range": [1, 1], "binning": 1, "units": "1"}, + } + }, + coordinate_system="horizontal", + ) + + grid = engine._generate_horizontal_grid() + + assert len(grid) > 0 + + azimuth_counts_by_zenith = {} + for point in grid: + zenith_value = round(point["zenith_angle"].to_value(u.deg), 6) + azimuth_counts_by_zenith.setdefault(zenith_value, set()).add( + round(point["azimuth"].to_value(u.deg), 6) + ) + + assert len(azimuth_counts_by_zenith[0.0]) == 1 + assert len(azimuth_counts_by_zenith[30.0]) > len(azimuth_counts_by_zenith[0.0]) + assert len(azimuth_counts_by_zenith[60.0]) > len(azimuth_counts_by_zenith[30.0]) + + +def test_create_circular_binning_treats_full_circle_range_as_full_span(): + engine = ProductionGridEngine(axes={"axes": {}}) + + binning = engine.create_circular_binning((0, 360), 4) + + assert np.allclose(binning, [0, 90, 180, 270]) + + def test_convert_altaz_to_radec_raises_without_time_of_observation(): engine = ProductionGridEngine(axes={"axes": {}}, time_of_observation=None) @@ -322,8 +371,8 @@ def test_generate_horizontal_grid_with_circular_azimuth_binning_uses_correct_ind } } ) - # Counterclockwise shortest path gives non-monotonic order: [10, 0, 350] - assert np.allclose(engine.target_values["azimuth"].value, np.array([10.0, 0.0, 350.0])) + # Directed span from start to end gives [10, 180, 350] + assert np.allclose(engine.target_values["azimuth"].value, np.array([10.0, 180.0, 350.0])) engine.interpolated_limits = { "lower_energy_threshold": np.array([[[0.1], [0.2], [0.3]]]), @@ -336,7 +385,7 @@ def test_generate_horizontal_grid_with_circular_azimuth_binning_uses_correct_ind } assert by_azimuth[10.0] == pytest.approx(0.1) - assert by_azimuth[0.0] == pytest.approx(0.2) + assert by_azimuth[180.0] == pytest.approx(0.2) assert by_azimuth[350.0] == pytest.approx(0.3) diff --git a/tests/unit_tests/production_configuration/test_simulation_jobs.py b/tests/unit_tests/production_configuration/test_simulation_jobs.py index 3ef6d3ebaa..570ce2ee3a 100644 --- a/tests/unit_tests/production_configuration/test_simulation_jobs.py +++ b/tests/unit_tests/production_configuration/test_simulation_jobs.py @@ -91,6 +91,139 @@ def test_build_axes_dict_from_cli_args_accepts_compact_axis_definitions(): } +def test_build_axes_dict_from_cli_args_derives_horizontal_binning_from_density(): + axes = build_axes_dict_from_cli_args( + { + "direction_grid_density": 1.0, + "axis": [ + ["azimuth", "310", "deg", "20", "deg", "3", "linear"], + ["zenith", "30", "deg", "40", "deg", "2"], + ["nsb", "4", "MHz", "5", "MHz", "2"], + ["offset", "0", "deg", "10", "deg", "2"], + ], + } + ) + + assert axes["azimuth"]["binning"] == 41 + assert axes["zenith_angle"]["binning"] == 10 + + +def test_build_axes_dict_from_cli_args_derives_radec_binning_from_density(): + axes = build_axes_dict_from_cli_args( + { + "direction_grid_density": 1.0, + "axis": [ + ["ra", "0", "deg", "360", "deg", "36", "linear"], + ["dec", "-90", "deg", "90", "deg", "18", "linear"], + ["nsb", "4", "MHz", "4", "MHz", "1", "linear"], + ["offset", "0", "deg", "10", "deg", "2", "linear"], + ], + } + ) + + assert axes["ra"]["binning"] == 230 + assert axes["dec"]["binning"] == 180 + + +def test_build_axes_dict_from_cli_args_reduces_ra_binning_towards_dec_poles(): + axes = build_axes_dict_from_cli_args( + { + "direction_grid_density": 1.0, + "axis": [ + ["ra", "0", "deg", "360", "deg", "36", "linear"], + ["dec", "80", "deg", "90", "deg", "10", "linear"], + ["nsb", "4", "MHz", "4", "MHz", "1", "linear"], + ["offset", "0", "deg", "10", "deg", "2", "linear"], + ], + } + ) + + assert axes["ra"]["binning"] == 32 + assert axes["dec"]["binning"] == 10 + + +def test_build_axes_dict_from_cli_args_reduces_az_binning_towards_zenith_pole(): + axes = build_axes_dict_from_cli_args( + { + "direction_grid_density": 1.0, + "axis": [ + ["azimuth", "0", "deg", "180", "deg", "36", "linear"], + ["zenith", "0", "deg", "10", "deg", "10", "linear"], + ["nsb", "4", "MHz", "4", "MHz", "1", "linear"], + ["offset", "0", "deg", "10", "deg", "2", "linear"], + ], + } + ) + + assert axes["azimuth"]["binning"] == 16 + assert axes["zenith_angle"]["binning"] == 10 + + +def test_build_axes_dict_from_cli_args_uses_directed_azimuth_span_for_density(): + axes = build_axes_dict_from_cli_args( + { + "direction_grid_density": 1.0, + "axis": [ + ["azimuth", "0", "deg", "240", "deg", "36", "linear"], + ["zenith", "0", "deg", "70", "deg", "2", "linear"], + ["nsb", "4", "MHz", "4", "MHz", "1", "linear"], + ["offset", "0", "deg", "10", "deg", "2", "linear"], + ], + } + ) + + assert axes["azimuth"]["binning"] == 130 + assert axes["zenith_angle"]["binning"] == 70 + + +def test_build_axes_dict_from_cli_args_sets_horizontal_density_metadata_for_adaptive_grid(): + axes = build_axes_dict_from_cli_args( + { + "direction_grid_density": 1.0, + "axis": [ + ["azimuth", "0", "deg", "240", "deg", "36", "linear"], + ["zenith", "0", "deg", "70", "deg", "2", "linear"], + ["nsb", "4", "MHz", "4", "MHz", "1", "linear"], + ["offset", "0", "deg", "10", "deg", "2", "linear"], + ], + } + ) + + assert axes["azimuth"]["direction_grid_density"] == pytest.approx(1.0) + + +def test_build_axes_dict_from_cli_args_uses_full_circle_span_for_azimuth_density(): + axes = build_axes_dict_from_cli_args( + { + "direction_grid_density": 1.0, + "axis": [ + ["azimuth", "0", "deg", "360", "deg", "36", "linear"], + ["zenith", "0", "deg", "70", "deg", "2", "linear"], + ["nsb", "4", "MHz", "4", "MHz", "1", "linear"], + ["offset", "0", "deg", "10", "deg", "2", "linear"], + ], + } + ) + + assert axes["azimuth"]["binning"] == 194 + assert axes["zenith_angle"]["binning"] == 70 + + +def test_build_axes_dict_from_cli_args_rejects_non_positive_density(): + with pytest.raises(ValueError, match="direction_grid_density must be strictly positive"): + build_axes_dict_from_cli_args( + { + "direction_grid_density": 0, + "axis": [ + ["azimuth", "310", "deg", "20", "deg", "3"], + ["zenith", "30", "deg", "40", "deg", "2"], + ["nsb", "4", "MHz", "5", "MHz", "2"], + ["offset", "0", "deg", "10", "deg", "2"], + ], + } + ) + + def test_build_axes_dict_from_cli_args_accepts_config_wrapped_axis_strings(): axes = build_axes_dict_from_cli_args( { @@ -172,6 +305,22 @@ def test_build_job_grid_metadata_includes_job_context(): assert metadata["corsika_limits"] == "limits.ecsv" +def test_build_job_grid_metadata_raises_for_radec_without_site(): + with pytest.raises(ValueError, match="site is required"): + build_job_grid_metadata( + { + "site": None, + "simulation_software": "corsika_sim_telarray", + "axis": [ + ["ra", "0", "deg", "1", "deg", "2"], + ["dec", "0", "deg", "1", "deg", "2"], + ], + "time_of_observation": "2017-09-16 00:00:00", + "corsika_limits": "limits.ecsv", + } + ) + + @patch("simtools.production_configuration.simulation_jobs.SiteModel") def test_build_observing_location_uses_site_model(mock_site_model): model = mock_site_model.return_value @@ -716,6 +865,8 @@ def test_build_simulation_jobs_expands_runs_from_observation_grid( { "azimuth": 180 * u.deg, "zenith_angle": 20 * u.deg, + "ra": 123 * u.deg, + "dec": -45 * u.deg, "lower_energy_threshold": 40 * u.GeV, "scatter_radius": 100 * u.m, "viewcone_radius": 2 * u.deg, @@ -747,6 +898,8 @@ def test_build_simulation_jobs_expands_runs_from_observation_grid( assert rows[0]["core_scatter_max"] == 100 * u.m assert rows[0]["view_cone_max"] == 2 * u.deg assert rows[0]["showers_per_run"] == 5 + assert rows[0]["ra"] == 123 * u.deg + assert rows[0]["dec"] == -45 * u.deg @patch("simtools.production_configuration.simulation_jobs._generate_observation_grids_per_layout") diff --git a/tests/unit_tests/visualization/test_plot_production_grid.py b/tests/unit_tests/visualization/test_plot_production_grid.py index 2591753870..3666477c48 100644 --- a/tests/unit_tests/visualization/test_plot_production_grid.py +++ b/tests/unit_tests/visualization/test_plot_production_grid.py @@ -2,37 +2,15 @@ from pathlib import Path -import astropy.units as u import matplotlib.pyplot as plt import pytest -from astropy.coordinates import AltAz, EarthLocation, SkyCoord from astropy.table import Table -from astropy.time import Time -from astropy.utils import iers from simtools.visualization.plot_production_grid import ( DEFAULT_OUTPUT_FILE_STEM, ProductionGridPlotter, ) -pytestmark = pytest.mark.filterwarnings("ignore::astropy.utils.iers.IERSWarning") - - -@pytest.fixture(autouse=True, scope="module") -def disable_iers_auto_download(): - """Disable IERS auto-download during tests to avoid network dependency.""" - previous_auto_download = iers.conf.auto_download - iers.conf.auto_download = False - try: - yield - finally: - iers.conf.auto_download = previous_auto_download - - -SITE_LOCATION_LAT = 28.76 -SITE_LOCATION_LON = -17.89 -SITE_LOCATION_HEIGHT = 2200.0 - def _write_grid_file(tmp_test_directory, file_name, grid_points): """Write grid points to a temporary ECSV file.""" @@ -60,53 +38,28 @@ def _write_grid_file(tmp_test_directory, file_name, grid_points): table["nsb_level"].unit = "MHz" if "offset" in table.colnames: table["offset"].unit = "deg" - table.meta["time_of_observation_utc"] = "2025-01-01T00:00:00.000" table.write(file_path, format="ascii.ecsv", overwrite=True) return file_path -def _build_radec_mesh_grid_points(location, observation_time): - """Build a small RA/Dec mesh around local sidereal time.""" - lst = observation_time.sidereal_time("apparent", longitude=location.lon).deg - ra_values = [(lst - 5.0) % 360.0, (lst + 5.0) % 360.0] - dec_values = [20.0, 30.0] - - return [ - {"ra": {"value": ra_value, "unit": "deg"}, "dec": {"value": dec_value, "unit": "deg"}} - for dec_value in dec_values - for ra_value in ra_values - ] - - -def _create_plotter(grid_file, observation_time, output_path): - """Create a plotter with the standard test site parameters.""" +def _create_plotter(grid_file, output_path): + """Create a plotter for file-driven plotting.""" return ProductionGridPlotter( grid_points_file=grid_file, - site_location_lat=SITE_LOCATION_LAT, - site_location_lon=SITE_LOCATION_LON, - site_location_height=SITE_LOCATION_HEIGHT, - observation_time=observation_time, output_path=output_path, ) -def test_normalize_altaz_point_creates_radec_coordinates(tmp_test_directory): - """Test that native Alt/Az points are converted to RA/Dec for the equatorial panel.""" +def test_normalize_altaz_point(tmp_test_directory): + """Keep native Alt/Az values and no RA/Dec when absent in file.""" grid_file = _write_grid_file( tmp_test_directory, "grid_altaz.ecsv", - [ - { - "azimuth": {"value": 180.0, "unit": "deg"}, - "zenith_angle": {"value": 20.0, "unit": "deg"}, - "nsb_level": {"value": 0.0, "unit": "MHz"}, - } - ], + [{"azimuth": 180.0, "zenith_angle": 20.0, "nsb_level": 0.0}], ) plotter = _create_plotter( grid_file=grid_file, - observation_time="2025-01-01 00:00:00", output_path=Path(tmp_test_directory) / "output", ) @@ -117,151 +70,128 @@ def test_normalize_altaz_point_creates_radec_coordinates(tmp_test_directory): assert normalized_points[0]["visible_in_altaz"] is True assert normalized_points[0]["azimuth"] == pytest.approx(180.0) assert normalized_points[0]["zenith"] == pytest.approx(20.0) - assert normalized_points[0]["ra"] is not None - assert normalized_points[0]["dec"] is not None - + assert normalized_points[0]["ra"] is None + assert normalized_points[0]["dec"] is None -def test_normalize_radec_point_projects_to_altaz(tmp_test_directory): - """Test that native RA/Dec points are converted to Alt/Az for the local panel.""" - location = EarthLocation( - lat=SITE_LOCATION_LAT * u.deg, - lon=SITE_LOCATION_LON * u.deg, - height=SITE_LOCATION_HEIGHT * u.m, - ) - observation_time = Time("2025-01-01 00:00:00") - source_altaz = SkyCoord( - AltAz( - alt=65.0 * u.deg, - az=210.0 * u.deg, - obstime=observation_time, - location=location, - ) - ) - source_radec = source_altaz.icrs +def test_normalize_flattened_job_grid_altaz_columns(tmp_test_directory): + """Handle job-grid style flattened coordinate columns.""" grid_file = _write_grid_file( tmp_test_directory, - "grid_radec.ecsv", + "grid_altaz_flattened.ecsv", [ { - "ra": {"value": source_radec.ra.deg, "unit": "deg"}, - "dec": {"value": source_radec.dec.deg, "unit": "deg"}, + "azimuth_angle_value": 0.0, + "azimuth_angle_unit": "deg", + "zenith_angle_value": 70.0, + "zenith_angle_unit": "deg", + "primary": "gamma", } ], ) plotter = _create_plotter( grid_file=grid_file, - observation_time=str(observation_time.value), output_path=Path(tmp_test_directory) / "output", ) normalized_points = plotter.normalize_grid_points() assert len(normalized_points) == 1 - assert normalized_points[0]["native_frame"] == "radec" - assert normalized_points[0]["visible_in_altaz"] is True - assert normalized_points[0]["azimuth"] == pytest.approx(210.0, abs=0.05) - assert normalized_points[0]["zenith"] == pytest.approx(25.0, abs=0.05) - assert normalized_points[0]["ra"] == pytest.approx(source_radec.ra.deg, abs=0.05) - assert normalized_points[0]["dec"] == pytest.approx(source_radec.dec.deg, abs=0.05) - - -def test_infer_radec_grid_tracks_from_native_points(tmp_test_directory): - """Infer RA and Dec grid tracks from a native RA/Dec mesh.""" - location = EarthLocation( - lat=SITE_LOCATION_LAT * u.deg, - lon=SITE_LOCATION_LON * u.deg, - height=SITE_LOCATION_HEIGHT * u.m, + assert normalized_points[0]["native_frame"] == "altaz" + assert normalized_points[0]["azimuth"] == pytest.approx(0.0) + assert normalized_points[0]["zenith"] == pytest.approx(70.0) + + +def test_normalize_altaz_keeps_explicit_radec_columns(tmp_test_directory): + """Use explicit RA/Dec columns from the same grid row when available.""" + grid_file = _write_grid_file( + tmp_test_directory, + "grid_altaz_with_radec.ecsv", + [ + { + "azimuth_angle_value": 10.0, + "azimuth_angle_unit": "deg", + "zenith_angle_value": 20.0, + "zenith_angle_unit": "deg", + "ra": 120.0, + "dec": -20.0, + } + ], ) - observation_time = Time("2025-01-01 00:00:00") - grid_points = _build_radec_mesh_grid_points(location, observation_time) - grid_file = _write_grid_file(tmp_test_directory, "grid_radec_mesh.ecsv", grid_points) plotter = _create_plotter( grid_file=grid_file, - observation_time=str(observation_time.value), output_path=Path(tmp_test_directory) / "output", ) normalized_points = plotter.normalize_grid_points() - track_groups = plotter.infer_radec_grid_tracks(normalized_points) - assert len(track_groups["declination_tracks"]) == 2 - assert len(track_groups["right_ascension_tracks"]) == 2 - assert all(len(group) == 2 for group in track_groups["declination_tracks"]) - assert all(len(group) == 2 for group in track_groups["right_ascension_tracks"]) + assert len(normalized_points) == 1 + assert normalized_points[0]["native_frame"] == "altaz" + assert normalized_points[0]["ra"] == pytest.approx(120.0) + assert normalized_points[0]["dec"] == pytest.approx(-20.0) -def test_plot_sky_projection_creates_outputs(tmp_test_directory): - """Test that the sky projection plot is written to disk.""" +def test_normalize_radec_point_without_altaz_projection(tmp_test_directory): + """Keep native RA/Dec points without deriving Alt/Az.""" grid_file = _write_grid_file( tmp_test_directory, - "grid_wrapped.ecsv", - [ - { - "azimuth": {"value": 310.0, "unit": "deg"}, - "zenith_angle": {"value": 30.0, "unit": "deg"}, - "nsb_level": {"value": 4.0, "unit": "MHz"}, - "offset": {"value": 0.0, "unit": "deg"}, - } - ], + "grid_radec.ecsv", + [{"ra": 155.0, "dec": 30.0}], ) - output_path = Path(tmp_test_directory) / "output" plotter = _create_plotter( grid_file=grid_file, - observation_time="2025-01-01 00:00:00", - output_path=output_path, + output_path=Path(tmp_test_directory) / "output", ) - plotter.plot_sky_projection(plot_ra_dec_tracks=True, dec_values=[20.0, 30.0]) + normalized_points = plotter.normalize_grid_points() - assert (output_path / f"{DEFAULT_OUTPUT_FILE_STEM}.png").exists() + assert len(normalized_points) == 1 + assert normalized_points[0]["native_frame"] == "radec" + assert normalized_points[0]["azimuth"] is None + assert normalized_points[0]["zenith"] is None + assert normalized_points[0]["ra"] == pytest.approx(155.0) + assert normalized_points[0]["dec"] == pytest.approx(30.0) -def test_plot_sky_projection_infers_radec_grid_tracks(tmp_test_directory): - """Plot inferred RA/Dec grid tracks.""" - location = EarthLocation( - lat=SITE_LOCATION_LAT * u.deg, - lon=SITE_LOCATION_LON * u.deg, - height=SITE_LOCATION_HEIGHT * u.m, +def test_plot_sky_projection_creates_output_altaz_only(tmp_test_directory): + """Write the sky projection plot with Alt/Az data only.""" + grid_file = _write_grid_file( + tmp_test_directory, + "grid_altaz_only.ecsv", + [{"azimuth": 120.0, "zenith_angle": 35.0}], ) - observation_time = Time("2025-01-01 00:00:00") - grid_points = _build_radec_mesh_grid_points(location, observation_time) - - grid_file = _write_grid_file(tmp_test_directory, "grid_radec_tracks.ecsv", grid_points) output_path = Path(tmp_test_directory) / "output" + plotter = _create_plotter( grid_file=grid_file, - observation_time=str(observation_time.value), output_path=output_path, ) - plotter.plot_sky_projection(plot_ra_dec_tracks=True) + plotter.plot_sky_projection(plot_ra_dec_tracks=True, dec_values=[20.0]) assert (output_path / f"{DEFAULT_OUTPUT_FILE_STEM}.png").exists() -def test_load_grid_points_from_ecsv(tmp_test_directory): - """Test that ECSV grid-point files are loaded and normalized.""" - grid_file = Path(tmp_test_directory) / "grid_points.ecsv" - table = Table(rows=[{"azimuth": 180.0, "zenith_angle": 20.0, "nsb_level": 1.0}]) - table["azimuth"].unit = "deg" - table["zenith_angle"].unit = "deg" - table["nsb_level"].unit = "MHz" - table.meta["time_of_observation_utc"] = "2025-01-01T00:00:00.000" - table.write(grid_file, format="ascii.ecsv", overwrite=True) +def test_plot_sky_projection_creates_output_with_radec_panel(tmp_test_directory): + """Write the sky projection plot with both Alt/Az and RA/Dec data.""" + grid_file = _write_grid_file( + tmp_test_directory, + "grid_altaz_radec.ecsv", + [{"azimuth": 100.0, "zenith_angle": 25.0, "ra": 180.0, "dec": -10.0}], + ) + output_path = Path(tmp_test_directory) / "output" plotter = _create_plotter( grid_file=grid_file, - observation_time=None, - output_path=Path(tmp_test_directory) / "output", + output_path=output_path, ) - normalized_points = plotter.normalize_grid_points() - assert len(normalized_points) == 1 - assert normalized_points[0]["native_frame"] == "altaz" + plotter.plot_sky_projection() + + assert (output_path / f"{DEFAULT_OUTPUT_FILE_STEM}.png").exists() def test_load_grid_points_file_not_found(tmp_test_directory): @@ -271,7 +201,6 @@ def test_load_grid_points_file_not_found(tmp_test_directory): with pytest.raises(FileNotFoundError, match="Grid points file not found"): _create_plotter( grid_file=missing_file, - observation_time="2025-01-01 00:00:00", output_path=Path(tmp_test_directory) / "output", ) @@ -284,7 +213,6 @@ def test_load_grid_points_wrong_suffix(tmp_test_directory): with pytest.raises(ValueError, match="must be ECSV"): _create_plotter( grid_file=wrong_file, - observation_time="2025-01-01 00:00:00", output_path=Path(tmp_test_directory) / "output", ) @@ -310,11 +238,10 @@ def test_configure_radec_axis_expands_flat_ranges(tmp_test_directory): grid_file = _write_grid_file( tmp_test_directory, "grid_single_radec.ecsv", - [{"ra": {"value": 40.0, "unit": "deg"}, "dec": {"value": 20.0, "unit": "deg"}}], + [{"ra": 40.0, "dec": 20.0}], ) plotter = _create_plotter( grid_file=grid_file, - observation_time="2025-01-01 00:00:00", output_path=Path(tmp_test_directory) / "output", ) @@ -334,7 +261,6 @@ def test_plot_frame_points_logs_no_valid_points(tmp_test_directory, caplog): grid_file = _write_grid_file(tmp_test_directory, "grid_empty.ecsv", []) plotter = _create_plotter( grid_file=grid_file, - observation_time="2025-01-01 00:00:00", output_path=Path(tmp_test_directory) / "output", ) @@ -346,8 +272,6 @@ def test_plot_frame_points_logs_no_valid_points(tmp_test_directory, caplog): plot_points=[], primary_frame="altaz", secondary_frame="radec", - primary_label="A", - secondary_label="B", primary_color="tab:blue", secondary_color="tab:orange", x_key="azimuth", @@ -361,11 +285,10 @@ def test_plot_frame_points_logs_no_valid_points(tmp_test_directory, caplog): def test_plot_altaz_points_logs_hidden_radec_points(tmp_test_directory, caplog): - """Log info when RA/Dec points are below horizon and skipped in Alt/Az panel.""" + """Log info when RA/Dec points are not visible in Alt/Az panel.""" grid_file = _write_grid_file(tmp_test_directory, "grid_empty_altaz.ecsv", []) plotter = _create_plotter( grid_file=grid_file, - observation_time="2025-01-01 00:00:00", output_path=Path(tmp_test_directory) / "output", ) @@ -389,71 +312,19 @@ def test_plot_altaz_points_logs_hidden_radec_points(tmp_test_directory, caplog): plt.close(figure) -def test_plot_inferred_radec_grid_logs_no_tracks(tmp_test_directory, caplog): - """Log info when no inferred RA/Dec tracks can be plotted.""" - grid_file = _write_grid_file(tmp_test_directory, "grid_empty_tracks.ecsv", []) - plotter = _create_plotter( - grid_file=grid_file, - observation_time="2025-01-01 00:00:00", - output_path=Path(tmp_test_directory) / "output", - ) - - figure = plt.figure() - axis = figure.add_subplot(1, 1, 1, projection="polar") - try: - with caplog.at_level("INFO"): - plotted = plotter._plot_inferred_radec_grid(axis, plot_points=[]) - assert plotted == 0 - assert "No inferred RA/Dec grid tracks available for plotting" in caplog.text - finally: - plt.close(figure) - - -def test_iers_disabled_with_env_plotter(monkeypatch, tmp_test_directory): - from simtools.application_control import _configure_iers_from_env - - iers.conf.auto_download = True - iers.conf.auto_max_age = 30 - - monkeypatch.setenv("SIMTOOLS_OFFLINE_IERS", "1") - - _configure_iers_from_env() - +def test_plot_sky_projection_logs_tracks_disabled(tmp_test_directory, caplog): + """Log that RA/Dec tracks are disabled in file-driven mode.""" grid_file = _write_grid_file( tmp_test_directory, - "grid.ecsv", - [{"azimuth": 0.0, "zenith_angle": 0.0}], + "grid_altaz_with_radec_for_tracks.ecsv", + [{"azimuth": 10.0, "zenith_angle": 20.0, "ra": 120.0, "dec": -20.0}], ) - - _create_plotter( + plotter = _create_plotter( grid_file=grid_file, - observation_time="2025-01-01 00:00:00", output_path=Path(tmp_test_directory) / "output", ) - assert iers.conf.auto_download is False - - -def test_iers_not_modified_without_env_plotter(monkeypatch, tmp_test_directory): - from simtools.application_control import _configure_iers_from_env - - iers.conf.auto_download = True - iers.conf.auto_max_age = 30 - - monkeypatch.delenv("SIMTOOLS_OFFLINE_IERS", raising=False) - - _configure_iers_from_env() - - grid_file = _write_grid_file( - tmp_test_directory, - "grid.ecsv", - [{"azimuth": 0.0, "zenith_angle": 0.0}], - ) - - _create_plotter( - grid_file=grid_file, - observation_time="2025-01-01 00:00:00", - output_path=Path(tmp_test_directory) / "output", - ) + with caplog.at_level("INFO"): + plotter.plot_sky_projection(plot_ra_dec_tracks=True) - assert iers.conf.auto_download is True + assert "RA/Dec tracks are disabled in file-driven plotting mode" in caplog.text From f8c13e5888a986b6d30f1ad77bdb07929ce2a4e4 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Thu, 28 May 2026 11:25:07 +0200 Subject: [PATCH 02/10] grid density --- .../applications/production_generate_grid.py | 49 +++++- .../observation_grid.py | 165 ++++++++++++++++++ .../simulation_jobs.py | 70 +++++++- .../visualization/plot_production_grid.py | 5 +- ...ction_generate_grid_horizontal_density.yml | 35 ++++ ...roduction_generate_grid_ra_dec_density.yml | 42 +++++ .../test_observation_grid.py | 117 +++++++++++++ .../test_simulation_jobs.py | 75 +++++++- .../test_plot_production_grid.py | 30 ++++ 9 files changed, 576 insertions(+), 12 deletions(-) create mode 100644 tests/integration_tests/config/production_generate_grid_horizontal_density.yml create mode 100644 tests/integration_tests/config/production_generate_grid_ra_dec_density.yml diff --git a/src/simtools/applications/production_generate_grid.py b/src/simtools/applications/production_generate_grid.py index 8bc6b50d74..a6074e31ae 100644 --- a/src/simtools/applications/production_generate_grid.py +++ b/src/simtools/applications/production_generate_grid.py @@ -26,6 +26,14 @@ Direction-grid density in points per deg^2. If provided, direction-axis binning (azimuth/zenith or ra/dec) is derived from range and density, while min/max values are kept from ``--axis`` definitions. + In ``ra_dec`` mode, local-sky constraints can be defined with + ``local_zenith_range`` and ``local_azimuth_range``. +local_zenith_range (quantity pair, optional) + Zenith range (deg) used to filter generated RA/Dec density nodes in local sky, + for example ``0 deg 70 deg``. +local_azimuth_range (quantity pair, optional) + Directed azimuth range (deg) used to filter generated RA/Dec density nodes in + local sky, for example ``300 deg 60 deg``. time_of_observation (str, optional) Time of the observation in UTC (format: 'YYYY-MM-DD HH:MM:SS'). Used only if RA/Dec axes are provided (for coordinate transforms and sidereal-time @@ -69,6 +77,23 @@ --axis offset 0 deg 10 deg 2 linear \ --time_of_observation "2017-09-16 00:00:00" \ --corsika_limits tests/resources/corsika_simulation_limits/merged_corsika_limits.ecsv + +To generate an RA/Dec density grid constrained to local sky ranges (for example +full zenith coverage from 0 to 70 deg and a directed azimuth window), execute: + +.. code-block:: console + + simtools-production-generate-grid --site North --model_version 6.0.2 \ + --array_layout_name alpha \ + --axis ra 0 deg 360 deg 1 linear \ + --axis dec -40 deg 80 deg 1 linear \ + --axis nsb 4 MHz 4 MHz 1 linear \ + --axis offset 0 deg 10 deg 2 linear \ + --direction_grid_density 0.25 \ + --local_zenith_range 0 deg 70 deg \ + --local_azimuth_range 300 deg 60 deg \ + --time_of_observation "2017-09-16 00:00:00" \ + --corsika_limits tests/resources/corsika_simulation_limits/merged_corsika_limits.ecsv """ from simtools.application_control import build_application @@ -108,7 +133,29 @@ def _add_arguments(parser): default=None, help=( "Direction-grid density in points per deg^2. If set, direction-axis binning is " - "derived from axis ranges and this density." + "derived from axis ranges and this density. In ra_dec mode, optional zenith and " + "azimuth axes are interpreted as local-sky constraints for filtering generated " + "points (for example zenith 0..70 deg)." + ), + ) + parser.add_argument( + "--local_zenith_range", + nargs="+", + required=False, + default=None, + help=( + "Local zenith range (quantity pair) used to filter RA/Dec density points, " + "for example: --local_zenith_range 0 deg 70 deg" + ), + ) + parser.add_argument( + "--local_azimuth_range", + nargs="+", + required=False, + default=None, + help=( + "Local azimuth range (quantity pair) used to filter RA/Dec density points, " + "for example: --local_azimuth_range 300 deg 60 deg" ), ) parser.add_argument( diff --git a/src/simtools/production_configuration/observation_grid.py b/src/simtools/production_configuration/observation_grid.py index dadd25166c..09d4647187 100644 --- a/src/simtools/production_configuration/observation_grid.py +++ b/src/simtools/production_configuration/observation_grid.py @@ -214,6 +214,9 @@ def _add_lookup_limits_to_point(self, point, zenith, azimuth): def _generate_grid_from_radec_axes(self, include_horizontal_coordinates=False): """Generate grid points from explicit RA/Dec axes definitions.""" + if self._is_adaptive_radec_density_enabled(): + return self._generate_adaptive_radec_grid(include_horizontal_coordinates) + time_of_observation = self._require_time_of_observation() axis_keys = [key for key in self.target_values if key not in ("zenith_angle", "azimuth")] @@ -334,6 +337,168 @@ def _directed_circular_span_degrees(azimuth_range): return 360.0 return float((azimuth_range[1] - azimuth_range[0]) % 360.0) + @staticmethod + def _is_in_directed_azimuth_range(azimuth_values_deg, azimuth_range): + """Return mask of azimuth values inside directed circular range [start -> end].""" + start, end = azimuth_range + raw_span = abs(float(end) - float(start)) + if raw_span > 0.0 and np.isclose(raw_span % 360.0, 0.0): + return np.ones_like(azimuth_values_deg, dtype=bool) + + start_norm = float(start) % 360.0 + span = float((end - start) % 360.0) + offsets = (np.asarray(azimuth_values_deg) - start_norm) % 360.0 + return offsets <= (span + 1e-12) + + def _is_adaptive_radec_density_enabled(self): + """Return whether RA/Dec grid generation should use per-declination adaptive RA bins.""" + return ( + self.coordinate_system == "ra_dec" + and "ra" in self.axes + and "dec" in self.axes + and self.axes["ra"].get("direction_grid_density") is not None + ) + + def _compute_visible_radec_strip(self, dec_deg, spacing, lst_deg, time_of_observation): + """Compute visible RA, zenith, and azimuth samples for one declination strip.""" + cos_dec = np.cos(np.deg2rad(dec_deg)) + if cos_dec <= 0.0: + return None + + ha_step = spacing / cos_dec + ha_values = np.arange(-180.0, 180.0, ha_step) + if len(ha_values) == 0: + ha_values = np.array([0.0]) + + ra_values = ((lst_deg - ha_values) % 360.0) * u.deg + skycoord = SkyCoord( + ra=ra_values.to(u.deg), + dec=np.full(len(ra_values), dec_deg) * u.deg, + frame="icrs", + ) + altaz = skycoord.transform_to( + AltAz(location=self.observing_location, obstime=time_of_observation) + ) + + visible_mask = altaz.alt.to_value(u.deg) >= 0.0 + + zenith_values = (90.0 * u.deg - altaz.alt).to(u.deg) + azimuth_values = altaz.az.to(u.deg) + + ra_axis = self.axes.get("ra", {}) + + zenith_range = ra_axis.get("local_zenith_range") + if zenith_range is None: + zenith_axis = self.axes.get("zenith_angle") + if zenith_axis and "range" in zenith_axis: + zenith_range = zenith_axis["range"] + + if zenith_range is not None: + zenith_min, zenith_max = sorted( + ( + float(zenith_range[0]), + float(zenith_range[1]), + ) + ) + zenith_values_deg = zenith_values.to_value(u.deg) + visible_mask &= (zenith_values_deg >= zenith_min) & (zenith_values_deg <= zenith_max) + + azimuth_range = ra_axis.get("local_azimuth_range") + if azimuth_range is None: + azimuth_axis = self.axes.get("azimuth") + if azimuth_axis and "range" in azimuth_axis: + azimuth_range = azimuth_axis["range"] + + if azimuth_range is not None: + azimuth_values_deg = azimuth_values.to_value(u.deg) + visible_mask &= self._is_in_directed_azimuth_range( + azimuth_values_deg, + azimuth_range, + ) + + if not np.any(visible_mask): + return None + + return ( + dec_deg * u.deg, + ra_values[visible_mask], + zenith_values[visible_mask], + azimuth_values[visible_mask], + ) + + def _generate_adaptive_radec_grid(self, include_horizontal_coordinates=False): + """Generate RA/Dec grid with Hour-Angle stepping adapted per declination strip. + + Pointings are arranged along declination lines evenly spaced in declination. + For each strip the telescope pointings are stepped through Hour Angle (equivalent + to stepping through time), tracing the actual arc a source at that declination + sweeps across the local sky. Hour-angle spacing is scaled by 1 / cos(dec) + to maintain uniform density per solid angle. + + Parameters + ---------- + include_horizontal_coordinates : bool, optional + If True, include zenith_angle and azimuth in each grid point. + + Returns + ------- + list of dict + Grid points with ra, dec and optionally zenith_angle, azimuth. + """ + time_of_observation = self._require_time_of_observation() + density = float(self.axes["ra"]["direction_grid_density"]) + spacing = 1.0 / np.sqrt(density) + dec_min, dec_max = sorted( + (float(self.axes["dec"]["range"][0]), float(self.axes["dec"]["range"][1])) + ) + dec_values_deg = np.arange(dec_min, dec_max + 0.5 * spacing, spacing) + + lst_deg = time_of_observation.sidereal_time( + "apparent", longitude=self.observing_location.lon + ).deg + + extra_keys, extra_units, extra_combinations = self._generate_extra_axis_combinations( + excluded_keys=("ra", "dec") + ) + + grid_points = [] + for dec_deg in dec_values_deg: + strip = self._compute_visible_radec_strip( + dec_deg=dec_deg, + spacing=spacing, + lst_deg=lst_deg, + time_of_observation=time_of_observation, + ) + if strip is None: + continue + dec, visible_ra_values, visible_zenith_values, visible_azimuth_values = strip + + for extra_combination in extra_combinations: + point_base = { + key: Quantity(extra_combination[i], extra_units[i]) + for i, key in enumerate(extra_keys) + } + for ra, zenith, azimuth in zip( + visible_ra_values, + visible_zenith_values, + visible_azimuth_values, + strict=True, + ): + point = {**point_base, "ra": ra, "dec": dec} + + if include_horizontal_coordinates: + point["zenith_angle"] = zenith + point["azimuth"] = azimuth + + self._add_lookup_limits_to_point( + point, + zenith=zenith.value, + azimuth=azimuth.value, + ) + grid_points.append(point) + + return grid_points + def _is_adaptive_horizontal_density_enabled(self): """Return whether horizontal grid generation should use row-wise adaptive azimuth bins.""" return ( diff --git a/src/simtools/production_configuration/simulation_jobs.py b/src/simtools/production_configuration/simulation_jobs.py index 1a3d5f84d1..1c8c467a70 100644 --- a/src/simtools/production_configuration/simulation_jobs.py +++ b/src/simtools/production_configuration/simulation_jobs.py @@ -77,6 +77,10 @@ _HORIZONTAL_AXES = ("azimuth", "zenith") _RADEC_AXES = ("ra", "dec") _REQUIRED_AXES = ("nsb", "offset") +_LOCAL_CONSTRAINT_ARGUMENTS = { + "local_zenith_range": ("zenith", "deg"), + "local_azimuth_range": ("azimuth", "deg"), +} def _parse_axis_range_tokens(range_tokens): @@ -174,6 +178,16 @@ def _resolve_axis_configs(args_dict): return axis_configs +def _parse_optional_range_argument(range_argument_value, default_unit): + """Parse optional quantity-pair CLI/config input into float values.""" + if range_argument_value is None: + return None + + tokens = _normalize_axis_spec_tokens(range_argument_value) + range_quantities = _parse_axis_range_tokens(tokens) + return [u.Quantity(value).to_value(default_unit) for value in range_quantities] + + def _circular_span_degrees(axis_range): """Return directed circular span (degrees) from start to end.""" start, end = axis_range @@ -240,12 +254,17 @@ def _apply_direction_grid_density(axis_configs, direction_axes, density): def _resolve_coordinate_system(axis_configs): - """Resolve the coordinate system from axis definitions.""" + """Resolve the coordinate system from axis definitions. + + Notes + ----- + If RA/Dec axes are present, ``ra_dec`` mode is selected. Optional horizontal + axes (azimuth/zenith) may still be provided and are interpreted as local-sky + constraints for RA/Dec density-based generation. + """ has_horizontal_axes = all(axis_name in axis_configs for axis_name in _HORIZONTAL_AXES) has_radec_axes = all(axis_name in axis_configs for axis_name in _RADEC_AXES) - if has_horizontal_axes and has_radec_axes: - raise ValueError("Cannot define both azimuth/zenith and ra/dec axes at the same time.") if has_radec_axes: return "ra_dec" if has_horizontal_axes: @@ -286,7 +305,12 @@ def build_observing_location(site, model_version): def build_axes_dict_from_cli_args(args_dict): - """Build ProductionGridEngine-compatible axes configuration from CLI arguments.""" + """Build ProductionGridEngine-compatible axes configuration from CLI arguments. + + In ``ra_dec`` mode, optional ``zenith`` and ``azimuth`` axis definitions are + carried through as local-sky constraints. This allows reusing the same + directional range parameters in both horizontal and RA/Dec density workflows. + """ axis_configs = _resolve_axis_configs(args_dict) coordinate_system = _resolve_coordinate_system(axis_configs) @@ -305,14 +329,49 @@ def build_axes_dict_from_cli_args(args_dict): axis_configs["azimuth"]["direction_grid_density"] = float( args_dict["direction_grid_density"] ) + if coordinate_system == "ra_dec" and args_dict.get("direction_grid_density") is not None: + axis_configs["ra"]["direction_grid_density"] = float(args_dict["direction_grid_density"]) _apply_direction_grid_density( axis_configs, direction_axes, args_dict.get("direction_grid_density"), ) + + if coordinate_system == "ra_dec": + local_constraints = {} + for constraint_argument, ( + legacy_axis_name, + default_unit, + ) in _LOCAL_CONSTRAINT_ARGUMENTS.items(): + parsed_argument_range = _parse_optional_range_argument( + args_dict.get(constraint_argument), + default_unit, + ) + legacy_axis = axis_configs.get(legacy_axis_name) + legacy_axis_range = legacy_axis.get("range") if legacy_axis else None + + if parsed_argument_range is not None and legacy_axis_range is not None: + raise ValueError( + f"Cannot define both '{constraint_argument}' and axis '{legacy_axis_name}' " + "for RA/Dec local-sky constraints." + ) + + if parsed_argument_range is not None: + local_constraints[constraint_argument] = parsed_argument_range + elif legacy_axis_range is not None: + logger.warning( + f"Using axis '{legacy_axis_name}' as RA/Dec local-sky constraint is " + f"deprecated. Use '{constraint_argument}' instead." + ) + local_constraints[constraint_argument] = legacy_axis_range + + axis_configs["ra"].update(local_constraints) + + axes_to_export = [*_REQUIRED_AXES, *direction_axes] + return { GRID_AXIS_ARGUMENTS[axis_name]["engine_axis"]: axis_configs[axis_name] - for axis_name in (*_REQUIRED_AXES, *direction_axes) + for axis_name in axes_to_export } @@ -360,6 +419,7 @@ def build_job_grid_metadata(args_dict): "site": args_dict.get("site"), "simulation_software": args_dict.get("simulation_software"), "coordinate_system": coordinate_system, + "direction_grid_density": args_dict.get("direction_grid_density"), "time_of_observation_utc": time_of_observation.isot if time_of_observation else None, "time_of_observation_scale": time_of_observation.scale if time_of_observation else None, "corsika_limits": ( diff --git a/src/simtools/visualization/plot_production_grid.py b/src/simtools/visualization/plot_production_grid.py index 8075d37785..3fa08f0b9b 100644 --- a/src/simtools/visualization/plot_production_grid.py +++ b/src/simtools/visualization/plot_production_grid.py @@ -299,6 +299,10 @@ def plot_sky_projection(self, plot_ra_dec_tracks=False, dec_values=None): subtitle_lines = [] if self.grid_metadata.get("site"): subtitle_lines.append(f"Site: {self.grid_metadata['site']}") + if self.grid_metadata.get("direction_grid_density") is not None: + subtitle_lines.append( + f"Grid density: {self.grid_metadata['direction_grid_density']} nodes/deg^2" + ) if self.grid_metadata.get("time_of_observation_utc"): subtitle_lines.append( f"Observation time (UTC): {self.grid_metadata['time_of_observation_utc']}" @@ -331,7 +335,6 @@ def _configure_altaz_axis(self, axis): axis.set_rmax(90) axis.set_rticks(np.arange(10, 91, 10)) axis.grid(True, color="gray", alpha=0.5, linestyle="--") - axis.text(0, 0, "Zenith", ha="center", va="center", fontsize=11, fontweight="bold") axis.set_title("Local Alt/Az", pad=18) def _configure_radec_axis(self, axis, plot_points): diff --git a/tests/integration_tests/config/production_generate_grid_horizontal_density.yml b/tests/integration_tests/config/production_generate_grid_horizontal_density.yml new file mode 100644 index 0000000000..79bf8565ab --- /dev/null +++ b/tests/integration_tests/config/production_generate_grid_horizontal_density.yml @@ -0,0 +1,35 @@ +--- +applications: +- application: simtools-production-generate-grid + configuration: + array_layout_name: + by_version: + "<7.0.0": 1lst + ">=7.0.0": LSTN-01 + axis: + - azimuth 0 deg 360 deg 0 linear + - zenith 0 deg 70 deg 0 linear + - nsb 4 MHz 5 MHz 2 linear + - offset 0 deg 10 deg 2 linear + direction_grid_density: 0.25 + core_scatter: 10 500 m + corsika_he_interaction: epos + corsika_le_interaction: urqmd + energy_range: 30 GeV 300 GeV + corsika_limits: tests/resources/corsika_simulation_limits/corsika_limits_for_test.ecsv + model_version: 7.0.0 + showers_per_run: 5 + showers_per_run_power_law: 0.0 1 TeV + number_of_runs: 2 + output_file: job_grid_horizontal_grid_density.ecsv + output_path: simtools-output + primary: gamma + simulation_software: corsika_sim_telarray + site: North + view_cone: 0 deg 10 deg + integration_tests: + - output_file: job_grid_horizontal_density.ecsv + model_version_use_current: true + test_name: production_generate_grid_horizontal +schema_name: application_workflow.metaschema +schema_version: 0.4.0 diff --git a/tests/integration_tests/config/production_generate_grid_ra_dec_density.yml b/tests/integration_tests/config/production_generate_grid_ra_dec_density.yml new file mode 100644 index 0000000000..f9e958d5b6 --- /dev/null +++ b/tests/integration_tests/config/production_generate_grid_ra_dec_density.yml @@ -0,0 +1,42 @@ +--- +applications: +- application: simtools-production-generate-grid + configuration: + array_layout_name: + by_version: + "<7.0.0": 1lst + ">=7.0.0": LSTN-01 + axis: + # RA/Dec define the sampling frame; binning placeholders are overridden by + # direction_grid_density in adaptive RA/Dec mode. + - ra 0 deg 360 deg 1 linear + - dec -40 deg 80 deg 1 linear + - nsb 4 MHz 4 MHz 1 linear + - offset 0 deg 10 deg 2 linear + # Explicit local-sky constraints for RA/Dec density filtering. + local_zenith_range: 0 deg 70 deg + # Directed azimuth interval (300 -> 60) across wrap-around. + local_azimuth_range: 300 deg 60 deg + direction_grid_density: 0.25 + core_scatter: 10 500 m + corsika_he_interaction: epos + corsika_le_interaction: urqmd + energy_range: 30 GeV 300 GeV + corsika_limits: tests/resources/corsika_simulation_limits/corsika_limits_for_test.ecsv + model_version: 7.0.0 + showers_per_run: 5 + showers_per_run_power_law: 0.0 1 TeV + number_of_runs: 2 + time_of_observation: '2017-09-16 00:00:00' + output_file: job_grid_ra_dec_density.ecsv + output_path: simtools-output + primary: gamma + simulation_software: corsika_sim_telarray + site: North + view_cone: 0 deg 10 deg + integration_tests: + - output_file: job_grid_ra_dec_density.ecsv + skip_integration_test: Requires model DB access for site-location lookup in ra_dec mode. + test_name: production_generate_grid_ra_dec +schema_name: application_workflow.metaschema +schema_version: 0.4.0 diff --git a/tests/unit_tests/production_configuration/test_observation_grid.py b/tests/unit_tests/production_configuration/test_observation_grid.py index c2b0dd191d..6e182cd641 100644 --- a/tests/unit_tests/production_configuration/test_observation_grid.py +++ b/tests/unit_tests/production_configuration/test_observation_grid.py @@ -203,6 +203,123 @@ def test_generate_horizontal_grid_uses_adaptive_azimuth_bins_per_zenith_row(): assert len(azimuth_counts_by_zenith[60.0]) > len(azimuth_counts_by_zenith[30.0]) +def test_generate_radec_grid_uses_adaptive_ra_bins_per_dec_strip(): + engine = ProductionGridEngine( + axes={ + "ra": {"range": [0, 360], "binning": 36, "units": "deg", "direction_grid_density": 1.0}, + "dec": {"range": [-60, 60], "binning": 13, "units": "deg"}, + "nsb_level": {"range": [4, 4], "binning": 1, "units": "MHz"}, + "offset": {"range": [0, 0], "binning": 1, "units": "deg"}, + }, + coordinate_system="ra_dec", + time_of_observation=Time("2020-01-01 00:00:00", scale="utc"), + ) + + assert engine._is_adaptive_radec_density_enabled() + grid = engine._generate_adaptive_radec_grid(include_horizontal_coordinates=False) + + assert len(grid) > 0 + + ra_counts_by_dec = {} + for point in grid: + dec_val = round(point["dec"].to_value(u.deg), 6) + ra_counts_by_dec.setdefault(dec_val, 0) + ra_counts_by_dec[dec_val] += 1 + + # Near equator (dec=0) should have more RA points than near the pole (dec=60 deg) + assert ra_counts_by_dec[0.0] > ra_counts_by_dec[60.0] + + +def test_generate_radec_grid_adaptive_density_keeps_only_visible_nodes(): + engine = ProductionGridEngine( + axes={ + "ra": { + "range": [0, 360], + "binning": 36, + "units": "deg", + "direction_grid_density": 0.25, + }, + "dec": {"range": [-40, 80], "binning": 10, "units": "deg"}, + "nsb_level": {"range": [4, 4], "binning": 1, "units": "MHz"}, + "offset": {"range": [0, 10], "binning": 2, "units": "deg"}, + }, + coordinate_system="ra_dec", + observing_location=EarthLocation( + lat=28.76 * u.deg, + lon=-17.89 * u.deg, + height=2200 * u.m, + ), + time_of_observation=Time("2020-01-01 00:00:00", scale="utc"), + ) + + grid = engine._generate_adaptive_radec_grid(include_horizontal_coordinates=True) + + assert len(grid) > 0 + assert all(point["zenith_angle"].to_value(u.deg) <= 90.0 for point in grid) + + +def test_generate_radec_grid_adaptive_density_applies_zenith_constraint(): + engine = ProductionGridEngine( + axes={ + "ra": { + "range": [0, 360], + "binning": 36, + "units": "deg", + "direction_grid_density": 0.25, + "local_zenith_range": [0, 70], + }, + "dec": {"range": [-40, 80], "binning": 10, "units": "deg"}, + "nsb_level": {"range": [4, 4], "binning": 1, "units": "MHz"}, + "offset": {"range": [0, 10], "binning": 2, "units": "deg"}, + }, + coordinate_system="ra_dec", + observing_location=EarthLocation( + lat=28.76 * u.deg, + lon=-17.89 * u.deg, + height=2200 * u.m, + ), + time_of_observation=Time("2020-01-01 00:00:00", scale="utc"), + ) + + grid = engine._generate_adaptive_radec_grid(include_horizontal_coordinates=True) + + assert len(grid) > 0 + assert all(point["zenith_angle"].to_value(u.deg) <= 70.0 for point in grid) + + +def test_generate_radec_grid_adaptive_density_applies_azimuth_constraint(): + engine = ProductionGridEngine( + axes={ + "ra": { + "range": [0, 360], + "binning": 36, + "units": "deg", + "direction_grid_density": 0.25, + "local_zenith_range": [0, 70], + "local_azimuth_range": [300, 60], + }, + "dec": {"range": [-40, 80], "binning": 10, "units": "deg"}, + "nsb_level": {"range": [4, 4], "binning": 1, "units": "MHz"}, + "offset": {"range": [0, 10], "binning": 2, "units": "deg"}, + }, + coordinate_system="ra_dec", + observing_location=EarthLocation( + lat=28.76 * u.deg, + lon=-17.89 * u.deg, + height=2200 * u.m, + ), + time_of_observation=Time("2020-01-01 00:00:00", scale="utc"), + ) + + grid = engine._generate_adaptive_radec_grid(include_horizontal_coordinates=True) + + assert len(grid) > 0 + assert all( + (point["azimuth"].to_value(u.deg) >= 300.0) or (point["azimuth"].to_value(u.deg) <= 60.0) + for point in grid + ) + + def test_create_circular_binning_treats_full_circle_range_as_full_span(): engine = ProductionGridEngine(axes={"axes": {}}) diff --git a/tests/unit_tests/production_configuration/test_simulation_jobs.py b/tests/unit_tests/production_configuration/test_simulation_jobs.py index 570ce2ee3a..e6d0502e0b 100644 --- a/tests/unit_tests/production_configuration/test_simulation_jobs.py +++ b/tests/unit_tests/production_configuration/test_simulation_jobs.py @@ -123,6 +123,7 @@ def test_build_axes_dict_from_cli_args_derives_radec_binning_from_density(): assert axes["ra"]["binning"] == 230 assert axes["dec"]["binning"] == 180 + assert axes["ra"]["direction_grid_density"] == pytest.approx(1.0) def test_build_axes_dict_from_cli_args_reduces_ra_binning_towards_dec_poles(): @@ -140,6 +141,46 @@ def test_build_axes_dict_from_cli_args_reduces_ra_binning_towards_dec_poles(): assert axes["ra"]["binning"] == 32 assert axes["dec"]["binning"] == 10 + assert axes["ra"]["direction_grid_density"] == pytest.approx(1.0) + + +def test_build_axes_dict_from_cli_args_sets_radec_density_metadata_for_adaptive_grid(): + axes = build_axes_dict_from_cli_args( + { + "direction_grid_density": 0.5, + "axis": [ + ["ra", "0", "deg", "360", "deg", "36", "linear"], + ["dec", "-30", "deg", "30", "deg", "6", "linear"], + ["nsb", "4", "MHz", "4", "MHz", "1", "linear"], + ["offset", "0", "deg", "10", "deg", "2", "linear"], + ], + } + ) + + assert axes["ra"]["direction_grid_density"] == pytest.approx(0.5) + + +def test_build_axes_dict_from_cli_args_keeps_horizontal_constraints_in_radec_mode(): + axes = build_axes_dict_from_cli_args( + { + "direction_grid_density": 0.25, + "local_zenith_range": ["0", "deg", "70", "deg"], + "local_azimuth_range": ["300", "deg", "60", "deg"], + "axis": [ + ["ra", "0", "deg", "360", "deg", "36", "linear"], + ["dec", "-40", "deg", "80", "deg", "10", "linear"], + ["nsb", "4", "MHz", "4", "MHz", "1", "linear"], + ["offset", "0", "deg", "10", "deg", "2", "linear"], + ], + } + ) + + assert "ra" in axes + assert "dec" in axes + assert "zenith_angle" not in axes + assert "azimuth" not in axes + assert axes["ra"]["local_zenith_range"] == pytest.approx([0.0, 70.0]) + assert axes["ra"]["local_azimuth_range"] == pytest.approx([300.0, 60.0]) def test_build_axes_dict_from_cli_args_reduces_az_binning_towards_zenith_pole(): @@ -260,18 +301,40 @@ def test_build_axes_dict_from_cli_args_accepts_merged_config_axis_list(): assert axes["offset"]["range"] == [0.0, 10.0] -def test_build_axes_dict_from_cli_args_rejects_multiple_direction_coordinate_systems(): - with pytest.raises(ValueError, match="Cannot define both azimuth/zenith and ra/dec axes"): +def test_build_axes_dict_from_cli_args_prefers_radec_with_horizontal_constraints(): + axes = build_axes_dict_from_cli_args( + { + "axis": [ + ["azimuth", "310", "deg", "20", "deg", "3"], + ["zenith", "30", "deg", "40", "deg", "2"], + ["ra", "0", "deg", "360", "deg", "36"], + ["dec", "-90", "deg", "90", "deg", "18"], + ["nsb", "4", "MHz", "5", "MHz", "2"], + ["offset", "0", "deg", "10", "deg", "2"], + ] + } + ) + + assert "ra" in axes + assert "dec" in axes + assert "azimuth" not in axes + assert "zenith_angle" not in axes + assert axes["ra"]["local_azimuth_range"] == pytest.approx([310.0, 20.0]) + assert axes["ra"]["local_zenith_range"] == pytest.approx([30.0, 40.0]) + + +def test_build_axes_dict_from_cli_args_rejects_duplicate_local_constraints(): + with pytest.raises(ValueError, match="Cannot define both 'local_zenith_range'"): build_axes_dict_from_cli_args( { + "local_zenith_range": ["0", "deg", "70", "deg"], "axis": [ - ["azimuth", "310", "deg", "20", "deg", "3"], - ["zenith", "30", "deg", "40", "deg", "2"], ["ra", "0", "deg", "360", "deg", "36"], ["dec", "-90", "deg", "90", "deg", "18"], + ["zenith", "0", "deg", "70", "deg", "2"], ["nsb", "4", "MHz", "5", "MHz", "2"], ["offset", "0", "deg", "10", "deg", "2"], - ] + ], } ) @@ -289,6 +352,7 @@ def test_build_job_grid_metadata_includes_job_context(): { "site": "North", "simulation_software": "corsika_sim_telarray", + "direction_grid_density": 0.25, "axis": [ ["ra", "0", "deg", "1", "deg", "2"], ["dec", "0", "deg", "1", "deg", "2"], @@ -301,6 +365,7 @@ def test_build_job_grid_metadata_includes_job_context(): assert metadata["site"] == "North" assert metadata["simulation_software"] == "corsika_sim_telarray" assert metadata["coordinate_system"] == "ra_dec" + assert metadata["direction_grid_density"] == pytest.approx(0.25) assert metadata["time_of_observation_utc"].startswith("2017-09-16T00:00:00") assert metadata["corsika_limits"] == "limits.ecsv" diff --git a/tests/unit_tests/visualization/test_plot_production_grid.py b/tests/unit_tests/visualization/test_plot_production_grid.py index 3666477c48..a43d0d5773 100644 --- a/tests/unit_tests/visualization/test_plot_production_grid.py +++ b/tests/unit_tests/visualization/test_plot_production_grid.py @@ -328,3 +328,33 @@ def test_plot_sky_projection_logs_tracks_disabled(tmp_test_directory, caplog): plotter.plot_sky_projection(plot_ra_dec_tracks=True) assert "RA/Dec tracks are disabled in file-driven plotting mode" in caplog.text + + +def test_plot_sky_projection_writes_grid_density_subtitle(tmp_test_directory, monkeypatch): + """Render direction grid density in subtitle when present in metadata.""" + grid_file = _write_grid_file( + tmp_test_directory, + "grid_with_density_meta.ecsv", + [{"azimuth": 100.0, "zenith_angle": 25.0, "ra": 180.0, "dec": -10.0}], + ) + table = Table.read(grid_file, format="ascii.ecsv") + table.meta["direction_grid_density"] = 0.25 + table.write(grid_file, format="ascii.ecsv", overwrite=True) + + plotter = _create_plotter( + grid_file=grid_file, + output_path=Path(tmp_test_directory) / "output", + ) + + recorded_text = [] + original_text = plt.Figure.text + + def _record_text(self, x, y, s, **kwargs): + recorded_text.append(s) + return original_text(self, x, y, s, **kwargs) + + monkeypatch.setattr(plt.Figure, "text", _record_text) + + plotter.plot_sky_projection() + + assert any(text == "Grid density: 0.25 nodes/deg^2" for text in recorded_text) From ff01f995e1abd5e2e4b94639104b00b106f211d9 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Thu, 28 May 2026 11:29:08 +0200 Subject: [PATCH 03/10] grid density units --- .../applications/production_generate_grid.py | 10 ++--- .../simulation_jobs.py | 41 +++++++++++++++---- .../visualization/plot_production_grid.py | 3 +- ...ction_generate_grid_horizontal_density.yml | 2 +- ...roduction_generate_grid_ra_dec_density.yml | 2 +- .../test_production_generate_grid.py | 11 ++++- .../test_simulation_jobs.py | 17 ++++++++ .../test_plot_production_grid.py | 3 +- 8 files changed, 71 insertions(+), 18 deletions(-) diff --git a/src/simtools/applications/production_generate_grid.py b/src/simtools/applications/production_generate_grid.py index a6074e31ae..29b85294a1 100644 --- a/src/simtools/applications/production_generate_grid.py +++ b/src/simtools/applications/production_generate_grid.py @@ -22,8 +22,8 @@ ``--axis [scaling]``. Example: ``--axis azimuth 310 deg 20 deg 3 linear``. Options for scaling are: linear, log, 1/cos. -direction_grid_density (float, optional) - Direction-grid density in points per deg^2. +direction_grid_density (float or quantity, optional) + Direction-grid density in ``1/deg^2``. If provided, direction-axis binning (azimuth/zenith or ra/dec) is derived from range and density, while min/max values are kept from ``--axis`` definitions. In ``ra_dec`` mode, local-sky constraints can be defined with @@ -89,7 +89,7 @@ --axis dec -40 deg 80 deg 1 linear \ --axis nsb 4 MHz 4 MHz 1 linear \ --axis offset 0 deg 10 deg 2 linear \ - --direction_grid_density 0.25 \ + --direction_grid_density 0.25 1/deg^2 \ --local_zenith_range 0 deg 70 deg \ --local_azimuth_range 300 deg 60 deg \ --time_of_observation "2017-09-16 00:00:00" \ @@ -128,11 +128,11 @@ def _add_arguments(parser): ) parser.add_argument( "--direction_grid_density", - type=float, + nargs="+", required=False, default=None, help=( - "Direction-grid density in points per deg^2. If set, direction-axis binning is " + "Direction-grid density in 1/deg^2. If set, direction-axis binning is " "derived from axis ranges and this density. In ra_dec mode, optional zenith and " "azimuth axes are interpreted as local-sky constraints for filtering generated " "points (for example zenith 0..70 deg)." diff --git a/src/simtools/production_configuration/simulation_jobs.py b/src/simtools/production_configuration/simulation_jobs.py index 1c8c467a70..bf4b17ed20 100644 --- a/src/simtools/production_configuration/simulation_jobs.py +++ b/src/simtools/production_configuration/simulation_jobs.py @@ -81,6 +81,7 @@ "local_zenith_range": ("zenith", "deg"), "local_azimuth_range": ("azimuth", "deg"), } +_DIRECTION_GRID_DENSITY_UNIT = 1 / u.deg**2 def _parse_axis_range_tokens(range_tokens): @@ -188,6 +189,29 @@ def _parse_optional_range_argument(range_argument_value, default_unit): return [u.Quantity(value).to_value(default_unit) for value in range_quantities] +def _parse_direction_grid_density(density_value): + """Parse direction-grid density and normalize it to ``1/deg^2`` units.""" + if density_value is None: + return None + + if isinstance(density_value, (int, float)): + return float(density_value) + + if isinstance(density_value, (str, list, tuple)): + tokens = _normalize_axis_spec_tokens(density_value) + if len(tokens) == 1: + return float(tokens[0]) + try: + quantity = u.Quantity(" ".join(tokens)) + return quantity.to_value(_DIRECTION_GRID_DENSITY_UNIT) + except (TypeError, ValueError, u.UnitConversionError) as exc: + raise ValueError( + "direction_grid_density must be a float or quantity in 1/deg^2." + ) from exc + + raise TypeError("direction_grid_density must be a number, string, or CLI-style token list.") + + def _circular_span_degrees(axis_range): """Return directed circular span (degrees) from start to end.""" start, end = axis_range @@ -324,17 +348,16 @@ def build_axes_dict_from_cli_args(args_dict): missing_axes = ", ".join(missing_required_axes) raise ValueError(f"Missing required shared axis definition(s): {missing_axes}.") + direction_grid_density = _parse_direction_grid_density(args_dict.get("direction_grid_density")) direction_axes = _RADEC_AXES if coordinate_system == "ra_dec" else _HORIZONTAL_AXES - if coordinate_system == "horizontal" and args_dict.get("direction_grid_density") is not None: - axis_configs["azimuth"]["direction_grid_density"] = float( - args_dict["direction_grid_density"] - ) - if coordinate_system == "ra_dec" and args_dict.get("direction_grid_density") is not None: - axis_configs["ra"]["direction_grid_density"] = float(args_dict["direction_grid_density"]) + if coordinate_system == "horizontal" and direction_grid_density is not None: + axis_configs["azimuth"]["direction_grid_density"] = direction_grid_density + if coordinate_system == "ra_dec" and direction_grid_density is not None: + axis_configs["ra"]["direction_grid_density"] = direction_grid_density _apply_direction_grid_density( axis_configs, direction_axes, - args_dict.get("direction_grid_density"), + direction_grid_density, ) if coordinate_system == "ra_dec": @@ -413,13 +436,15 @@ def build_job_grid_metadata(args_dict): args_dict, ) coordinate_system = _resolve_coordinate_system_from_args(args_dict) + direction_grid_density = _parse_direction_grid_density(args_dict.get("direction_grid_density")) if coordinate_system == "ra_dec" and not args_dict.get("site"): raise ValueError("site is required when using RA/Dec axes.") return { "site": args_dict.get("site"), "simulation_software": args_dict.get("simulation_software"), "coordinate_system": coordinate_system, - "direction_grid_density": args_dict.get("direction_grid_density"), + "direction_grid_density": direction_grid_density, + "direction_grid_density_unit": "1/deg^2" if direction_grid_density is not None else None, "time_of_observation_utc": time_of_observation.isot if time_of_observation else None, "time_of_observation_scale": time_of_observation.scale if time_of_observation else None, "corsika_limits": ( diff --git a/src/simtools/visualization/plot_production_grid.py b/src/simtools/visualization/plot_production_grid.py index 3fa08f0b9b..1ba1592249 100644 --- a/src/simtools/visualization/plot_production_grid.py +++ b/src/simtools/visualization/plot_production_grid.py @@ -300,8 +300,9 @@ def plot_sky_projection(self, plot_ra_dec_tracks=False, dec_values=None): if self.grid_metadata.get("site"): subtitle_lines.append(f"Site: {self.grid_metadata['site']}") if self.grid_metadata.get("direction_grid_density") is not None: + density_unit = self.grid_metadata.get("direction_grid_density_unit") or "1/deg^2" subtitle_lines.append( - f"Grid density: {self.grid_metadata['direction_grid_density']} nodes/deg^2" + f"Grid density: {self.grid_metadata['direction_grid_density']} {density_unit}" ) if self.grid_metadata.get("time_of_observation_utc"): subtitle_lines.append( diff --git a/tests/integration_tests/config/production_generate_grid_horizontal_density.yml b/tests/integration_tests/config/production_generate_grid_horizontal_density.yml index 79bf8565ab..2d9608b60e 100644 --- a/tests/integration_tests/config/production_generate_grid_horizontal_density.yml +++ b/tests/integration_tests/config/production_generate_grid_horizontal_density.yml @@ -11,7 +11,7 @@ applications: - zenith 0 deg 70 deg 0 linear - nsb 4 MHz 5 MHz 2 linear - offset 0 deg 10 deg 2 linear - direction_grid_density: 0.25 + direction_grid_density: 0.25 1/deg^2 core_scatter: 10 500 m corsika_he_interaction: epos corsika_le_interaction: urqmd diff --git a/tests/integration_tests/config/production_generate_grid_ra_dec_density.yml b/tests/integration_tests/config/production_generate_grid_ra_dec_density.yml index f9e958d5b6..ff5e47c0ca 100644 --- a/tests/integration_tests/config/production_generate_grid_ra_dec_density.yml +++ b/tests/integration_tests/config/production_generate_grid_ra_dec_density.yml @@ -17,7 +17,7 @@ applications: local_zenith_range: 0 deg 70 deg # Directed azimuth interval (300 -> 60) across wrap-around. local_azimuth_range: 300 deg 60 deg - direction_grid_density: 0.25 + direction_grid_density: 0.25 1/deg^2 core_scatter: 10 500 m corsika_he_interaction: epos corsika_le_interaction: urqmd diff --git a/tests/unit_tests/applications/test_production_generate_grid.py b/tests/unit_tests/applications/test_production_generate_grid.py index 4cf70cda50..ea96a4fe47 100644 --- a/tests/unit_tests/applications/test_production_generate_grid.py +++ b/tests/unit_tests/applications/test_production_generate_grid.py @@ -83,4 +83,13 @@ def test_add_arguments_accepts_direction_grid_density(): args = parser.parse_args(["--direction_grid_density", "1.5"]) - assert args.direction_grid_density == pytest.approx(1.5) + assert args.direction_grid_density == ["1.5"] + + +def test_add_arguments_accepts_direction_grid_density_with_unit(): + parser = CommandLineParser() + app._add_arguments(parser) + + args = parser.parse_args(["--direction_grid_density", "0.25", "1/deg^2"]) + + assert args.direction_grid_density == ["0.25", "1/deg^2"] diff --git a/tests/unit_tests/production_configuration/test_simulation_jobs.py b/tests/unit_tests/production_configuration/test_simulation_jobs.py index e6d0502e0b..0b6d9000d7 100644 --- a/tests/unit_tests/production_configuration/test_simulation_jobs.py +++ b/tests/unit_tests/production_configuration/test_simulation_jobs.py @@ -108,6 +108,22 @@ def test_build_axes_dict_from_cli_args_derives_horizontal_binning_from_density() assert axes["zenith_angle"]["binning"] == 10 +def test_build_axes_dict_from_cli_args_accepts_density_with_unit(): + axes = build_axes_dict_from_cli_args( + { + "direction_grid_density": "1.0 1/deg^2", + "axis": [ + ["azimuth", "310", "deg", "20", "deg", "3", "linear"], + ["zenith", "30", "deg", "40", "deg", "2"], + ["nsb", "4", "MHz", "5", "MHz", "2"], + ["offset", "0", "deg", "10", "deg", "2"], + ], + } + ) + + assert axes["azimuth"]["direction_grid_density"] == pytest.approx(1.0) + + def test_build_axes_dict_from_cli_args_derives_radec_binning_from_density(): axes = build_axes_dict_from_cli_args( { @@ -366,6 +382,7 @@ def test_build_job_grid_metadata_includes_job_context(): assert metadata["simulation_software"] == "corsika_sim_telarray" assert metadata["coordinate_system"] == "ra_dec" assert metadata["direction_grid_density"] == pytest.approx(0.25) + assert metadata["direction_grid_density_unit"] == "1/deg^2" assert metadata["time_of_observation_utc"].startswith("2017-09-16T00:00:00") assert metadata["corsika_limits"] == "limits.ecsv" diff --git a/tests/unit_tests/visualization/test_plot_production_grid.py b/tests/unit_tests/visualization/test_plot_production_grid.py index a43d0d5773..b463e194d3 100644 --- a/tests/unit_tests/visualization/test_plot_production_grid.py +++ b/tests/unit_tests/visualization/test_plot_production_grid.py @@ -339,6 +339,7 @@ def test_plot_sky_projection_writes_grid_density_subtitle(tmp_test_directory, mo ) table = Table.read(grid_file, format="ascii.ecsv") table.meta["direction_grid_density"] = 0.25 + table.meta["direction_grid_density_unit"] = "1/deg^2" table.write(grid_file, format="ascii.ecsv", overwrite=True) plotter = _create_plotter( @@ -357,4 +358,4 @@ def _record_text(self, x, y, s, **kwargs): plotter.plot_sky_projection() - assert any(text == "Grid density: 0.25 nodes/deg^2" for text in recorded_text) + assert any(text == "Grid density: 0.25 1/deg^2" for text in recorded_text) From e744a8f942411fcd194f075a007cb6fc82cd4b46 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Thu, 28 May 2026 11:45:26 +0200 Subject: [PATCH 04/10] unit tests --- .../applications/production_generate_grid.py | 5 +- .../observation_grid.py | 10 --- .../simulation_jobs.py | 42 ++--------- .../test_observation_grid.py | 34 +++++++++ .../test_simulation_jobs.py | 70 +++++++++++++------ 5 files changed, 90 insertions(+), 71 deletions(-) diff --git a/src/simtools/applications/production_generate_grid.py b/src/simtools/applications/production_generate_grid.py index 29b85294a1..e75331d735 100644 --- a/src/simtools/applications/production_generate_grid.py +++ b/src/simtools/applications/production_generate_grid.py @@ -133,9 +133,8 @@ def _add_arguments(parser): default=None, help=( "Direction-grid density in 1/deg^2. If set, direction-axis binning is " - "derived from axis ranges and this density. In ra_dec mode, optional zenith and " - "azimuth axes are interpreted as local-sky constraints for filtering generated " - "points (for example zenith 0..70 deg)." + "derived from axis ranges and this density. In ra_dec mode, use " + "local_zenith_range/local_azimuth_range to filter generated points in local sky." ), ) parser.add_argument( diff --git a/src/simtools/production_configuration/observation_grid.py b/src/simtools/production_configuration/observation_grid.py index 09d4647187..cd08395aa9 100644 --- a/src/simtools/production_configuration/observation_grid.py +++ b/src/simtools/production_configuration/observation_grid.py @@ -388,11 +388,6 @@ def _compute_visible_radec_strip(self, dec_deg, spacing, lst_deg, time_of_observ ra_axis = self.axes.get("ra", {}) zenith_range = ra_axis.get("local_zenith_range") - if zenith_range is None: - zenith_axis = self.axes.get("zenith_angle") - if zenith_axis and "range" in zenith_axis: - zenith_range = zenith_axis["range"] - if zenith_range is not None: zenith_min, zenith_max = sorted( ( @@ -404,11 +399,6 @@ def _compute_visible_radec_strip(self, dec_deg, spacing, lst_deg, time_of_observ visible_mask &= (zenith_values_deg >= zenith_min) & (zenith_values_deg <= zenith_max) azimuth_range = ra_axis.get("local_azimuth_range") - if azimuth_range is None: - azimuth_axis = self.axes.get("azimuth") - if azimuth_axis and "range" in azimuth_axis: - azimuth_range = azimuth_axis["range"] - if azimuth_range is not None: azimuth_values_deg = azimuth_values.to_value(u.deg) visible_mask &= self._is_in_directed_azimuth_range( diff --git a/src/simtools/production_configuration/simulation_jobs.py b/src/simtools/production_configuration/simulation_jobs.py index bf4b17ed20..97c774b48a 100644 --- a/src/simtools/production_configuration/simulation_jobs.py +++ b/src/simtools/production_configuration/simulation_jobs.py @@ -78,8 +78,8 @@ _RADEC_AXES = ("ra", "dec") _REQUIRED_AXES = ("nsb", "offset") _LOCAL_CONSTRAINT_ARGUMENTS = { - "local_zenith_range": ("zenith", "deg"), - "local_azimuth_range": ("azimuth", "deg"), + "local_zenith_range": "deg", + "local_azimuth_range": "deg", } _DIRECTION_GRID_DENSITY_UNIT = 1 / u.deg**2 @@ -278,17 +278,12 @@ def _apply_direction_grid_density(axis_configs, direction_axes, density): def _resolve_coordinate_system(axis_configs): - """Resolve the coordinate system from axis definitions. - - Notes - ----- - If RA/Dec axes are present, ``ra_dec`` mode is selected. Optional horizontal - axes (azimuth/zenith) may still be provided and are interpreted as local-sky - constraints for RA/Dec density-based generation. - """ + """Resolve the coordinate system from axis definitions.""" has_horizontal_axes = all(axis_name in axis_configs for axis_name in _HORIZONTAL_AXES) has_radec_axes = all(axis_name in axis_configs for axis_name in _RADEC_AXES) + if has_horizontal_axes and has_radec_axes: + raise ValueError("Cannot define both azimuth/zenith and ra/dec axes at the same time.") if has_radec_axes: return "ra_dec" if has_horizontal_axes: @@ -329,12 +324,7 @@ def build_observing_location(site, model_version): def build_axes_dict_from_cli_args(args_dict): - """Build ProductionGridEngine-compatible axes configuration from CLI arguments. - - In ``ra_dec`` mode, optional ``zenith`` and ``azimuth`` axis definitions are - carried through as local-sky constraints. This allows reusing the same - directional range parameters in both horizontal and RA/Dec density workflows. - """ + """Build ProductionGridEngine-compatible axes configuration from CLI arguments.""" axis_configs = _resolve_axis_configs(args_dict) coordinate_system = _resolve_coordinate_system(axis_configs) @@ -362,31 +352,13 @@ def build_axes_dict_from_cli_args(args_dict): if coordinate_system == "ra_dec": local_constraints = {} - for constraint_argument, ( - legacy_axis_name, - default_unit, - ) in _LOCAL_CONSTRAINT_ARGUMENTS.items(): + for constraint_argument, default_unit in _LOCAL_CONSTRAINT_ARGUMENTS.items(): parsed_argument_range = _parse_optional_range_argument( args_dict.get(constraint_argument), default_unit, ) - legacy_axis = axis_configs.get(legacy_axis_name) - legacy_axis_range = legacy_axis.get("range") if legacy_axis else None - - if parsed_argument_range is not None and legacy_axis_range is not None: - raise ValueError( - f"Cannot define both '{constraint_argument}' and axis '{legacy_axis_name}' " - "for RA/Dec local-sky constraints." - ) - if parsed_argument_range is not None: local_constraints[constraint_argument] = parsed_argument_range - elif legacy_axis_range is not None: - logger.warning( - f"Using axis '{legacy_axis_name}' as RA/Dec local-sky constraint is " - f"deprecated. Use '{constraint_argument}' instead." - ) - local_constraints[constraint_argument] = legacy_axis_range axis_configs["ra"].update(local_constraints) diff --git a/tests/unit_tests/production_configuration/test_observation_grid.py b/tests/unit_tests/production_configuration/test_observation_grid.py index 6e182cd641..93b323d117 100644 --- a/tests/unit_tests/production_configuration/test_observation_grid.py +++ b/tests/unit_tests/production_configuration/test_observation_grid.py @@ -230,6 +230,31 @@ def test_generate_radec_grid_uses_adaptive_ra_bins_per_dec_strip(): assert ra_counts_by_dec[0.0] > ra_counts_by_dec[60.0] +def test_generate_grid_from_radec_axes_routes_to_adaptive_density_path(): + engine = ProductionGridEngine( + axes={ + "ra": { + "range": [0, 360], + "binning": 36, + "units": "deg", + "direction_grid_density": 0.25, + }, + "dec": {"range": [-40, 80], "binning": 10, "units": "deg"}, + "nsb_level": {"range": [4, 4], "binning": 1, "units": "MHz"}, + "offset": {"range": [0, 10], "binning": 2, "units": "deg"}, + }, + coordinate_system="ra_dec", + time_of_observation=Time("2020-01-01 00:00:00", scale="utc"), + ) + expected_grid = [{"ra": 1 * u.deg, "dec": 2 * u.deg}] + engine._generate_adaptive_radec_grid = Mock(return_value=expected_grid) + + grid = engine._generate_grid_from_radec_axes(include_horizontal_coordinates=True) + + assert grid == expected_grid + engine._generate_adaptive_radec_grid.assert_called_once_with(True) + + def test_generate_radec_grid_adaptive_density_keeps_only_visible_nodes(): engine = ProductionGridEngine( axes={ @@ -320,6 +345,15 @@ def test_generate_radec_grid_adaptive_density_applies_azimuth_constraint(): ) +def test_is_in_directed_azimuth_range_returns_all_true_for_full_circle(): + mask = ProductionGridEngine._is_in_directed_azimuth_range( + azimuth_values_deg=np.array([0.0, 90.0, 180.0, 270.0]), + azimuth_range=(0.0, 360.0), + ) + + assert np.all(mask) + + def test_create_circular_binning_treats_full_circle_range_as_full_span(): engine = ProductionGridEngine(axes={"axes": {}}) diff --git a/tests/unit_tests/production_configuration/test_simulation_jobs.py b/tests/unit_tests/production_configuration/test_simulation_jobs.py index 0b6d9000d7..4916d13e8e 100644 --- a/tests/unit_tests/production_configuration/test_simulation_jobs.py +++ b/tests/unit_tests/production_configuration/test_simulation_jobs.py @@ -124,6 +124,22 @@ def test_build_axes_dict_from_cli_args_accepts_density_with_unit(): assert axes["azimuth"]["direction_grid_density"] == pytest.approx(1.0) +def test_build_axes_dict_from_cli_args_accepts_density_as_single_token_list(): + axes = build_axes_dict_from_cli_args( + { + "direction_grid_density": ["1.0"], + "axis": [ + ["azimuth", "310", "deg", "20", "deg", "3", "linear"], + ["zenith", "30", "deg", "40", "deg", "2"], + ["nsb", "4", "MHz", "5", "MHz", "2"], + ["offset", "0", "deg", "10", "deg", "2"], + ], + } + ) + + assert axes["azimuth"]["direction_grid_density"] == pytest.approx(1.0) + + def test_build_axes_dict_from_cli_args_derives_radec_binning_from_density(): axes = build_axes_dict_from_cli_args( { @@ -317,37 +333,45 @@ def test_build_axes_dict_from_cli_args_accepts_merged_config_axis_list(): assert axes["offset"]["range"] == [0.0, 10.0] -def test_build_axes_dict_from_cli_args_prefers_radec_with_horizontal_constraints(): - axes = build_axes_dict_from_cli_args( - { - "axis": [ - ["azimuth", "310", "deg", "20", "deg", "3"], - ["zenith", "30", "deg", "40", "deg", "2"], - ["ra", "0", "deg", "360", "deg", "36"], - ["dec", "-90", "deg", "90", "deg", "18"], - ["nsb", "4", "MHz", "5", "MHz", "2"], - ["offset", "0", "deg", "10", "deg", "2"], - ] - } - ) +def test_build_axes_dict_from_cli_args_rejects_multiple_direction_coordinate_systems(): + with pytest.raises(ValueError, match="Cannot define both azimuth/zenith and ra/dec axes"): + build_axes_dict_from_cli_args( + { + "axis": [ + ["azimuth", "310", "deg", "20", "deg", "3"], + ["zenith", "30", "deg", "40", "deg", "2"], + ["ra", "0", "deg", "360", "deg", "36"], + ["dec", "-90", "deg", "90", "deg", "18"], + ["nsb", "4", "MHz", "5", "MHz", "2"], + ["offset", "0", "deg", "10", "deg", "2"], + ] + } + ) - assert "ra" in axes - assert "dec" in axes - assert "azimuth" not in axes - assert "zenith_angle" not in axes - assert axes["ra"]["local_azimuth_range"] == pytest.approx([310.0, 20.0]) - assert axes["ra"]["local_zenith_range"] == pytest.approx([30.0, 40.0]) + +def test_build_axes_dict_from_cli_args_rejects_invalid_density_unit(): + with pytest.raises(ValueError, match="direction_grid_density must be a float or quantity"): + build_axes_dict_from_cli_args( + { + "direction_grid_density": "1.0 1/s", + "axis": [ + ["ra", "0", "deg", "360", "deg", "36"], + ["dec", "-90", "deg", "90", "deg", "18"], + ["nsb", "4", "MHz", "5", "MHz", "2"], + ["offset", "0", "deg", "10", "deg", "2"], + ], + } + ) -def test_build_axes_dict_from_cli_args_rejects_duplicate_local_constraints(): - with pytest.raises(ValueError, match="Cannot define both 'local_zenith_range'"): +def test_build_axes_dict_from_cli_args_rejects_invalid_density_type(): + with pytest.raises(TypeError, match="direction_grid_density must be a number"): build_axes_dict_from_cli_args( { - "local_zenith_range": ["0", "deg", "70", "deg"], + "direction_grid_density": {"value": 1.0}, "axis": [ ["ra", "0", "deg", "360", "deg", "36"], ["dec", "-90", "deg", "90", "deg", "18"], - ["zenith", "0", "deg", "70", "deg", "2"], ["nsb", "4", "MHz", "5", "MHz", "2"], ["offset", "0", "deg", "10", "deg", "2"], ], From fbc370b4299071103b487b4c6bb25a2e50ee8939 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Thu, 28 May 2026 12:12:53 +0200 Subject: [PATCH 05/10] simplification --- .../applications/production_generate_grid.py | 2 +- .../production_configuration/angle_ranges.py | 42 +++++++++++++++++++ .../observation_grid.py | 26 +++--------- .../simulation_jobs.py | 25 +++-------- .../test_angle_ranges.py | 27 ++++++++++++ 5 files changed, 82 insertions(+), 40 deletions(-) create mode 100644 src/simtools/production_configuration/angle_ranges.py create mode 100644 tests/unit_tests/production_configuration/test_angle_ranges.py diff --git a/src/simtools/applications/production_generate_grid.py b/src/simtools/applications/production_generate_grid.py index e75331d735..eae91bf5e2 100644 --- a/src/simtools/applications/production_generate_grid.py +++ b/src/simtools/applications/production_generate_grid.py @@ -23,7 +23,7 @@ Example: ``--axis azimuth 310 deg 20 deg 3 linear``. Options for scaling are: linear, log, 1/cos. direction_grid_density (float or quantity, optional) - Direction-grid density in ``1/deg^2``. + Direction-grid density (typically in ``1/deg^2``). If provided, direction-axis binning (azimuth/zenith or ra/dec) is derived from range and density, while min/max values are kept from ``--axis`` definitions. In ``ra_dec`` mode, local-sky constraints can be defined with diff --git a/src/simtools/production_configuration/angle_ranges.py b/src/simtools/production_configuration/angle_ranges.py new file mode 100644 index 0000000000..92d74d1297 --- /dev/null +++ b/src/simtools/production_configuration/angle_ranges.py @@ -0,0 +1,42 @@ +"""Helpers for directed circular angle ranges used in production grids.""" + +import numpy as np + + +def directed_circular_span_degrees(axis_range): + """Return directed circular span (degrees) from start to end. + + Parameters + ---------- + axis_range : tuple or list + Two-element range ``(start, end)`` in degrees. + + Returns + ------- + float + Directed circular span in degrees, where full-circle ranges are represented as ``360.0``. + """ + start, end = axis_range + raw_span = abs(float(end) - float(start)) + if raw_span > 0.0 and np.isclose(raw_span % 360.0, 0.0): + return 360.0 + return float((end - start) % 360.0) + + +def ceil_with_tolerance(value): + """Ceil a float while avoiding near-integer floating-point artifacts. + + Parameters + ---------- + value : float + Input value. + + Returns + ------- + int + Ceiled integer with near-integer tolerance handling. + """ + nearest_integer = round(value) + if np.isclose(value, nearest_integer): + return int(nearest_integer) + return int(np.ceil(value)) diff --git a/src/simtools/production_configuration/observation_grid.py b/src/simtools/production_configuration/observation_grid.py index cd08395aa9..d8bc83d19c 100644 --- a/src/simtools/production_configuration/observation_grid.py +++ b/src/simtools/production_configuration/observation_grid.py @@ -13,6 +13,10 @@ from astropy.coordinates import AltAz, EarthLocation, SkyCoord from astropy.units import Quantity +from simtools.production_configuration.angle_ranges import ( + ceil_with_tolerance, + directed_circular_span_degrees, +) from simtools.production_configuration.corsika_limits_lookup import ( CorsikaLimitsLookup, attach_lookup_limits_to_point, @@ -321,22 +325,6 @@ def create_circular_binning(self, azimuth_range, num_bins): % 360.0 ) - @staticmethod - def _ceil_with_tolerance(value): - """Ceil a float while avoiding near-integer floating-point artifacts.""" - nearest_integer = round(value) - if np.isclose(value, nearest_integer): - return int(nearest_integer) - return int(np.ceil(value)) - - @staticmethod - def _directed_circular_span_degrees(azimuth_range): - """Return directed circular span (degrees) from start to end.""" - raw_span = abs(float(azimuth_range[1]) - float(azimuth_range[0])) - if raw_span > 0.0 and np.isclose(raw_span % 360.0, 0.0): - return 360.0 - return float((azimuth_range[1] - azimuth_range[0]) % 360.0) - @staticmethod def _is_in_directed_azimuth_range(azimuth_values_deg, azimuth_range): """Return mask of azimuth values inside directed circular range [start -> end].""" @@ -502,7 +490,7 @@ def _generate_adaptive_horizontal_grid(self): """Generate horizontal grid with azimuth binning adapted per zenith row.""" density = float(self.axes["azimuth"]["direction_grid_density"]) azimuth_range = self.axes["azimuth"]["range"] - azimuth_span = self._directed_circular_span_degrees(azimuth_range) + azimuth_span = directed_circular_span_degrees(azimuth_range) zenith_values = self.target_values["zenith_angle"] zenith_step = 1.0 / np.sqrt(density) @@ -528,9 +516,7 @@ def _generate_adaptive_horizontal_grid(self): if azimuth_span > 0.0 and altitude_cosine > 0.0: azimuth_bins = max( 1, - self._ceil_with_tolerance( - azimuth_span * density * zenith_step * altitude_cosine - ), + ceil_with_tolerance(azimuth_span * density * zenith_step * altitude_cosine), ) azimuth_values = self.create_circular_binning(azimuth_range, azimuth_bins) * u.deg diff --git a/src/simtools/production_configuration/simulation_jobs.py b/src/simtools/production_configuration/simulation_jobs.py index 97c774b48a..106ca6c3f9 100644 --- a/src/simtools/production_configuration/simulation_jobs.py +++ b/src/simtools/production_configuration/simulation_jobs.py @@ -17,6 +17,10 @@ from simtools.configuration.commandline_parser import CommandLineParser from simtools.layout.array_layout_utils import resolve_array_layout_name from simtools.model.site_model import SiteModel +from simtools.production_configuration.angle_ranges import ( + ceil_with_tolerance, + directed_circular_span_degrees, +) from simtools.production_configuration.corsika_limits_lookup import ( CorsikaLimitsLookup, attach_lookup_limits_to_point, @@ -212,23 +216,6 @@ def _parse_direction_grid_density(density_value): raise TypeError("direction_grid_density must be a number, string, or CLI-style token list.") -def _circular_span_degrees(axis_range): - """Return directed circular span (degrees) from start to end.""" - start, end = axis_range - raw_span = abs(float(end) - float(start)) - if raw_span > 0.0 and np.isclose(raw_span % 360.0, 0.0): - return 360.0 - return float((end - start) % 360.0) - - -def _ceil_with_tolerance(value): - """Ceil a float while avoiding near-integer floating-point artifacts.""" - nearest_integer = round(value) - if np.isclose(value, nearest_integer): - return int(nearest_integer) - return int(np.ceil(value)) - - def _mean_cosine_over_dec_span(dec_range): """Return mean cos(dec) over a declination interval (degrees).""" dec_min, dec_max = sorted((float(dec_range[0]), float(dec_range[1]))) @@ -264,7 +251,7 @@ def _apply_direction_grid_density(axis_configs, direction_axes, density): for axis_name in direction_axes: axis_range = axis_configs[axis_name]["range"] if axis_name == "azimuth": - span_degrees = _circular_span_degrees(axis_range) + span_degrees = directed_circular_span_degrees(axis_range) else: span_degrees = abs(axis_range[1] - axis_range[0]) @@ -273,7 +260,7 @@ def _apply_direction_grid_density(axis_configs, direction_axes, density): axis_configs[axis_name]["binning"] = max( 1, - _ceil_with_tolerance(span_degrees * density_sqrt), + ceil_with_tolerance(span_degrees * density_sqrt), ) diff --git a/tests/unit_tests/production_configuration/test_angle_ranges.py b/tests/unit_tests/production_configuration/test_angle_ranges.py new file mode 100644 index 0000000000..313c9aa009 --- /dev/null +++ b/tests/unit_tests/production_configuration/test_angle_ranges.py @@ -0,0 +1,27 @@ +import pytest + +from simtools.production_configuration.angle_ranges import ( + ceil_with_tolerance, + directed_circular_span_degrees, +) + + +def test_directed_circular_span_degrees_returns_directed_span_for_non_wrapping_range(): + assert directed_circular_span_degrees((10, 110)) == pytest.approx(100.0) + + +def test_directed_circular_span_degrees_returns_directed_span_for_wrapping_range(): + assert directed_circular_span_degrees((350, 10)) == pytest.approx(20.0) + + +def test_directed_circular_span_degrees_treats_full_circle_as_360_degrees(): + assert directed_circular_span_degrees((0, 360)) == pytest.approx(360.0) + assert directed_circular_span_degrees((15, 375)) == pytest.approx(360.0) + + +def test_ceil_with_tolerance_rounds_near_integer_without_overshooting(): + assert ceil_with_tolerance(10.00000000001) == 10 + + +def test_ceil_with_tolerance_uses_ceiling_for_non_integer_values(): + assert ceil_with_tolerance(10.2) == 11 From cb81c90b3c68acc50e92f73a0a87beeb65db568f Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Thu, 28 May 2026 12:16:48 +0200 Subject: [PATCH 06/10] fix integration test for plotting --- tests/integration_tests/config/plot_production_grid.yml | 4 ---- 1 file changed, 4 deletions(-) diff --git a/tests/integration_tests/config/plot_production_grid.yml b/tests/integration_tests/config/plot_production_grid.yml index a9b0096507..b15b03eb95 100644 --- a/tests/integration_tests/config/plot_production_grid.yml +++ b/tests/integration_tests/config/plot_production_grid.yml @@ -4,10 +4,6 @@ applications: configuration: grid_points_file: tests/resources/production_grid_points_horizontal.ecsv output_path: simtools-output - site: North - model_version: 6.0.2 - observation_time: '2017-09-16 00:00:00' - plot_ra_dec_tracks: true integration_tests: - output_file: production_grid_sky_projection.png test_name: default From b8a8f88ba5b60154eed8e1503b74775e15526b94 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Thu, 28 May 2026 12:17:11 +0200 Subject: [PATCH 07/10] sonar --- .../visualization/plot_production_grid.py | 104 ++++++--- .../test_observation_grid.py | 214 +++++++----------- 2 files changed, 152 insertions(+), 166 deletions(-) diff --git a/src/simtools/visualization/plot_production_grid.py b/src/simtools/visualization/plot_production_grid.py index 1ba1592249..ceae842cbb 100644 --- a/src/simtools/visualization/plot_production_grid.py +++ b/src/simtools/visualization/plot_production_grid.py @@ -249,53 +249,58 @@ def _split_visible_segments(mask): segments = np.split(visible_indices, split_indices) return [segment for segment in segments if len(segment) >= 2] - def plot_sky_projection(self, plot_ra_dec_tracks=False, dec_values=None): - """ - Create sky projection plots with Alt/Az and RA/Dec grid points. - - Parameters - ---------- - plot_ra_dec_tracks : bool - Whether to plot RA/Dec coordinate tracks. - dec_values : list of float, optional - List of declination values to plot as tracks. - """ - plot_points = self.normalize_grid_points() - has_radec_points = any( - point["ra"] is not None and point["dec"] is not None for point in plot_points - ) - show_radec_panel = self.has_radec_columns and has_radec_points + @staticmethod + def _has_plottable_radec_points(plot_points): + """Return whether normalized points include plottable RA/Dec coordinates.""" + return any(point["ra"] is not None and point["dec"] is not None for point in plot_points) + @staticmethod + def _create_projection_axes(show_radec_panel): + """Create figure and projection axes based on available coordinate panels.""" figure = plt.figure(figsize=(15, 7) if show_radec_panel else (8, 7)) if show_radec_panel: altaz_axis = figure.add_subplot(1, 2, 1, projection="polar") radec_axis = figure.add_subplot(1, 2, 2) - else: - altaz_axis = figure.add_subplot(1, 1, 1, projection="polar") - radec_axis = None + return figure, altaz_axis, radec_axis - altaz_count = self._plot_altaz_points(altaz_axis, plot_points) - radec_count = self._plot_radec_points(radec_axis, plot_points) if radec_axis else 0 - if plot_ra_dec_tracks: - if dec_values: - logger.info( - "RA/Dec tracks are disabled in file-driven plotting mode " - "(ignoring manual dec_values)" - ) - else: - logger.info("RA/Dec tracks are disabled in file-driven plotting mode") + altaz_axis = figure.add_subplot(1, 1, 1, projection="polar") + return figure, altaz_axis, None + @staticmethod + def _add_axis_legend_if_present(axis, location_kwargs): + """Add a legend to an axis only when visible labels are present.""" + _, labels = axis.get_legend_handles_labels() + if any(label and not label.startswith("_") for label in labels): + axis.legend(**location_kwargs) + + def _add_panel_legends(self, altaz_axis, radec_axis, altaz_count, radec_count): + """Add legends to Alt/Az and RA/Dec panels when needed.""" if altaz_count > 0: - _, altaz_labels = altaz_axis.get_legend_handles_labels() - if any(label and not label.startswith("_") for label in altaz_labels): - altaz_axis.legend(loc="upper left", bbox_to_anchor=(1.0, 1.15)) + self._add_axis_legend_if_present( + altaz_axis, + {"loc": "upper left", "bbox_to_anchor": (1.0, 1.15)}, + ) if radec_axis and radec_count > 0: - _, radec_labels = radec_axis.get_legend_handles_labels() - if any(label and not label.startswith("_") for label in radec_labels): - radec_axis.legend(loc="upper right") + self._add_axis_legend_if_present(radec_axis, {"loc": "upper right"}) - figure.suptitle("Production Grid Points", fontsize=14, y=0.98) + @staticmethod + def _log_track_request_status(plot_ra_dec_tracks, dec_values): + """Log status for manual RA/Dec track requests in file-driven plotting mode.""" + if not plot_ra_dec_tracks: + return + + if dec_values: + logger.info( + "RA/Dec tracks are disabled in file-driven plotting mode " + "(ignoring manual dec_values)" + ) + return + + logger.info("RA/Dec tracks are disabled in file-driven plotting mode") + + def _build_subtitle_lines(self): + """Build subtitle lines from available grid metadata.""" subtitle_lines = [] if self.grid_metadata.get("site"): subtitle_lines.append(f"Site: {self.grid_metadata['site']}") @@ -308,7 +313,11 @@ def plot_sky_projection(self, plot_ra_dec_tracks=False, dec_values=None): subtitle_lines.append( f"Observation time (UTC): {self.grid_metadata['time_of_observation_utc']}" ) + return subtitle_lines + @staticmethod + def _render_subtitle_lines(figure, subtitle_lines): + """Render subtitle lines below the figure title.""" for index, subtitle_line in enumerate(subtitle_lines): figure.text( 0.5, @@ -319,6 +328,29 @@ def plot_sky_projection(self, plot_ra_dec_tracks=False, dec_values=None): fontsize=10, ) + def plot_sky_projection(self, plot_ra_dec_tracks=False, dec_values=None): + """ + Create sky projection plots with Alt/Az and RA/Dec grid points. + + Parameters + ---------- + plot_ra_dec_tracks : bool + Whether to plot RA/Dec coordinate tracks. + dec_values : list of float, optional + List of declination values to plot as tracks. + """ + plot_points = self.normalize_grid_points() + show_radec_panel = self.has_radec_columns and self._has_plottable_radec_points(plot_points) + figure, altaz_axis, radec_axis = self._create_projection_axes(show_radec_panel) + + altaz_count = self._plot_altaz_points(altaz_axis, plot_points) + radec_count = self._plot_radec_points(radec_axis, plot_points) if radec_axis else 0 + self._log_track_request_status(plot_ra_dec_tracks, dec_values) + self._add_panel_legends(altaz_axis, radec_axis, altaz_count, radec_count) + + figure.suptitle("Production Grid Points", fontsize=14, y=0.98) + self._render_subtitle_lines(figure, self._build_subtitle_lines()) + figure.tight_layout(rect=(0.0, 0.0, 1.0, 0.92)) output_file = self.output_path / ( diff --git a/tests/unit_tests/production_configuration/test_observation_grid.py b/tests/unit_tests/production_configuration/test_observation_grid.py index 93b323d117..ea2561b7ad 100644 --- a/tests/unit_tests/production_configuration/test_observation_grid.py +++ b/tests/unit_tests/production_configuration/test_observation_grid.py @@ -9,18 +9,57 @@ from simtools.production_configuration.observation_grid import ProductionGridEngine +DEFAULT_OBSERVING_LOCATION = EarthLocation( + lat=28.76 * u.deg, + lon=-17.89 * u.deg, + height=2200 * u.m, +) +DEFAULT_OBSERVATION_TIME = Time("2020-01-01 00:00:00", scale="utc") -def test_generate_simulation_grid_keeps_horizontal_coordinates_for_radec_axes(): + +def _make_engine(axes=None, **kwargs): + """Build a production-grid engine with default wrapped axes input.""" + return ProductionGridEngine(axes={"axes": axes or {}}, **kwargs) + + +def _adaptive_radec_axes(local_zenith_range=None, local_azimuth_range=None): + """Return baseline RA/Dec adaptive-density axes with optional local constraints.""" axes = { - "axes": { + "ra": { + "range": [0, 360], + "binning": 36, + "units": "deg", + "direction_grid_density": 0.25, + }, + "dec": {"range": [-40, 80], "binning": 10, "units": "deg"}, + "nsb_level": {"range": [4, 4], "binning": 1, "units": "MHz"}, + "offset": {"range": [0, 10], "binning": 2, "units": "deg"}, + } + if local_zenith_range is not None: + axes["ra"]["local_zenith_range"] = local_zenith_range + if local_azimuth_range is not None: + axes["ra"]["local_azimuth_range"] = local_azimuth_range + return axes + + +def _single_bin_horizontal_axes(azimuth_range=None): + """Return single-bin horizontal axes used for interpolation-limit tests.""" + az_range = [0, 0] if azimuth_range is None else azimuth_range + return { + "zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}, + "azimuth": {"range": az_range, "binning": 1 if az_range == [0, 0] else 3, "units": "deg"}, + "nsb_level": {"range": [1, 1], "binning": 1, "units": "1"}, + } + + +def test_generate_simulation_grid_keeps_horizontal_coordinates_for_radec_axes(): + engine = _make_engine( + axes={ "ra": {"range": [0, 0], "binning": 1, "scaling": "linear", "units": "deg"}, "dec": {"range": [0, 0], "binning": 1, "scaling": "linear", "units": "deg"}, - } - } - engine = ProductionGridEngine( - axes=axes, + }, coordinate_system="ra_dec", - observing_location=EarthLocation(lat=28.76 * u.deg, lon=-17.89 * u.deg, height=2200 * u.m), + observing_location=DEFAULT_OBSERVING_LOCATION, time_of_observation=Time("2017-09-16 00:00:00", scale="utc"), lookup_table=None, ) @@ -34,22 +73,20 @@ def test_generate_simulation_grid_keeps_horizontal_coordinates_for_radec_axes(): def test_require_time_of_observation_raises_without_time(): - engine = ProductionGridEngine(axes={"axes": {}}, time_of_observation=None) + engine = _make_engine(time_of_observation=None) with pytest.raises(ValueError, match="Observing time"): engine._require_time_of_observation() def test_get_max_zenith_for_radec_mode_reads_axis_range(): - engine = ProductionGridEngine( - axes={"axes": {"zenith_angle": {"range": [10, 60], "binning": 2, "units": "deg"}}} - ) + engine = _make_engine(axes={"zenith_angle": {"range": [10, 60], "binning": 2, "units": "deg"}}) assert engine._get_max_zenith_for_radec_mode() == 60 def test_get_max_zenith_for_radec_mode_raises_for_invalid_axis_definition(): - engine = ProductionGridEngine(axes={"axes": {}}) + engine = _make_engine() with pytest.raises(ValueError, match="valid two-element 'range'"): engine._get_max_zenith_for_radec_mode() @@ -75,7 +112,7 @@ def test_generate_target_values_supports_log_and_inverse_cos_scaling(): def test_add_lookup_limits_to_point_uses_quantity_nsb_level(): - engine = ProductionGridEngine(axes={"axes": {}}, lookup_table=None) + engine = _make_engine(lookup_table=None) engine.lookup_table = "limits.ecsv" engine._interpolate_limits_for_point = Mock( return_value={ @@ -97,7 +134,7 @@ def test_add_lookup_limits_to_point_uses_quantity_nsb_level(): def test_interpolate_limits_for_point_delegates_to_lookup_helper(): - engine = ProductionGridEngine(axes={"axes": {}}, lookup_table=None) + engine = _make_engine(lookup_table=None) engine._limits_lookup = Mock() engine._limits_lookup.interpolate_point.return_value = {"lower_energy_threshold": 0.02} @@ -108,7 +145,7 @@ def test_interpolate_limits_for_point_delegates_to_lookup_helper(): def test_prepare_lookup_table_limits_for_point_interpolation_prepares(): - engine = ProductionGridEngine(axes={"axes": {}}, lookup_table=None) + engine = _make_engine(lookup_table=None) engine._limits_lookup = Mock() engine._prepare_lookup_table_limits_for_point_interpolation() @@ -117,9 +154,7 @@ def test_prepare_lookup_table_limits_for_point_interpolation_prepares(): def test_generate_extra_axis_combinations_returns_empty_placeholder_for_no_extra_axes(): - engine = ProductionGridEngine( - axes={"axes": {"zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}}} - ) + engine = _make_engine(axes={"zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}}) keys, units, combinations = engine._generate_extra_axis_combinations(("zenith_angle",)) @@ -147,7 +182,7 @@ def test_generate_extra_axis_combinations_returns_mesh_for_remaining_axes(): def test_create_circular_binning_uses_directed_span_from_start_to_end(): - engine = ProductionGridEngine(axes={"axes": {}}) + engine = _make_engine() binning = engine.create_circular_binning((10, 350), 3) @@ -155,7 +190,7 @@ def test_create_circular_binning_uses_directed_span_from_start_to_end(): def test_create_circular_binning_uses_clockwise_path_when_shorter(): - engine = ProductionGridEngine(axes={"axes": {}}) + engine = _make_engine() binning = engine.create_circular_binning((350, 10), 3) @@ -163,7 +198,7 @@ def test_create_circular_binning_uses_clockwise_path_when_shorter(): def test_create_circular_binning_covers_requested_range_0_to_240(): - engine = ProductionGridEngine(axes={"axes": {}}) + engine = _make_engine() binning = engine.create_circular_binning((0, 240), 3) @@ -232,19 +267,9 @@ def test_generate_radec_grid_uses_adaptive_ra_bins_per_dec_strip(): def test_generate_grid_from_radec_axes_routes_to_adaptive_density_path(): engine = ProductionGridEngine( - axes={ - "ra": { - "range": [0, 360], - "binning": 36, - "units": "deg", - "direction_grid_density": 0.25, - }, - "dec": {"range": [-40, 80], "binning": 10, "units": "deg"}, - "nsb_level": {"range": [4, 4], "binning": 1, "units": "MHz"}, - "offset": {"range": [0, 10], "binning": 2, "units": "deg"}, - }, + axes=_adaptive_radec_axes(), coordinate_system="ra_dec", - time_of_observation=Time("2020-01-01 00:00:00", scale="utc"), + time_of_observation=DEFAULT_OBSERVATION_TIME, ) expected_grid = [{"ra": 1 * u.deg, "dec": 2 * u.deg}] engine._generate_adaptive_radec_grid = Mock(return_value=expected_grid) @@ -257,24 +282,10 @@ def test_generate_grid_from_radec_axes_routes_to_adaptive_density_path(): def test_generate_radec_grid_adaptive_density_keeps_only_visible_nodes(): engine = ProductionGridEngine( - axes={ - "ra": { - "range": [0, 360], - "binning": 36, - "units": "deg", - "direction_grid_density": 0.25, - }, - "dec": {"range": [-40, 80], "binning": 10, "units": "deg"}, - "nsb_level": {"range": [4, 4], "binning": 1, "units": "MHz"}, - "offset": {"range": [0, 10], "binning": 2, "units": "deg"}, - }, + axes=_adaptive_radec_axes(), coordinate_system="ra_dec", - observing_location=EarthLocation( - lat=28.76 * u.deg, - lon=-17.89 * u.deg, - height=2200 * u.m, - ), - time_of_observation=Time("2020-01-01 00:00:00", scale="utc"), + observing_location=DEFAULT_OBSERVING_LOCATION, + time_of_observation=DEFAULT_OBSERVATION_TIME, ) grid = engine._generate_adaptive_radec_grid(include_horizontal_coordinates=True) @@ -285,25 +296,10 @@ def test_generate_radec_grid_adaptive_density_keeps_only_visible_nodes(): def test_generate_radec_grid_adaptive_density_applies_zenith_constraint(): engine = ProductionGridEngine( - axes={ - "ra": { - "range": [0, 360], - "binning": 36, - "units": "deg", - "direction_grid_density": 0.25, - "local_zenith_range": [0, 70], - }, - "dec": {"range": [-40, 80], "binning": 10, "units": "deg"}, - "nsb_level": {"range": [4, 4], "binning": 1, "units": "MHz"}, - "offset": {"range": [0, 10], "binning": 2, "units": "deg"}, - }, + axes=_adaptive_radec_axes(local_zenith_range=[0, 70]), coordinate_system="ra_dec", - observing_location=EarthLocation( - lat=28.76 * u.deg, - lon=-17.89 * u.deg, - height=2200 * u.m, - ), - time_of_observation=Time("2020-01-01 00:00:00", scale="utc"), + observing_location=DEFAULT_OBSERVING_LOCATION, + time_of_observation=DEFAULT_OBSERVATION_TIME, ) grid = engine._generate_adaptive_radec_grid(include_horizontal_coordinates=True) @@ -314,26 +310,13 @@ def test_generate_radec_grid_adaptive_density_applies_zenith_constraint(): def test_generate_radec_grid_adaptive_density_applies_azimuth_constraint(): engine = ProductionGridEngine( - axes={ - "ra": { - "range": [0, 360], - "binning": 36, - "units": "deg", - "direction_grid_density": 0.25, - "local_zenith_range": [0, 70], - "local_azimuth_range": [300, 60], - }, - "dec": {"range": [-40, 80], "binning": 10, "units": "deg"}, - "nsb_level": {"range": [4, 4], "binning": 1, "units": "MHz"}, - "offset": {"range": [0, 10], "binning": 2, "units": "deg"}, - }, - coordinate_system="ra_dec", - observing_location=EarthLocation( - lat=28.76 * u.deg, - lon=-17.89 * u.deg, - height=2200 * u.m, + axes=_adaptive_radec_axes( + local_zenith_range=[0, 70], + local_azimuth_range=[300, 60], ), - time_of_observation=Time("2020-01-01 00:00:00", scale="utc"), + coordinate_system="ra_dec", + observing_location=DEFAULT_OBSERVING_LOCATION, + time_of_observation=DEFAULT_OBSERVATION_TIME, ) grid = engine._generate_adaptive_radec_grid(include_horizontal_coordinates=True) @@ -355,7 +338,7 @@ def test_is_in_directed_azimuth_range_returns_all_true_for_full_circle(): def test_create_circular_binning_treats_full_circle_range_as_full_span(): - engine = ProductionGridEngine(axes={"axes": {}}) + engine = _make_engine() binning = engine.create_circular_binning((0, 360), 4) @@ -363,16 +346,15 @@ def test_create_circular_binning_treats_full_circle_range_as_full_span(): def test_convert_altaz_to_radec_raises_without_time_of_observation(): - engine = ProductionGridEngine(axes={"axes": {}}, time_of_observation=None) + engine = _make_engine(time_of_observation=None) with pytest.raises(ValueError, match="time_of_observation"): engine.convert_altaz_to_radec(45 * u.deg, 180 * u.deg) def test_convert_altaz_to_radec_returns_icrs_coordinates(): - engine = ProductionGridEngine( - axes={"axes": {}}, - observing_location=EarthLocation(lat=28.76 * u.deg, lon=-17.89 * u.deg, height=2200 * u.m), + engine = _make_engine( + observing_location=DEFAULT_OBSERVING_LOCATION, time_of_observation=Time("2017-09-16 00:00:00", scale="utc"), ) @@ -383,7 +365,7 @@ def test_convert_altaz_to_radec_returns_icrs_coordinates(): def test_convert_coordinates_drops_horizontal_coordinates_when_requested(): - engine = ProductionGridEngine(axes={"axes": {}}, coordinate_system="ra_dec") + engine = _make_engine(coordinate_system="ra_dec") engine.convert_altaz_to_radec = Mock(return_value=Mock(ra=Mock(deg=10), dec=Mock(deg=-20))) points = [{"zenith_angle": 20 * u.deg, "azimuth": 180 * u.deg}] @@ -396,7 +378,7 @@ def test_convert_coordinates_drops_horizontal_coordinates_when_requested(): def test_convert_coordinates_keeps_points_without_horizontal_values(): - engine = ProductionGridEngine(axes={"axes": {}}, coordinate_system="ra_dec") + engine = _make_engine(coordinate_system="ra_dec") points = [{"ra": 10 * u.deg, "dec": -20 * u.deg}] converted = engine.convert_coordinates(points, keep_horizontal_coordinates=False) @@ -432,8 +414,8 @@ def test_init_with_horizontal_lookup_applies_interpolated_limits(mock_corsika_li "lower_energy_threshold": np.array([[[0.02]]]) } - engine = ProductionGridEngine( - axes={"axes": {"zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}}}, + engine = _make_engine( + axes={"zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}}, lookup_table="limits.ecsv", array_layout_name="alpha", ) @@ -442,15 +424,7 @@ def test_init_with_horizontal_lookup_applies_interpolated_limits(mock_corsika_li def test_generate_horizontal_grid_uses_interpolated_limit_arrays(): - engine = ProductionGridEngine( - axes={ - "axes": { - "zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}, - "azimuth": {"range": [0, 0], "binning": 1, "units": "deg"}, - "nsb_level": {"range": [1, 1], "binning": 1, "units": "1"}, - } - } - ) + engine = _make_engine(axes=_single_bin_horizontal_axes()) engine.interpolated_limits = { "lower_energy_threshold": np.array([[[0.02]]]), "upper_scatter_radius": np.array([[[120.0]]]), @@ -491,15 +465,7 @@ def test_generate_radec_grid_direction_points_filters_by_max_zenith( def test_generate_horizontal_grid_handles_partial_interpolated_limit_arrays(): - engine = ProductionGridEngine( - axes={ - "axes": { - "zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}, - "azimuth": {"range": [0, 0], "binning": 1, "units": "deg"}, - "nsb_level": {"range": [1, 1], "binning": 1, "units": "1"}, - } - } - ) + engine = _make_engine(axes=_single_bin_horizontal_axes()) engine.interpolated_limits = { "upper_scatter_radius": np.array([[[120.0]]]), "viewcone_radius": np.array([[[3.0]]]), @@ -513,15 +479,7 @@ def test_generate_horizontal_grid_handles_partial_interpolated_limit_arrays(): def test_generate_horizontal_grid_with_circular_azimuth_binning_uses_correct_indices(): - engine = ProductionGridEngine( - axes={ - "axes": { - "zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}, - "azimuth": {"range": [10, 350], "binning": 3, "units": "deg"}, - "nsb_level": {"range": [1, 1], "binning": 1, "units": "1"}, - } - } - ) + engine = _make_engine(axes=_single_bin_horizontal_axes(azimuth_range=[10, 350])) # Directed span from start to end gives [10, 180, 350] assert np.allclose(engine.target_values["azimuth"].value, np.array([10.0, 180.0, 350.0])) @@ -541,9 +499,7 @@ def test_generate_horizontal_grid_with_circular_azimuth_binning_uses_correct_ind def test_generate_grid_radec_mode_adds_extra_axis_quantities(): - engine = ProductionGridEngine( - axes={"axes": {"zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}}} - ) + engine = _make_engine(axes={"zenith_angle": {"range": [20, 20], "binning": 1, "units": "deg"}}) engine.coordinate_system = "ra_dec" engine._generate_radec_grid_direction_points = Mock( return_value=[{"zenith_angle": 20 * u.deg, "azimuth": 180 * u.deg}] @@ -562,9 +518,7 @@ def test_generate_grid_radec_mode_adds_extra_axis_quantities(): def test_generate_grid_radec_mode_uses_direction_points_when_ra_dec_axes_missing(): - engine = ProductionGridEngine( - axes={"axes": {"zenith_angle": {"range": [20, 20], "binning": 1}}} - ) + engine = _make_engine(axes={"zenith_angle": {"range": [20, 20], "binning": 1}}) engine.coordinate_system = "ra_dec" engine._generate_radec_grid_direction_points = Mock( return_value=[{"zenith_angle": 20 * u.deg, "azimuth": 180 * u.deg}] @@ -579,7 +533,7 @@ def test_generate_grid_radec_mode_uses_direction_points_when_ra_dec_axes_missing def test_generate_simulation_grid_uses_horizontal_grid_for_horizontal_mode(): - engine = ProductionGridEngine(axes={"axes": {}}) + engine = _make_engine() engine.coordinate_system = "horizontal" engine._generate_horizontal_grid = Mock(return_value=[{"azimuth": 0 * u.deg}]) From b368cdc83fa1eeebf4315728c74e19e311f1ae4a Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Thu, 28 May 2026 12:19:36 +0200 Subject: [PATCH 08/10] changelog --- docs/changes/2213.feature.md | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/changes/2213.feature.md diff --git a/docs/changes/2213.feature.md b/docs/changes/2213.feature.md new file mode 100644 index 0000000000..e4cbb84716 --- /dev/null +++ b/docs/changes/2213.feature.md @@ -0,0 +1 @@ +Add functionality to to define simulation production grids using a grid density parameter. Improve plotting. From b1863a22e0507d35580608b3bb5518fdf68b6af8 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Thu, 28 May 2026 12:22:13 +0200 Subject: [PATCH 09/10] docs --- docs/source/api-reference/production_configuration.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/docs/source/api-reference/production_configuration.md b/docs/source/api-reference/production_configuration.md index b8cf346fb2..db79af5130 100644 --- a/docs/source/api-reference/production_configuration.md +++ b/docs/source/api-reference/production_configuration.md @@ -53,6 +53,16 @@ the calculation of the number of events to be simulated given a pre-determined m ``` +(angle-ranges)= + +## angle_ranges + +```{eval-rst} +.. automodule:: production_configuration.angle_ranges + :members: +``` + + (interpolation-handler)= ## interpolation_handler From e9f45ce2714a322bb462bd47316c672a3bbbd2e3 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Thu, 28 May 2026 13:17:27 +0200 Subject: [PATCH 10/10] integration tests --- docs/changes/2213.feature.md | 2 +- src/simtools/visualization/plot_production_grid.py | 10 ++++++---- .../production_generate_grid_horizontal_density.yml | 2 +- .../test_observation_grid.py | 12 ++++++++++++ 4 files changed, 20 insertions(+), 6 deletions(-) diff --git a/docs/changes/2213.feature.md b/docs/changes/2213.feature.md index e4cbb84716..8de8cc12f0 100644 --- a/docs/changes/2213.feature.md +++ b/docs/changes/2213.feature.md @@ -1 +1 @@ -Add functionality to to define simulation production grids using a grid density parameter. Improve plotting. +Add functionality to define simulation production grids using a grid density parameter. Improve plotting. diff --git a/src/simtools/visualization/plot_production_grid.py b/src/simtools/visualization/plot_production_grid.py index ceae842cbb..e5d2058313 100644 --- a/src/simtools/visualization/plot_production_grid.py +++ b/src/simtools/visualization/plot_production_grid.py @@ -175,7 +175,7 @@ def _normalize_grid_point(self, point): "zenith": None, "ra": float(ra % 360.0), "dec": float(dec), - "visible_in_altaz": False, + "visible_in_altaz": None, } logger.warning(f"Skipping point without supported coordinates: {point}") @@ -335,9 +335,11 @@ def plot_sky_projection(self, plot_ra_dec_tracks=False, dec_values=None): Parameters ---------- plot_ra_dec_tracks : bool - Whether to plot RA/Dec coordinate tracks. + Kept for backward-compatible CLI/API usage. + In file-driven plotting mode, RA/Dec tracks are not rendered and this flag is ignored. dec_values : list of float, optional - List of declination values to plot as tracks. + Kept for backward-compatible CLI/API usage. + In file-driven plotting mode, this argument is ignored. """ plot_points = self.normalize_grid_points() show_radec_panel = self.has_radec_columns and self._has_plottable_radec_points(plot_points) @@ -499,7 +501,7 @@ def _plot_altaz_points(self, axis, plot_points): x_transform=np.radians, ) - hidden_points = sum(not point["visible_in_altaz"] for point in plot_points) + hidden_points = sum(point["visible_in_altaz"] is False for point in plot_points) if hidden_points > 0: logger.info(f"Skipping {hidden_points} RA/Dec points below the horizon in Alt/Az panel") return plotted_points diff --git a/tests/integration_tests/config/production_generate_grid_horizontal_density.yml b/tests/integration_tests/config/production_generate_grid_horizontal_density.yml index 2d9608b60e..a8def9e973 100644 --- a/tests/integration_tests/config/production_generate_grid_horizontal_density.yml +++ b/tests/integration_tests/config/production_generate_grid_horizontal_density.yml @@ -28,7 +28,7 @@ applications: site: North view_cone: 0 deg 10 deg integration_tests: - - output_file: job_grid_horizontal_density.ecsv + - output_file: job_grid_horizontal_grid_density.ecsv model_version_use_current: true test_name: production_generate_grid_horizontal schema_name: application_workflow.metaschema diff --git a/tests/unit_tests/production_configuration/test_observation_grid.py b/tests/unit_tests/production_configuration/test_observation_grid.py index ea2561b7ad..e9cc5ae9b1 100644 --- a/tests/unit_tests/production_configuration/test_observation_grid.py +++ b/tests/unit_tests/production_configuration/test_observation_grid.py @@ -6,6 +6,7 @@ from astropy.coordinates import EarthLocation, SkyCoord from astropy.tests.helper import assert_quantity_allclose from astropy.time import Time +from astropy.utils import iers from simtools.production_configuration.observation_grid import ProductionGridEngine @@ -17,6 +18,17 @@ DEFAULT_OBSERVATION_TIME = Time("2020-01-01 00:00:00", scale="utc") +@pytest.fixture(autouse=True, scope="module") +def disable_iers_auto_download(): + """Disable IERS auto-download to keep tests deterministic and offline-safe.""" + previous_auto_download = iers.conf.auto_download + iers.conf.auto_download = False + try: + yield + finally: + iers.conf.auto_download = previous_auto_download + + def _make_engine(axes=None, **kwargs): """Build a production-grid engine with default wrapped axes input.""" return ProductionGridEngine(axes={"axes": axes or {}}, **kwargs)