Code Flow: Powered Grounding & Allision (Cat I + II)¶

This chapter walks OMRAT’s IWRAP Category I and Category II powered-grounding and powered-allision calculations one function at a time, in the order the calls fire when a user presses Run model. Pair it with Powered Grounding and Allision, which derives the Cat I formula \(N_I = P_c Q \cdot \mathrm{mass}\) and the Cat II formula \(N_{II} = P_c Q \cdot \mathrm{mass} \cdot \exp(-d_\mathrm{mean}/(a_i V))\) and explains the shadow model.

Entry points ¶

Powered grounding is phase 3 and powered allision is phase 4 of CalculationTask (see Code Flow: From “Run model” to Results). They share the same ray-cast engine and the same obstacle projection step, but the obstacle sources differ:

Grounding uses the depth polygons and filters by ship draught: only depths shallower than the ship’s draught are hazards. One computation is run per depth bin – a bin is the set of depths a draught-class “sees”, so all draughts in the same bin share one shadow-aware result.
Allision uses the object polygons and filters by ship height: if a ship is shorter than a structure it passes under (bridge case).

Both phases call geometries.get_powered_overlap._run_all_computations() to produce the shadow-aware ray-cast summaries, then iterate traffic to multiply by per-cell frequency.

Grounding entry: compute/powered_model.py:417 – run_powered_grounding_model()

Allision entry: compute/powered_model.py:458 – run_powered_allision_model()

Top-level call tree ¶

run_powered_grounding_model(data)
  |
  +-- _parse_point(first_seg["Start_Point"])
  +-- SimpleProjector(lon0, lat0)                                     # equirect projection origin
  +-- for each ship-draught, pick a depth bin
  +-- for each bin:
  |     _build_legs_and_obstacles(data, proj, mode="grounding", max_draft=bin)
  |       +-- sw.loads(wkt_str) per depth polygon
  |       +-- _project_wkt_geom(geom, proj)                           # shapely.ops.transform
  |       +-- _weighted_avg_speed_knots(seg_traffic[d_name])          # per direction avg speed
  |     _run_all_computations(legs, all_obstacles)
  |       +-- per (leg, dir):
  |             _leg_vectors(start, end)
  |             _compute_cat2_with_shadows(turn_pt, ext_dir, perp,    # vectorised ray cast
  |                                         mean, std, ai, V, obstacles)
  |               +-- _extract_edges_local(obs_geom, origin, u, u_perp)
  |               +-- numpy broadcast: (n_rays x n_edges) crossings
  |               +-- per-obstacle argmin + summary dict
  +-- for each traffic cell x direction:
        lookup comps by draught-bin; sum per-obstacle mass * exp(-d_mean/(ai*V))

run_powered_allision_model(data)
  |
  +-- SimpleProjector + _build_legs_and_obstacles(mode="allision")    # max_draft=0 flags objects only
  +-- _run_all_computations(legs, all_obstacles)                      # same as grounding
  +-- for each comp x traffic cell:
        filter by ship_height < obstacle_height (bridge clearance)
        sum * pc_allision

Shared inputs ¶

Both phases consume data:

Field	Use
`segment_data[leg]['Start_Point'] / 'End_Point'`	Leg endpoints in lon/lat. The first leg’s `Start_Point` defines the projector origin.
`segment_data[leg]['ai1'] / 'ai2'`	Position check interval (seconds) per direction. Plugged into the Cat II exponential: \(\exp(-d / (a_i V))\).
`segment_data[leg]['mean1_1'] / 'std1_1' / 'mean2_1' / 'std2_1'`	Per-direction lateral distribution. A single Gaussian is assumed (the first superposed component); multi-component distributions would need the analytical probability path (see Code Flow: Drifting Model).
`traffic_data[leg][dir]['Frequency (ships/year)'] / 'Speed (knots)'`	Per-cell frequency and speed; speed is averaged if the cell holds a list. Frequency is read after `compute.data_preparation.apply_traffic_scaling()` has applied the per-cell `Scaling (%)` multiplier (see the .omrat data format chapter).
`traffic_data[leg][dir]['Draught (meters)']` (grounding only)	Selects the depth bin.
`traffic_data[leg][dir]['Ship heights (meters)']` (allision only)	Filters out ships short enough to pass under the structure.
`objects` / `depths`	Obstacle geometries (WKT) + their `height` / `depth` attribute.
`pc.grounding` / `pc.allision`	Causation factors. \(N_{II} = P_c \cdot Q \cdot \ldots\).

Projection: `SimpleProjector`¶

The powered-model geometry is all done in a local equirectangular frame centred on the first leg’s start point. This is cheaper than UTM (no zone selection, no pyproj dependency) and fine for per-leg-scale distances since the rays are only cast a few tens of kilometres.

