"""
Create QGIS layers showing grounding and allision results with probability attributes.

Each layer contains lines representing obstacles that were hit, with attributes:
- total_probability: Sum of all leg contributions
- leg_X_prob: Contribution from each leg (aggregated across all directions)

Layers are styled with graduated symbology:
- Red = highest probability
- Yellow = medium probability
- Green = lowest probability

Per-segment probability calculation:
Each segment's probability is calculated based on its orientation relative to drift directions.
A segment can only be "hit" by drift from directions that face the segment's exposed side.
For example, a north-facing edge (normal pointing north) can only be hit by southward drift.
"""
from typing import Any
import numpy as np
from qgis.PyQt.QtCore import QVariant
from qgis.core import (
    QgsVectorLayer,
    QgsFeature,
    QgsGeometry,
    QgsField,
    QgsProject,
    QgsGraduatedSymbolRenderer,
    QgsRendererRange,
    QgsLineSymbol,
    QgsMarkerSymbol,
    QgsPointXY,
)
from qgis.PyQt.QtGui import QColor
import logging

logger = logging.getLogger(__name__)

# Human-readable aliases for the result-layer fields. Keys are the QGIS field
# names (snake_case, used in code and attribute lookups); values are the
# aliases shown in the attribute table and identify-feature popup.
_FIELD_ALIASES = {
    'total_edge_probability': 'Total edge probability',
    'object_probability': 'Object probability',
    'obstacle_id': 'Obstacle id',
    'segment_idx': 'Segment idx',
    'segment_id': 'Edge fid',
    'normal_deg': 'Normal (deg)',
    'value': 'Value',
}


def _set_result_field_aliases(layer: QgsVectorLayer) -> None:
    """Apply user-friendly aliases to known result-layer field names.

    Field aliases only affect display in the attribute table / identify
    panel; the underlying snake_case field names stay stable so code,
    expressions, and saved QGIS projects keep working.
    """
    fields = layer.fields()
    for name, alias in _FIELD_ALIASES.items():
        idx = fields.indexOf(name)
        if idx >= 0:
            layer.setFieldAlias(idx, alias)


def _enable_fid_labels(layer: QgsVectorLayer, field_name: str) -> None:
    """Turn on text labels for ``layer`` showing each feature's ``field_name``.

    Used on the result layers so the user can see at a glance which
    obstacle / leg fid every feature corresponds to without opening the
    attribute table.
    """
    if layer is None:
        return
    if layer.fields().indexOf(field_name) < 0:
        return
    try:
        from qgis.core import (
            QgsPalLayerSettings,
            QgsTextFormat,
            QgsVectorLayerSimpleLabeling,
        )
        from qgis.PyQt.QtGui import QColor as _QColor
    except Exception:  # nosec B110 B112
        return

    settings = QgsPalLayerSettings()
    settings.fieldName = field_name
    settings.isExpression = False
    settings.enabled = True

    # A small but readable text format with a white halo so the label
    # stays legible regardless of the underlying symbology colour.
    fmt = QgsTextFormat()
    fmt.setSize(8)
    fmt.setColor(_QColor('black'))
    try:
        buffer = fmt.buffer()
        buffer.setEnabled(True)
        buffer.setSize(0.6)
        buffer.setColor(_QColor('white'))
        fmt.setBuffer(buffer)
    except Exception:  # nosec B110 B112
        pass
    settings.setFormat(fmt)

    # Place lines along the line; centre points / polygons.
    geom_type = layer.geometryType()
    try:
        from qgis.core import Qgis
        if geom_type == Qgis.GeometryType.Line:
            settings.placement = QgsPalLayerSettings.Placement.Line
        else:
            settings.placement = QgsPalLayerSettings.Placement.OverPoint
    except Exception:  # nosec B110 B112
        pass

    layer.setLabeling(QgsVectorLayerSimpleLabeling(settings))
    layer.setLabelsEnabled(True)
    layer.triggerRepaint()


# Drift directions used in the model (standard nautical convention: 0=N, 90=E, CW)
DRIFT_DIRECTIONS = {
    'North': 0,
    'NorthEast': 45,
    'East': 90,
    'SouthEast': 135,
    'South': 180,
    'SouthWest': 225,
    'West': 270,
    'NorthWest': 315,
}


def _segment_normal_angle(x1: float, y1: float, x2: float, y2: float) -> float:
    """
    Calculate the outward normal angle of a segment.

    For a polygon with counter-clockwise vertices (the standard for
    WGS84 GeoJSON/WKT), the outward normal points to the RIGHT of the
    segment direction when walking along the boundary.

    Args:
        x1, y1: Start point of segment.
        x2, y2: End point of segment.

    Returns:
        Standard nautical bearing in degrees (0°=N, 90°=E, 180°=S,
        270°=W; clockwise from North).  Matches the compass convention
        used by the rest of the drifting model, so the angle can be
        compared directly with ``drift_direction`` values from
        :func:`_parse_angle_from_key` (which are produced as
        ``d_idx * 45`` in the same convention).
    """
    dx = x2 - x1
    dy = y2 - y1

    # Right-perpendicular of (dx, dy) in UTM is (dy, -dx) -- the outward
    # normal for a CCW-wound polygon.
    nx = dy
    ny = -dx

    # Nautical bearing of (nx, ny) in UTM (+x=East, +y=North):
    # bearing = atan2(nx, ny) measured clockwise from +y.
    return float(np.degrees(np.arctan2(nx, ny)) % 360)


