Allow to define production grids using a grid density parameter

docs/changes/2213.feature.md

@@ -0,0 +1 @@
+Add functionality to define simulation production grids using a grid density parameter. Improve plotting.

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

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(),
     )
 

src/simtools/applications/production_generate_grid.py

@@ -24,6 +24,18 @@
     ``--axis <name> <min> <unit> <max> <unit> <binning> [scaling]``.
     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 (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
+    ``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
@@ -72,6 +84,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 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" \
+            --corsika_limits tests/resources/corsika_simulation_limits/merged_corsika_limits.ecsv
 """
 
 from simtools.application_control import build_application
@@ -104,6 +133,37 @@ 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",
+        nargs="+",
+        required=False,
+        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, use "
+            "local_zenith_range/local_azimuth_range to filter generated points in local sky."
+        ),
+    )
+    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(
         "--output_file",
         type=str,

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))

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)