Source: geometries/get_powered_overlap.py:71 – SimpleProjector

class SimpleProjector:
    def __init__(self, lon_ref, lat_ref):
        self.mx = 111_320.0 * cos(radians(lat_ref))
        self.my = 110_540.0

    def transform(self, lon, lat):
        return ((lon - self.lon_ref) * self.mx,
                (lat - self.lat_ref) * self.my)

Used by _build_legs_and_obstacles() to project legs and every obstacle once.

`run_powered_grounding_model()`¶

Source: compute/powered_model.py:417 – run_powered_grounding_model()

Flow:

Short-circuit if traffic / segments / depths are empty.
Build a SimpleProjector from the first leg’s Start_Point.
Collect every unique non-zero draught across traffic. Defaults to {5.0} if no traffic has draught set.
Depth binning. For each unique draught, find the largest depth_value below it (_depth_bin_key). Two ships whose draughts fall in the same bin “see” the same obstacle set, so we avoid recomputing the shadow geometry. Bins are sorted so the progress bar advances monotonically; None is included if some draughts are below every depth polygon (no obstacles).
Per bin:
- _build_legs_and_obstacles(data, proj, mode='grounding', max_draft=bin) - returns (legs, all_obstacles, depth_geoms, ...). Only depths \(\le\) max_draft become obstacles.
- bin_results[bin] = _run_all_computations(legs, all_obstacles) - per (leg, direction) shadow-aware ray cast.
- Progress reported via the top-level progress callback with phase "grounding".

Iterate the traffic matrix:

for leg_key, leg_dirs in traffic_data.items():
    for dir_idx, (dir_key, dir_data) in enumerate(leg_dirs.items()):
        ai_seconds = ai_per_dir[min(dir_idx, 1)]
        for loa_i, freq_row in enumerate(dir_data['Frequency (ships/year)']):
            for type_j, q in enumerate(freq_row):
                if q <= 0: continue
                draught, speed = ...  # per-cell
                comps = bin_results[_depth_bin_key(draught)]
                for comp in comps:
                    if comp['seg_id'] != leg_key: continue
                    if comp['dir_idx'] != dir_idx: continue
                    for key, s in comp['summaries'].items():
                        mass, d_mean = s['mass'], s['mean_dist']
                        recovery = ai_seconds * speed_ms
                        total += pc_grounding * q * mass * exp(-d_mean / recovery)

Write LEPPoweredGrounding.setText(f"{total:.3e}") and return total.

`run_powered_allision_model()`¶

Source: compute/powered_model.py:458 – run_powered_allision_model()

Structurally identical to grounding but:

No depth binning. The obstacle set is fixed (objects), so _run_all_computations() is called exactly once.
Height filter. For each obstacle summary, the inner loop skips ship cells whose ship_height < obj_height: those pass under the structure (bridge case).
Uses pc.allision as the causation factor.

Everything else (projector build, leg + obstacle build, ray cast, cell iteration) is the same pattern as grounding.

`_build_legs_and_obstacles()`¶

Source: geometries/get_powered_overlap.py:615 – _build_legs_and_obstacles()

Purpose: given data + projector + mode (+ max_draft), produce all the geometry the ray cast needs in one local-frame dict.

Returns (legs, all_obstacles, depth_geoms, depth_geoms_deep, object_geoms):

legs[seg_id] is a dict with:

{'start', 'end', 'name', 'start_wkt', 'end_wkt', 'dirs': [{'name', 'speed_kn', 'speed_ms', 'ai', 'mean', 'std'}, ...]}. start / end are numpy (x, y) in metres.
all_obstacles is a list of (obstacle_dict, kind) where kind is 'depth' or 'object'. Obstacle dicts carry 'id', 'geom' (shapely Polygon in local frame), and either 'depth' or 'height'.
depth_geoms / depth_geoms_deep split the depths at max_draft: only shallow ones become obstacles; the “deep” list is used by the visualiser to colour the safe-water background.
object_geoms is the projected object list; used by the visualiser to draw structures separately from ray-hit summaries.

Per-direction speed is computed by _weighted_avg_speed_knots(traffic_dir)() which sums freq * speed over every cell with non-zero frequency and divides by the total frequency.

`_run_all_computations()`¶

Source: geometries/get_powered_overlap.py:641 – _run_all_computations()

Loops every (seg_id, dir_idx) with non-zero speed and returns a list of computation dicts:

{
    'seg_id', 'leg', 'dir_idx', 'dir_info',
    'turn_pt', 'ext_dir', 'perp',
    'summaries': {(kind, obs_id): {'mass', 'mean_dist', 'p_integral',
                                    'p_approx', 'n_rays',
                                    'ray_offsets', 'ray_dists',
                                    'obs', 'kind'}},
    'ray_data', 'offsets', 'pdf_vals',
    'start', 'end',
    'cat1': {'origin', 'along_dir', 'length',
             'summaries': {(kind, obs_id): {...same keys...}},
             'ray_data'},
}

One entry per (leg, dir) that has obstacle hits in either category. summaries / ray_data are the Cat II result (rays from the turning point onwards); the cat1 block holds the Cat I result (rays along the leg itself).

Setup per (leg, dir):

u, n, L = _leg_vectors(start, end) – unit leg direction, unit perpendicular, leg length.
origin = start if dir_idx == 0 else end – the waypoint the ships come from; Cat I rays start here.
turn_pt = end if dir_idx == 0 else start – ships fail to turn at the downstream waypoint, so the Cat II ray origin is the turning point.
ext_dir = u for direction 0, -u for direction 1 (shared by both categories).

Then _compute_cat2_with_shadows(turn_pt, ext_dir, n, d['mean'], d['std'], d['ai'], d['speed_ms'], all_obstacles) and _compute_cat1_in_lane(origin, ext_dir, n, d['mean'], d['std'], L, all_obstacles) do the work.

`_compute_cat1_in_lane()`¶

Source: geometries/get_powered_overlap.py – _compute_cat1_in_lane()

The Category-I counterpart of the function below. It reuses _build_hit_matrix() and _accumulate_obs_hits() with two differences:

The rays start at the origin waypoint and max_range is the leg length, so only obstacles between the two waypoints are hit. Anything past the turning point is left to Cat II, which keeps the two categories disjoint.
_build_cat1_summaries() sets p_approx = p_integral = mass: there is no exp(-d / (a_i V)) factor because the ship is already on a collision course (Hansen eq. 4.15). mean_dist is still recorded (distance from the origin waypoint) for the visualiser.

In compute.powered_model, _iter_hit_probs() yields ('cat2', key, mass * exp(-d_mean / recovery)) for the Cat II summaries and ('cat1', key, mass) for the Cat I ones; the caller multiplies by Q and by the category’s causation factor (pc['grounding'] / pc['allision'] for Cat II, pc['grounding_cat1'] / pc['allision_cat1'] for Cat I) and the _ContribAccumulator keeps the split under by_category.

`_compute_cat2_with_shadows()`¶

Source: geometries/get_powered_overlap.py:415 – _compute_cat2_with_shadows()

This is the core of the powered model. It casts N_RAYS = 1000 parallel rays across the lateral distribution (mean +/- 6 sigma, compute/ray_sampling.py) and for each ray keeps the first obstacle it hits (that’s what produces the shadowing: closer obstacles block farther ones automatically).

The implementation is fully vectorised in numpy – no per-ray Python loop over shapely.

Steps:

Build the ray offsets: offsets = linspace(mean - 4 sigma, mean + 4 sigma, 500). masses[i] = pdf(offsets[i]) * dx is the mass of the lateral distribution attributable to that ray.
For each obstacle call _extract_edges_local(geom, turn_pt, ext_dir, perp)(). This transforms the polygon’s boundary into the local (along, lateral) frame where every ray is the horizontal line y = offset_i and returns an (M, 2, 2) array of edges.

For each obstacle compute a (n_rays, n_edges) matrix of along-track crossings in one numpy expression:

# Edge (y0,x0) -> (y1,x1); ray at y = ray_ys
crosses = (ray_ys >= y_min) & (ray_ys <= y_max) & (dy != 0)
t = (ray_ys - y0) / dy
along = x0 + t * (x1 - x0)
valid = crosses & (along > 0) & (along < MAX_RANGE)
along = where(valid, along, inf)
hit_matrix[:, obs_idx] = along.min(axis=1)   # nearest edge per ray

Per ray pick the first-hit obstacle: best_obs_idx = argmin(hit_matrix, axis=1), best_dists = hit_matrix[arange(n_rays), best_obs_idx]. Rays that miss every obstacle keep inf.
Accumulate per-obstacle stats (mass, weighted_dist, p_integral, ray_offsets, ray_dists) by walking the 500 rays once. Everything is O(N_RAYS) Python work because the heavy numerical work already happened.
Each obstacle summary:
- mean_dist = weighted_dist / mass
- p_approx = mass * _powered_na(mean_dist, ai, speed_ms) – the closed-form Cat II probability at the mean distance, equivalent for small dispersions.
- p_integral = sum_i m_i * exp(-dist_i / (ai * V)) – the properly-integrated Cat II probability (used for tighter bounds / visualisation).

This replaces an earlier ray-by-ray shapely loop that dominated end-to-end runtime; the vectorised form is ~74 x faster on the same 500-ray x many-obstacle problem.