def _exposure_factor(segment_normal: float, drift_direction: float) -> float:
    """
    Calculate how exposed a segment is to drift from a given direction.

    A segment is fully exposed (factor=1.0) when drift comes directly toward
    the segment's face (i.e., drift direction is OPPOSITE to the segment normal).
    Exposure decreases with cos(angle) and is 0 when the drift is perpendicular
    or coming from behind the segment.

    Example:
    - A segment facing north (normal=0°) is hit by ships drifting SOUTH (toward 180°)
    - Ships drifting south are moving toward the north-facing edge

    Args:
        segment_normal: Outward normal angle of segment (OMRAT degrees: 0=N, 90=W, 180=S, 270=E)
        drift_direction: Direction ships drift TOWARD (OMRAT degrees)

    Returns:
        Exposure factor between 0.0 and 1.0
    """
    # The drift direction is where ships are drifting TO.
    # A segment facing north (normal=0) is hit by ships drifting SOUTH (toward 180).
    # The drift direction that hits a segment is OPPOSITE to the segment's normal.
    #
    # So we compare drift_direction with (segment_normal + 180) % 360

    # Direction that would hit this segment head-on
    hit_direction = (segment_normal + 180) % 360

    # Angular difference between actual drift and the hit direction
    angle_diff = abs(drift_direction - hit_direction)
    if angle_diff > 180:
        angle_diff = 360 - angle_diff

    # Exposure is cos(angle_diff), clamped to [0, 1]
    # Full exposure when angle_diff = 0 (drift comes directly at segment face)
    # Zero exposure when angle_diff >= 90 (drift perpendicular or from behind)
    if angle_diff >= 90:
        return 0.0

    return np.cos(np.radians(angle_diff))


def _parse_angle_from_key(leg_dir_key: str) -> float | None:
    """
    Parse drift angle from leg-direction key.

    Key format: "seg_id:direction_label:angle" (e.g., "1:NNW:45")
    The angle is the third component (d_idx * 45 from run_calculations.py).

    Returns:
        Drift angle in degrees (0, 45, 90, etc.) or None if parsing fails
    """
    parts = leg_dir_key.split(':')
    if len(parts) >= 3:
        try:
            return float(parts[2])
        except (ValueError, TypeError):
            pass
    return None


def _aggregate_by_leg(leg_contributions: dict[str, float]) -> dict[str, float]:
    """
    Aggregate leg-direction contributions to per-leg totals.

    The leg keys are in format "seg_id:direction:angle" (e.g., "1:North:0").
    We want to sum all directions for each leg to get per-leg totals.

    Args:
        leg_contributions: Dict of {leg_dir_key: probability}

    Returns:
        Dict of {leg_id: total_probability}
    """
    by_leg: dict[str, float] = {}
    for leg_dir_key, contrib in leg_contributions.items():
        # Parse leg key: "seg_id:direction:angle"
        parts = leg_dir_key.split(':')
        if parts:
            leg_id = parts[0]  # Just the segment ID
            by_leg[leg_id] = by_leg.get(leg_id, 0.0) + contrib
    return by_leg


def _build_initial_obstacle_dicts(
    by_object: dict, struct_lookup: dict, depth_lookup: dict,
) -> tuple[dict, dict]:
    allision_data: dict = {}
    grounding_data: dict = {}
    for obj_key, obj_vals in by_object.items():
        if not isinstance(obj_vals, dict):
            continue
        allision_prob = obj_vals.get('allision', 0.0)
        grounding_prob = obj_vals.get('grounding', 0.0)
        if obj_key.startswith('Structure - '):
            obj_id = obj_key.replace('Structure - ', '')
            if obj_id in struct_lookup and allision_prob > 0:
                meta = struct_lookup[obj_id]
                allision_data[obj_id] = {
                    'total_probability': allision_prob,
                    'geometry': meta.get('wkt_wgs84') or meta.get('wkt'),
                    'value': meta['height'],
                    'leg_contributions': {}, 'leg_dir_contributions': {}, 'segment_contributions': {},
                }
        elif obj_key.startswith('Depth - '):
            obj_id = obj_key.replace('Depth - ', '')
            if obj_id in depth_lookup and grounding_prob > 0:
                meta = depth_lookup[obj_id]
                grounding_data[obj_id] = {
                    'total_probability': grounding_prob,
                    'geometry': meta.get('wkt_wgs84') or meta.get('wkt'),
                    'value': meta['depth'],
                    'leg_contributions': {}, 'leg_dir_contributions': {}, 'segment_contributions': {},
                }
    return allision_data, grounding_data


def _fill_legdir_contributions(obstacle_data: dict, legdir_map: dict, prefix: str) -> None:
    for obs_key, leg_contribs in legdir_map.items():
        if not isinstance(leg_contribs, dict):
            continue
        obj_id = obs_key.replace(prefix, '')
        if obj_id in obstacle_data:
            obstacle_data[obj_id]['leg_dir_contributions'] = dict(leg_contribs)
            obstacle_data[obj_id]['leg_contributions'] = _aggregate_by_leg(leg_contribs)


def _fill_segment_contributions(obstacle_data: dict, seg_map: dict, prefix: str) -> None:
    for obs_key, seg_data in seg_map.items():
        if not isinstance(seg_data, dict):
            continue
        obj_id = obs_key.replace(prefix, '')
        if obj_id in obstacle_data:
            obstacle_data[obj_id]['segment_contributions'] = dict(seg_data)


def _fill_segment_edge_meta(obstacle_data: dict, edge_meta: dict, prefix: str) -> None:
    """Attach per-segment edge_calc_meta to the obstacle dicts.

    Structure of ``edge_meta`` is
    ``{obstacle_key: {seg_key: {leg_dir_key: {edge_dist_m, reachable_width_m,
    edge_p_nr, edge_h_eff, contrib_sum, ...}}}}``.  For each obstacle we
    condense every (seg, leg_dir) entry into a single per-segment attribute
    set using contribution-weighted averages (so segments hit from many
    directions still report one representative value per attribute).
    """
    for obs_key, seg_map in (edge_meta or {}).items():
        if not obs_key.startswith(prefix):
            continue
        obj_id = obs_key.replace(prefix, '')
        if obj_id not in obstacle_data:
            continue
        per_seg: dict[str, dict[str, float]] = {}
        for seg_key, ld_map in (seg_map or {}).items():
            total_w = 0.0
            wsum = {'edge_dist_m': 0.0, 'reachable_width_m': 0.0,
                    'edge_p_nr': 0.0, 'edge_h_eff': 0.0}
            min_dist = None
            max_reach = 0.0
            for ld_key, entry in (ld_map or {}).items():
                w = float(entry.get('contrib_sum', 0.0))
                if w <= 0.0:
                    continue
                for k in wsum:
                    wsum[k] += float(entry.get(k, 0.0)) * w
                total_w += w
                ed = float(entry.get('edge_dist_m', 0.0))
                if min_dist is None or ed < min_dist:
                    min_dist = ed
                rw = float(entry.get('reachable_width_m', 0.0))
                if rw > max_reach:
                    max_reach = rw
            if total_w > 0.0:
                per_seg[seg_key] = {
                    'edge_dist_m': wsum['edge_dist_m'] / total_w,
                    'reachable_width_m': wsum['reachable_width_m'] / total_w,
                    'edge_p_nr': wsum['edge_p_nr'] / total_w,
                    'edge_h_eff': wsum['edge_h_eff'] / total_w,
                    'min_edge_dist_m': float(min_dist) if min_dist is not None else 0.0,
                    'max_reachable_width_m': max_reach,
                    'contrib_sum': total_w,
                }
        if per_seg:
            obstacle_data[obj_id]['segment_edge_meta'] = per_seg


def extract_obstacle_probabilities(
    report: dict[str, Any],
    structures: list[dict[str, Any]],
    depths: list[dict[str, Any]],
) -> tuple[dict[str, dict[str, Any]], dict[str, dict[str, Any]]]:
    """
    Extract probability contributions per obstacle from the drifting report.

    Args:
        report: The drifting_report from Calculation
        structures: List of structure metadata dicts with 'id', 'height', 'wkt', 'wkt_wgs84'
        depths: List of depth metadata dicts with 'id', 'depth', 'wkt', 'wkt_wgs84'

    Returns:
        (allision_data, grounding_data) where each is a dict:
        {
            'obstacle_id': {
                'total_probability': float,
                'geometry': BaseGeometry (in WGS84),
                'value': float (height or depth),
                'leg_contributions': {'leg_id': float, ...},  # Aggregated by leg
                'leg_dir_contributions': {'leg_dir_key': float, ...},  # Raw leg-direction data
                'segment_contributions': {'seg_idx': {'leg_dir_key': float, ...}, ...}  # Per-segment
            }
        }
    """
    struct_lookup = {s['id']: s for s in structures}
    depth_lookup = {d['id']: d for d in depths}
    allision_data, grounding_data = _build_initial_obstacle_dicts(
        report.get('by_object', {}), struct_lookup, depth_lookup,
    )
    _fill_legdir_contributions(allision_data, report.get('by_structure_legdir', {}), 'Structure - ')
    _fill_legdir_contributions(grounding_data, report.get('by_depth_legdir', {}), 'Depth - ')
    _fill_segment_contributions(allision_data, report.get('by_structure_segment_legdir', {}), 'Structure - ')
    _fill_segment_contributions(grounding_data, report.get('by_depth_segment_legdir', {}), 'Depth - ')
    edge_meta = report.get('edge_calc_meta', {})
    _fill_segment_edge_meta(allision_data, edge_meta, 'Structure - ')
    _fill_segment_edge_meta(grounding_data, edge_meta, 'Depth - ')
    return allision_data, grounding_data


def _extract_line_segments_with_normals(geom) -> list[tuple[float, float, float, float, float]]:
    """
    Extract individual line segments from a polygon boundary with their normal angles.

    For a rectangular polygon, this returns 4 separate LineStrings with normals.
    Each segment connects two consecutive vertices.

    IMPORTANT: This function normalizes polygon orientation to CCW (counter-clockwise)
    before extracting segments. This ensures consistent outward normal calculation.

    For CCW polygons:
    - Exterior ring goes counter-clockwise
    - Interior (hole) rings go clockwise
    - Outward normal = rotate segment vector 90° clockwise (right-hand rule)

    Args:
        geom: A shapely geometry (Polygon, MultiPolygon, LineString, etc.)

    Returns:
        List of (x1, y1, x2, y2, normal_angle) tuples representing line segments
        where normal_angle is the outward-facing direction in compass degrees
    """
    from shapely.geometry import Polygon, MultiPolygon, LineString, LinearRing
    from shapely.geometry import polygon as shapely_polygon

    segments: list[tuple[float, float, float, float, float]] = []

    def extract_from_ring(ring_coords):
        """Extract segments from a ring's coordinates."""
        coords = list(ring_coords)
        for i in range(len(coords) - 1):
            x1, y1 = coords[i][:2]
            x2, y2 = coords[i + 1][:2]
            # Skip zero-length segments
            if (x1, y1) != (x2, y2):
                normal = _segment_normal_angle(x1, y1, x2, y2)
                segments.append((x1, y1, x2, y2, normal))

    if isinstance(geom, Polygon):
        # Normalize polygon to CCW exterior, CW holes using shapely's orient()
        # This ensures consistent outward normal calculation
        oriented_geom = shapely_polygon.orient(geom, sign=1.0)  # 1.0 = CCW exterior
        # Extract exterior ring segments
        extract_from_ring(oriented_geom.exterior.coords)
        # Extract interior rings (holes) if any
        for interior in oriented_geom.interiors:
            extract_from_ring(interior.coords)
    elif isinstance(geom, MultiPolygon):
        for poly in geom.geoms:
            segments.extend(_extract_line_segments_with_normals(poly))
    elif isinstance(geom, (LineString, LinearRing)):
        extract_from_ring(geom.coords)
    elif hasattr(geom, 'boundary'):
        # For other geometry types, try to get boundary
        boundary = geom.boundary
        if hasattr(boundary, 'coords'):
            extract_from_ring(boundary.coords)
        elif hasattr(boundary, 'geoms'):
            # MultiLineString boundary
            for line in boundary.geoms:
                extract_from_ring(line.coords)

    return segments


def _calculate_segment_probability(
    segment_normal: float,
    leg_dir_contributions: dict[str, float],
) -> tuple[float, dict[str, float]]:
    """
    Calculate the probability for a specific segment based on its orientation.

    Each segment's probability is the sum of contributions from each leg-direction,
    weighted by how exposed the segment is to drift from that direction.

    Args:
        segment_normal: Outward normal angle of segment (OMRAT degrees: 0=N, 90=W, 180=S, 270=E)
        leg_dir_contributions: Dict of {leg_dir_key: probability}
            where leg_dir_key is "seg_id:direction_label:angle" (e.g., "1:NNW:45")

    Returns:
        (total_probability, per_leg_contributions)
        where per_leg_contributions is {leg_id: probability} aggregated by leg
    """
    segment_total = 0.0
    per_leg: dict[str, float] = {}

    for leg_dir_key, contrib in leg_dir_contributions.items():
        # Parse drift angle from key (e.g., "1:NNW:45" -> 45.0)
        drift_angle = _parse_angle_from_key(leg_dir_key)
        if drift_angle is None:
            continue

        # Calculate exposure factor for this segment to this drift direction
        exposure = _exposure_factor(segment_normal, drift_angle)

        # Segment contribution = original contribution × exposure factor
        segment_contrib = contrib * exposure
        segment_total += segment_contrib

        # Aggregate by leg
        parts = leg_dir_key.split(':')
        if parts:
            leg_id = parts[0]
            per_leg[leg_id] = per_leg.get(leg_id, 0.0) + segment_contrib

    return segment_total, per_leg