`_extract_edges_local()`¶

Source: geometries/get_powered_overlap.py:225 – _extract_edges_local()

Helper used by _compute_cat2_with_shadows(). Walks the geometry (Polygon, MultiPolygon, LineString, MultiLineString, GeometryCollection) collecting the exterior + interior rings of every component. For every ring it subtracts turn_pt and projects onto along_dir and perp_dir (both unit vectors in world frame) using two matrix multiplications, producing an (M, 2, 2) edge array in the local frame.

Returns None for Point / MultiPoint geometries – a zero-width ray can never hit a point. Returns None for empty geometries.

Output: the traffic accumulation loop ¶

After _run_all_computations() returns, the mixin method walks traffic_data to apply per-cell frequency and causation factors. Each computation['summaries'] entry gives the mass fraction and mean distance for one obstacle in one (leg, direction); the traffic loop multiplies by the cell’s frequency q and adds P_c \cdot q \cdot mass \cdot \exp(-d/(a_i V)) to the running total.

Only two filters can exclude a contribution:

Grounding: the depth polygon’s depth must be \(\le\) ship’s draught. This is implicitly enforced by the depth-bin selection (max_draft passed to _build_legs_and_obstacles()).
Allision: ship_height >= obj_height. Explicit check in the traffic loop; lets short ships pass under tall structures (e.g., a tug under a high-clearance bridge).

The result is written directly to self.p.main_widget.LEPPoweredGrounding / LEPPoweredAllision and returned as a float.

Interactive visualiser ¶

The View button of the powered rows re-uses the same computations list to build an interactive matplotlib view. The same _run_all_computations() call – with the project’s Cat I start zone – produces the ray summaries; the visualiser walks them to draw the overview (selected leg, heading arrow and the outlines of the two local windows; legend and how-to beside it) on top, and Cat I (the leg in its own frame) | Cat II geometry over the Cat II probability mass x exp(-d / (ai V)) below. The local panels are true to scale: local_view sizes each window around the largest hits (the next ones and the ray start / turning point join while the window stays compact) and to the panel’s shape after a layout pass, and the lateral distribution is an inset at their left edge. The (leg, direction) is picked in the table only.

The table next to the figure comes from compute/powered_breakdown.py:build_sidebar_rows: one row per (leg, direction) with the route-table leg name, the annual frequency (sum of by_obstacle_leg[obs]["<seg>:<dir>"]) and its % of total. The frequencies come from <stem>.powered.json, written next to the run’s GeoPackage by RunHistoryMixin._write_powered_report_sidecar and handed to the dialog as calc.powered_view_reports. For a run saved before that sidecar existed, allision is rebuilt exactly from the View’s own ray casts (allision_report_from_computations, the same accumulation as run_powered_allision_model); grounding is ray cast once per draught bin in the model, so it is not rebuilt and its columns show --.

Visualiser: geometries/get_powered_overlap.py:869 – PoweredOverlapVisualizer

Function reference ¶

compute/powered_model.py (PoweredModelMixin)

run_powered_grounding_model – entry point, draught binning.
run_powered_allision_model – entry point, height filter.

geometries/get_powered_overlap.py

SimpleProjector – equirectangular lon/lat -> metres.
_parse_point – parse “x y” -> (float, float).
_project_wkt_geom – projects a shapely geometry via shapely.ops.transform().
_weighted_avg_speed_knots – per-cell freq-weighted speed.
_leg_vectors – leg unit vector, perpendicular, length.
_build_legs_and_obstacles – construct local-frame legs + obstacles.
_run_all_computations – per (leg, dir) ray-cast.
_compute_cat2_with_shadows – vectorised shadow-aware ray cast.
_extract_edges_local – polygon boundaries in the leg’s local frame.
_get_all_coords / _ray_hit_distance – kept for the interactive visualiser and tests (not on the hot path anymore).
_powered_na – \(\exp(-d/(a_i V))\) scalar helper.
local_view – true-scale window of a Cat I / Cat II panel, zoomed to the hits.
PoweredOverlapVisualizer – interactive matplotlib panel.

compute/powered_breakdown.py (QGIS-free)

build_sidebar_rows – (leg, direction) rows: name, frequency, % of total, per-ship hit probability.
leg_direction_frequencies / report_total – per leg-direction frequencies and the total from a model report or a sidecar block.
summarize_report / sidecar_payload – the .powered.json content.
allision_report_from_computations – rebuilds the allision report from the View’s ray casts for runs without a sidecar.

compute/basic_equations.py

powered_na – same formula as _powered_na, exposed as a module-level function for tests and the drifting cascade.

OMRAT

Navigation

Related Topics