def _collect_all_leg_keys(obstacle_data: dict[str, dict[str, Any]]) -> set[str]:
    all_leg_keys: set[str] = set()
    for obs_data in obstacle_data.values():
        all_leg_keys.update(obs_data.get('leg_contributions', {}).keys())
    return all_leg_keys


def _init_result_layer(
    name: str, all_leg_keys: set[str],
) -> tuple[QgsVectorLayer, Any, dict[str, str]]:
    layer = QgsVectorLayer("LineString?crs=epsg:4326", name, "memory")
    provider = layer.dataProvider()
    if provider is None:
        logger.error(f"Failed to get data provider for layer {name}")
        return layer, None, {}
    fields = [
        QgsField("obstacle_id", QVariant.String),
        QgsField("segment_idx", QVariant.Int),
        QgsField("total_edge_probability", QVariant.Double),
        QgsField("object_probability", QVariant.Double),
        QgsField("normal_deg", QVariant.Double),
        QgsField("value", QVariant.Double),
        # Per-edge back-tracing attributes (contribution-weighted across
        # leg-directions and wind angles that struck this segment).
        QgsField("edge_dist_m", QVariant.Double),
        QgsField("reach_width_m", QVariant.Double),
        QgsField("edge_p_nr", QVariant.Double),
        QgsField("edge_h_eff", QVariant.Double),
        # Extremes across leg-directions for filtering / worst-case colouring.
        QgsField("min_edge_dist_m", QVariant.Double),
        QgsField("max_reach_width_m", QVariant.Double),
    ]
    leg_key_to_field: dict[str, str] = {}
    for leg_key in sorted(all_leg_keys):
        safe_name = f"leg_{leg_key}"
        leg_key_to_field[leg_key] = safe_name
        fields.append(QgsField(safe_name, QVariant.Double))
    provider.addAttributes(fields)
    layer.updateFields()
    _set_result_field_aliases(layer)
    _enable_fid_labels(layer, 'obstacle_id')
    return layer, provider, leg_key_to_field


def _resolve_seg_probs(
    seg_idx: int,
    segment_contributions: dict,
    obs_total_prob: float,
    obs_leg_contribs: dict,
) -> tuple[float, dict[str, float]]:
    if not segment_contributions:
        return obs_total_prob, obs_leg_contribs
    seg_legdir_contribs = segment_contributions.get(f"seg_{seg_idx}", {})
    seg_prob = sum(seg_legdir_contribs.values()) if seg_legdir_contribs else 0.0
    seg_leg_contribs: dict[str, float] = {}
    for leg_dir_key, contrib in seg_legdir_contribs.items():
        parts = leg_dir_key.split(':')
        if parts:
            leg_id = parts[0]
            seg_leg_contribs[leg_id] = seg_leg_contribs.get(leg_id, 0.0) + contrib
    return seg_prob, seg_leg_contribs


def _make_obstacle_segment_feat(
    layer: QgsVectorLayer,
    obs_id: str,
    obs_data: dict,
    seg_idx: int,
    x1: float, y1: float, x2: float, y2: float,
    normal_angle: float,
    seg_prob: float,
    seg_leg_contribs: dict[str, float],
    leg_key_to_field: dict[str, str],
) -> QgsFeature:
    feat = QgsFeature(layer.fields())
    feat.setGeometry(QgsGeometry.fromWkt(f"LINESTRING({x1} {y1}, {x2} {y2})"))
    feat.setAttribute("obstacle_id", obs_id)
    feat.setAttribute("segment_idx", seg_idx)
    feat.setAttribute("total_edge_probability", seg_prob)
    feat.setAttribute("object_probability", obs_data.get('total_probability', 0.0))
    feat.setAttribute("normal_deg", normal_angle)
    feat.setAttribute("value", obs_data.get('value', 0.0))
    # Per-edge back-tracing attributes from edge_calc_meta (or 0 if the
    # segment never entered the per-edge branch -- e.g. simple-branch obstacles).
    seg_edge_meta = obs_data.get('segment_edge_meta', {}) or {}
    em = seg_edge_meta.get(f"seg_{seg_idx}", {})
    feat.setAttribute("edge_dist_m", float(em.get('edge_dist_m', 0.0)))
    feat.setAttribute("reach_width_m", float(em.get('reachable_width_m', 0.0)))
    feat.setAttribute("edge_p_nr", float(em.get('edge_p_nr', 0.0)))
    feat.setAttribute("edge_h_eff", float(em.get('edge_h_eff', 0.0)))
    feat.setAttribute("min_edge_dist_m", float(em.get('min_edge_dist_m', 0.0)))
    feat.setAttribute("max_reach_width_m", float(em.get('max_reachable_width_m', 0.0)))
    for leg_key, field_name in leg_key_to_field.items():
        feat.setAttribute(field_name, seg_leg_contribs.get(leg_key, 0.0))
    return feat


def _build_features_for_obstacle(
    layer: QgsVectorLayer,
    obs_id: str,
    obs_data: dict[str, Any],
    leg_key_to_field: dict[str, str],
) -> list[QgsFeature]:
    geom = obs_data.get('geometry')
    if geom is None:
        return []
    try:
        segments = _extract_line_segments_with_normals(geom)
    except Exception as e:
        logger.warning(f"Failed to extract segments for {obs_id}: {e}")
        return []
    if not segments:
        logger.warning(f"No segments extracted for {obs_id}")
        return []
    obs_total_prob = obs_data.get('total_probability', 0.0)
    segment_contributions = obs_data.get('segment_contributions', {})
    obs_leg_contribs = obs_data.get('leg_contributions', {})
    feats = []
    for seg_idx, (x1, y1, x2, y2, normal_angle) in enumerate(segments):
        seg_prob, seg_leg_contribs = _resolve_seg_probs(
            seg_idx, segment_contributions, obs_total_prob, obs_leg_contribs,
        )
        feats.append(_make_obstacle_segment_feat(
            layer, obs_id, obs_data, seg_idx,
            x1, y1, x2, y2, normal_angle,
            seg_prob, seg_leg_contribs, leg_key_to_field,
        ))
    return feats


def _log_result_layer_stats(
    features: list, layer_type: str, obstacle_data: dict,
) -> None:
    probs = [
        f["total_edge_probability"] for f in features
        if f["total_edge_probability"] is not None
    ]
    if probs:
        logger.info(
            f"Created {layer_type} layer with {len(features)} segments from "
            f"{len(obstacle_data)} obstacles. Probability range: {min(probs):.2e} - {max(probs):.2e}"
        )
    else:
        logger.info(
            f"Created {layer_type} layer with {len(features)} segments from {len(obstacle_data)} obstacles"
        )


def create_result_layer(
    name: str,
    obstacle_data: dict[str, dict[str, Any]],
    layer_type: str = 'allision',
) -> QgsVectorLayer | None:
    """
    Create a QGIS vector layer with obstacle geometries and probability attributes.

    Each obstacle's boundary is split into individual line segments. All segments
    of an obstacle share the same probability values since the underlying calculation
    computes probability per obstacle (not per segment).

    Geometries are expected to already be in WGS84.

    Args:
        name: Layer name
        obstacle_data: Dict from extract_obstacle_probabilities (geometries in WGS84)
        layer_type: 'allision' or 'grounding'

    Returns:
        QgsVectorLayer with features and attributes in WGS84
    """
    if not obstacle_data:
        logger.info(f"No {layer_type} data to create layer")
        return None
    all_leg_keys = _collect_all_leg_keys(obstacle_data)
    layer, provider, leg_key_to_field = _init_result_layer(name, all_leg_keys)
    if provider is None:
        return None
    features: list[QgsFeature] = []
    for obs_id, obs_data in obstacle_data.items():
        features.extend(_build_features_for_obstacle(layer, obs_id, obs_data, leg_key_to_field))
    provider.addFeatures(features)
    _log_result_layer_stats(features, layer_type, obstacle_data)
    return layer


def _make_graduated_ranges(
    min_val: float, max_val: float, num_classes: int,
) -> list[QgsRendererRange]:
    colors = [
        QColor(0, 255, 0),
        QColor(128, 255, 0),
        QColor(255, 255, 0),
        QColor(255, 128, 0),
        QColor(255, 0, 0),
    ]
    step = (max_val - min_val) / num_classes
    ranges: list[QgsRendererRange] = []
    for i in range(num_classes):
        lower = min_val + i * step
        upper = min_val + (i + 1) * step
        color = colors[min(i, len(colors) - 1)]
        symbol = QgsLineSymbol.createSimple({
            'color': color.name(),
            'width': str(0.5 + 0.3 * i),
        })
        ranges.append(QgsRendererRange(lower, upper, symbol, f"{lower:.2e} - {upper:.2e}"))
    return ranges


def apply_graduated_symbology(
    layer: QgsVectorLayer,
    attribute: str = "total_edge_probability",
    num_classes: int = 5,
) -> None:
    """
    Apply graduated symbology to a layer: green (low) -> yellow (medium) -> red (high).

    A grey "Not hit" class [0, 1e-18] is prepended so zero-contribution
    obstacle edges (e.g. segments shielded by shallower obstacles) are
    visually distinct from segments with actual probability.

    Args:
        layer: The vector layer to style
        attribute: The attribute field to classify by
        num_classes: Number of classification bins
    """
    if layer is None or layer.featureCount() == 0:
        return
    if layer.fields().indexOf(attribute) < 0:
        logger.warning(f"Attribute {attribute} not found in layer")
        return

    not_hit_symbol = QgsLineSymbol.createSimple({
        'color': '160,160,160,140',
        'width': '0.3',
    })
    not_hit_range = QgsRendererRange(0.0, 1e-18, not_hit_symbol, "Not hit")

    values = [f[attribute] for f in layer.getFeatures() if f[attribute] is not None and f[attribute] > 0]
    if not values:
        renderer = QgsGraduatedSymbolRenderer(attribute, [not_hit_range])
        renderer.setMode(QgsGraduatedSymbolRenderer.Mode.Custom)
        layer.setRenderer(renderer)
        layer.triggerRepaint()
        return
    min_val, max_val = min(values), max(values)
    if min_val >= max_val:
        symbol = QgsLineSymbol.createSimple({'color': 'yellow', 'width': '0.8'})
        from qgis.core import QgsSingleSymbolRenderer
        layer.setRenderer(QgsSingleSymbolRenderer(symbol))
        layer.triggerRepaint()
        return
    ranges = [not_hit_range] + _make_graduated_ranges(min_val, max_val, num_classes)
    renderer = QgsGraduatedSymbolRenderer(attribute, ranges)
    renderer.setMode(QgsGraduatedSymbolRenderer.Mode.Custom)
    layer.setRenderer(renderer)
    layer.triggerRepaint()


def create_result_layers(
    report: dict[str, Any],
    structures: list[dict[str, Any]],
    depths: list[dict[str, Any]],
    add_to_project: bool = True,
) -> tuple[QgsVectorLayer | None, QgsVectorLayer | None]:
    """
    Create grounding and allision result layers from calculation results.

    The geometries are expected to be in WGS84 (from 'wkt_wgs84' field in metadata).

    Args:
        report: The drifting_report from Calculation
        structures: List of structure metadata dicts with 'wkt_wgs84'
        depths: List of depth metadata dicts with 'wkt_wgs84'
        add_to_project: If True, add layers to the current QGIS project

    Returns:
        (allision_layer, grounding_layer)
    """
    if report is None:
        logger.warning("No report data provided")
        return None, None

    # Extract probability data (uses WGS84 geometries)
    allision_data, grounding_data = extract_obstacle_probabilities(
        report, structures, depths
    )

    # Create layers (geometries are already in WGS84)
    allision_layer = create_result_layer(
        "Allision Results",
        allision_data,
        "allision",
    )

    grounding_layer = create_result_layer(
        "Grounding Results",
        grounding_data,
        "grounding",
    )

    # Apply symbology
    if allision_layer is not None:
        apply_graduated_symbology(allision_layer)
        if add_to_project:
            QgsProject.instance().addMapLayer(allision_layer)
            logger.info("Added Allision Results layer to project")

    if grounding_layer is not None:
        apply_graduated_symbology(grounding_layer)
        if add_to_project:
            QgsProject.instance().addMapLayer(grounding_layer)
            logger.info("Added Grounding Results layer to project")

    return allision_layer, grounding_layer


# ---------------------------------------------------------------------------
# Ship-collision result layers (lines on legs, points at waypoints)
# ---------------------------------------------------------------------------

def _color_ramp(num_classes: int = 5) -> list[QColor]:
    """Standard green->red ramp shared by every result layer."""
    return [
        QColor(0, 255, 0),     # Green - lowest
        QColor(128, 255, 0),
        QColor(255, 255, 0),   # Yellow - medium
        QColor(255, 128, 0),
        QColor(255, 0, 0),     # Red - highest
    ][:num_classes]


def _apply_graduated(
    layer: QgsVectorLayer,
    attribute: str,
    symbol_factory,  # callable: (QColor, fraction in [0,1]) -> QgsSymbol
    num_classes: int = 5,
) -> None:
    """Generic graduated-symbology helper for line / point / polygon layers.

    ``symbol_factory(color, frac)`` is called once per class to build the
    range symbol; ``frac`` runs from 0 (lowest class) to 1 (highest) so
    callers can vary line width / marker size by class.
    """
    if layer is None or layer.featureCount() == 0:
        return
    if layer.fields().indexOf(attribute) < 0:
        logger.warning(f"Attribute {attribute} not found on {layer.name()}")
        return
    values = [
        f[attribute] for f in layer.getFeatures()
        if f[attribute] is not None and f[attribute] > 0
    ]
    if not values:
        return
    min_val, max_val = min(values), max(values)
    if min_val >= max_val:
        symbol = symbol_factory(QColor(255, 255, 0), 0.5)
        from qgis.core import QgsSingleSymbolRenderer
        layer.setRenderer(QgsSingleSymbolRenderer(symbol))
        layer.triggerRepaint()
        return
    colors = _color_ramp(num_classes)
    step = (max_val - min_val) / num_classes
    ranges = []
    for i in range(num_classes):
        lower = min_val + i * step
        upper = min_val + (i + 1) * step
        frac = (i / max(num_classes - 1, 1)) if num_classes > 1 else 0.0
        symbol = symbol_factory(colors[min(i, len(colors) - 1)], frac)
        ranges.append(QgsRendererRange(
            lower, upper, symbol, f"{lower:.2e} - {upper:.2e}",
        ))
    renderer = QgsGraduatedSymbolRenderer(attribute, ranges)
    renderer.setMode(QgsGraduatedSymbolRenderer.Mode.Custom)
    layer.setRenderer(renderer)
    layer.triggerRepaint()


def _line_symbol_factory(color: QColor, frac: float) -> QgsLineSymbol:
    """Wider line symbol with transparency, scaled by class fraction."""
    width = 1.4 + 1.6 * frac      # 1.4mm at low-class, 3.0mm at top-class
    rgba = f"{color.red()},{color.green()},{color.blue()},140"  # ~55% alpha
    return QgsLineSymbol.createSimple({
        'color': rgba,
        'width': f'{width}',
        'capstyle': 'round',
    })


def _marker_symbol_factory(color: QColor, frac: float) -> QgsMarkerSymbol:
    """Bigger circle for higher-class points."""
    size = 3.0 + 4.0 * frac
    return QgsMarkerSymbol.createSimple({
        'name': 'circle',
        'color': color.name(),
        'outline_color': '50,50,50,255',
        'outline_width': '0.4',
        'size': f'{size}',
    })


def _build_collision_line_layer(
    by_leg: dict, segment_data: dict,
) -> QgsVectorLayer:
    layer = QgsVectorLayer('LineString?crs=EPSG:4326', 'Ship-Ship Collision (per leg)', 'memory')
    lp = layer.dataProvider()
    lp.addAttributes([
        QgsField('leg_id', QVariant.String),
        QgsField('head_on', QVariant.Double),
        QgsField('overtaking', QVariant.Double),
        QgsField('combined', QVariant.Double),
    ])
    layer.updateFields()
    feats: list[QgsFeature] = []
    for leg_id, contribs in by_leg.items():
        seg = segment_data.get(leg_id, {})
        sp_str, ep_str = seg.get('Start_Point', ''), seg.get('End_Point', '')
        if not sp_str or not ep_str:
            continue
        try:
            sp = [float(x) for x in str(sp_str).split()][:2]
            ep = [float(x) for x in str(ep_str).split()][:2]
        except (TypeError, ValueError):
            continue
        head_on = float(contribs.get('head_on', 0.0) or 0.0)
        overtaking = float(contribs.get('overtaking', 0.0) or 0.0)
        combined = head_on + overtaking
        if combined <= 0.0:
            continue
        feat = QgsFeature(layer.fields())
        feat.setGeometry(QgsGeometry.fromPolylineXY([QgsPointXY(sp[0], sp[1]), QgsPointXY(ep[0], ep[1])]))
        feat.setAttribute('leg_id', str(leg_id))
        feat.setAttribute('head_on', head_on)
        feat.setAttribute('overtaking', overtaking)
        feat.setAttribute('combined', combined)
        feats.append(feat)
    lp.addFeatures(feats)
    _apply_graduated(layer, 'combined', _line_symbol_factory)
    return layer


def _build_collision_point_layer(by_waypoint: dict) -> QgsVectorLayer:
    layer = QgsVectorLayer('Point?crs=EPSG:4326', 'Ship-Ship Collision (waypoints)', 'memory')
    pp = layer.dataProvider()
    pp.addAttributes([
        QgsField('waypoint', QVariant.String),
        QgsField('crossing', QVariant.Double),
        QgsField('merging', QVariant.Double),
        QgsField('bend', QVariant.Double),
        QgsField('combined', QVariant.Double),
    ])
    layer.updateFields()
    feats: list[QgsFeature] = []
    for wp_key, rec in by_waypoint.items():
        try:
            lon_str, lat_str = wp_key.split()
            lon, lat = float(lon_str), float(lat_str)
        except Exception:  # nosec B110 B112
            continue
        crossing = float(rec.get('crossing', 0.0) or 0.0)
        merging = float(rec.get('merging', 0.0) or 0.0)
        bend = float(rec.get('bend', 0.0) or 0.0)
        combined = crossing + merging + bend
        if combined <= 0.0:
            continue
        feat = QgsFeature(layer.fields())
        feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(lon, lat)))
        feat.setAttribute('waypoint', wp_key)
        feat.setAttribute('crossing', crossing)
        feat.setAttribute('merging', merging)
        feat.setAttribute('bend', bend)
        feat.setAttribute('combined', combined)
        feats.append(feat)
    pp.addFeatures(feats)
    _apply_graduated(layer, 'combined', _marker_symbol_factory)
    return layer


def _add_collision_layers_to_project(
    line_layer: QgsVectorLayer, point_layer: QgsVectorLayer,
) -> None:
    if line_layer.featureCount() > 0:
        QgsProject.instance().addMapLayer(line_layer)
        logger.info(f"Added Ship-Ship Collision (per leg) layer with {line_layer.featureCount()} features")
    if point_layer.featureCount() > 0:
        QgsProject.instance().addMapLayer(point_layer)
        logger.info(f"Added Ship-Ship Collision (waypoints) layer with {point_layer.featureCount()} features")


def create_collision_layers(
    collision_report: dict[str, Any] | None,
    segment_data: dict[str, Any],
    add_to_project: bool = True,
) -> tuple[QgsVectorLayer | None, QgsVectorLayer | None]:
    """Build the two ship-ship collision result layers.

    Returns ``(line_layer, point_layer)``:

    * ``line_layer`` is a per-leg line layer with one feature per leg,
      attributes ``head_on``, ``overtaking``, and ``combined``
      (head-on + overtaking).
    * ``point_layer`` is a point layer with one feature per shared
      waypoint, attributes ``crossing``, ``bend``, ``combined``
      (crossing + bend).

    Both layers are styled graduated red->green by ``combined``.  The
    line symbol is intentionally wider than the digitised route line
    and semi-transparent so the underlying route stays visible.
    """
    if collision_report is None:
        return None, None
    by_leg = collision_report.get('by_leg', {}) or {}
    by_waypoint = collision_report.get('by_waypoint', {}) or {}
    line_layer = _build_collision_line_layer(by_leg, segment_data)
    point_layer = _build_collision_point_layer(by_waypoint)
    if add_to_project:
        _add_collision_layers_to_project(line_layer, point_layer)
    return line_layer, point_layer


# ---------------------------------------------------------------------------
# Powered (Cat II) result layers (obstacle boundary segments by probability)
# ---------------------------------------------------------------------------


def _parse_lonlat(text: str | None) -> tuple[float, float] | None:
    """Parse a ``"lon lat"`` or ``"Point (lon lat)"`` coordinate string."""
    if not isinstance(text, str):
        return None
    s = text.strip()
    if '(' in s and ')' in s:
        s = s.split('(', 1)[1].split(')', 1)[0]
    parts = s.replace(',', ' ').split()
    if len(parts) < 2:
        return None
    try:
        return float(parts[0]), float(parts[1])
    except ValueError:
        return None


def _seg_bearing(seg: dict) -> float | None:
    """Compass bearing (0=N, 90=E, CW) of a route leg from Start_Point to End_Point."""
    sp = _parse_lonlat(str(seg.get('Start_Point', '') or seg.get('Start Point', '')))
    ep = _parse_lonlat(str(seg.get('End_Point', '') or seg.get('End Point', '')))
    if sp is None or ep is None:
        return None
    dx = ep[0] - sp[0]   # eastward (longitude difference)
    dy = ep[1] - sp[1]   # northward (latitude difference)
    if dx == 0.0 and dy == 0.0:
        return None
    return float(np.degrees(np.arctan2(dx, dy)) % 360)


def _build_powered_layer(
    name: str,
    by_obstacle: dict[str, float],
    obstacle_records: list[dict[str, Any]],
    value_key: str,            # 'depth' or 'height'
    add_to_project: bool,
    by_obstacle_leg: dict[str, dict[str, float]] | None = None,
    segment_data: dict[str, dict[str, Any]] | None = None,
) -> QgsVectorLayer | None:
    """Obstacle boundary-segment layer coloured by powered probability.

    Each feature is a boundary segment of a depth contour or structure
    polygon, coloured by that obstacle's total powered probability.
    Mirrors the drifting grounding/allision result layer layout.
    """
    if not by_obstacle:
        return None
    return _build_powered_layer_per_obstacle_segments(
        name=name,
        by_obstacle=by_obstacle,
        obstacle_records=obstacle_records,
        value_key=value_key,
        by_obstacle_leg=by_obstacle_leg,
        add_to_project=add_to_project,
        segment_data=segment_data,
    )


def _build_powered_layer_per_obstacle_segments(
    *,
    name: str,
    by_obstacle: dict[str, float],
    obstacle_records: list[dict[str, Any]],
    value_key: str,
    by_obstacle_leg: dict[str, dict[str, float]] | None,
    add_to_project: bool,
    segment_data: dict[str, dict[str, Any]] | None = None,
) -> QgsVectorLayer | None:
    """Per-obstacle boundary-segment layer for powered grounding/allision.

    Each boundary segment of the obstacle polygon becomes one LineString
    feature.  When ``segment_data`` is available, the segment probability
    is the obstacle's total leg contribution weighted by the cosine exposure
    factor between the leg bearing and the segment's outward normal —
    matching the approach used by the drifting grounding layer.  Without
    ``segment_data`` every segment gets the obstacle's total probability.
    """
    import shapely.wkt as shapely_wkt

    all_leg_ids, leg_to_field = _collect_leg_fields(by_obstacle_leg)
    layer = QgsVectorLayer('LineString?crs=EPSG:4326', name, 'memory')
    pr = layer.dataProvider()
    fields = [
        QgsField('obstacle_id', QVariant.String),
        QgsField('segment_idx', QVariant.Int),
        QgsField('total_edge_probability', QVariant.Double),
        QgsField('object_probability', QVariant.Double),
        QgsField('value', QVariant.Double),
    ]
    for _leg_id, field_name in leg_to_field.items():
        fields.append(QgsField(field_name, QVariant.Double))
    pr.addAttributes(fields)
    layer.updateFields()
    _set_result_field_aliases(layer)
    _enable_fid_labels(layer, 'obstacle_id')

    lookup = {str(rec.get('id', '')): rec for rec in obstacle_records}
    feats: list[QgsFeature] = []
    for obs_id, prob in by_obstacle.items():
        if prob <= 0.0:
            continue
        meta = lookup.get(str(obs_id))
        if meta is None:
            continue
        geom_raw = meta.get('wkt_wgs84') or meta.get('wkt')
        if geom_raw is None:
            continue
        try:
            geom_str = geom_raw.wkt if hasattr(geom_raw, 'wkt') else str(geom_raw)
            shapely_geom = shapely_wkt.loads(geom_str)
        except Exception:  # nosec B110 B112
            continue
        segments = _extract_line_segments_with_normals(shapely_geom)
        if not segments:
            continue
        try:
            value = float(meta.get(value_key, 0.0) or 0.0)
        except (TypeError, ValueError):
            value = 0.0
        obs_leg_contribs = (by_obstacle_leg or {}).get(str(obs_id), {})

        # Pre-compute bearings for each leg-direction key ("seg_id:dir_idx").
        # dir_idx 0 = forward (Start→End); dir_idx 1 = reverse (+180°).
        leg_bearings: dict[str, float | None] = {}
        if segment_data and obs_leg_contribs:
            for leg_dir_key in obs_leg_contribs:
                parts = leg_dir_key.rsplit(':', 1)
                seg_key = parts[0]
                dir_idx_val = int(parts[1]) if len(parts) == 2 and parts[1].isdigit() else 0
                seg = segment_data.get(seg_key) or segment_data.get(str(seg_key))
                bearing = _seg_bearing(seg) if seg is not None else None
                if bearing is not None and dir_idx_val == 1:
                    bearing = (bearing + 180.0) % 360.0
                leg_bearings[leg_dir_key] = bearing

        for seg_idx, (x1, y1, x2, y2, seg_normal) in enumerate(segments):
            # Per-segment probability: sum each leg's contribution weighted by
            # how directly the leg bearing faces this segment edge.
            if segment_data and obs_leg_contribs:
                seg_prob = 0.0
                for leg_dir_key, leg_contrib in obs_leg_contribs.items():
                    bearing = leg_bearings.get(leg_dir_key)
                    if bearing is None:
                        seg_prob += leg_contrib  # unknown bearing: add uniformly
                    else:
                        seg_prob += leg_contrib * _exposure_factor(seg_normal, bearing)
            else:
                seg_prob = float(prob)

            feat = QgsFeature(layer.fields())
            feat.setGeometry(QgsGeometry.fromWkt(f'LINESTRING({x1} {y1}, {x2} {y2})'))
            feat.setAttribute('obstacle_id', str(obs_id))
            feat.setAttribute('segment_idx', seg_idx)
            feat.setAttribute('total_edge_probability', seg_prob)
            feat.setAttribute('object_probability', float(prob))
            feat.setAttribute('value', value)
            for leg_id, field_name in leg_to_field.items():
                feat.setAttribute(field_name, float(obs_leg_contribs.get(leg_id, 0.0) or 0.0))
            feats.append(feat)

    if not feats:
        return None
    pr.addFeatures(feats)
    _apply_graduated(layer, 'total_edge_probability', _line_symbol_factory)
    if add_to_project and layer.featureCount() > 0:
        QgsProject.instance().addMapLayer(layer)
        logger.info(f"Added {name} layer with {layer.featureCount()} features (per obstacle segment)")
    return layer


def _collect_leg_fields(
    by_obstacle_leg: dict | None,
) -> tuple[list[str], dict[str, str]]:
    all_leg_ids: list[str] = []
    leg_to_field: dict[str, str] = {}
    if by_obstacle_leg:
        seen: set[str] = set()
        for legs_for_obs in by_obstacle_leg.values():
            for leg_id in legs_for_obs:
                if leg_id not in seen:
                    seen.add(leg_id)
                    all_leg_ids.append(str(leg_id))
        all_leg_ids.sort()
        for leg_id in all_leg_ids:
            leg_to_field[leg_id] = f'leg_{leg_id.replace(":", "_")}'
    return all_leg_ids, leg_to_field


def create_powered_grounding_layer(
    powered_grounding_report: dict[str, Any] | None,
    depths: list[dict[str, Any]],
    add_to_project: bool = True,
    segment_data: dict[str, dict[str, Any]] | None = None,
) -> QgsVectorLayer | None:
    """Depth-contour boundary-segment layer coloured by powered-grounding probability.

    Each feature is a boundary segment of a depth contour polygon,
    coloured by that contour's total powered-grounding probability.
    Per-leg contributions are kept as ``leg_<id>`` columns.
    """
    if not powered_grounding_report:
        return None
    by_obstacle = powered_grounding_report.get('by_obstacle', {}) or {}
    if not by_obstacle:
        return None
    by_obstacle_leg = powered_grounding_report.get('by_obstacle_leg', {}) or {}
    return _build_powered_layer(
        'Powered Grounding Results',
        by_obstacle, depths, 'depth', add_to_project,
        by_obstacle_leg=by_obstacle_leg,
        segment_data=segment_data,
    )


def create_powered_allision_layer(
    powered_allision_report: dict[str, Any] | None,
    structures: list[dict[str, Any]],
    add_to_project: bool = True,
    segment_data: dict[str, dict[str, Any]] | None = None,
) -> QgsVectorLayer | None:
    """Structure boundary-segment layer coloured by powered-allision probability.

    Each feature is a boundary segment of a structure polygon, coloured
    by that structure's total powered-allision probability. Per-leg
    contributions are kept as ``leg_<id>`` columns.
    """
    if not powered_allision_report:
        return None
    by_obstacle = powered_allision_report.get('by_obstacle', {}) or {}
    if not by_obstacle:
        return None
    by_obstacle_leg = powered_allision_report.get('by_obstacle_leg', {}) or {}
    return _build_powered_layer(
        'Powered Allision Results',
        by_obstacle, structures, 'height', add_to_project,
        by_obstacle_leg=by_obstacle_leg,
        segment_data=segment_data,
    )
