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

This chapter walks OMRAT’s IWRAP 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 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:28 – run_powered_grounding_model()

Allision entry: compute/powered_model.py:204 – 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.
`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:51 – 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:28 – 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:204 – 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:322 – _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:424 – _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',
}

One entry per (leg, dir) that has obstacle hits.

Setup per (leg, dir):

u, n, L = _leg_vectors(start, end) – unit leg direction, unit perpendicular, leg length.
turn_pt = end if dir_idx == 0 else start – ships fail to turn at the downstream waypoint, so the ray origin is the turning point.
ext_dir = u for direction 0, -u for direction 1.

Then _compute_cat2_with_shadows(turn_pt, ext_dir, n, d['mean'], d['std'], d['ai'], d['speed_ms'], all_obstacles) does the work.

`_compute_cat2_with_shadows()`¶

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

This is the core of the powered model. It casts N_RAYS = 500 parallel rays across the lateral distribution 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:121 – _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 ¶

For “Show analysis” (previously called “Run analysis”) the plugin re-uses the same computations list to build an interactive matplotlib view. The same _run_all_computations() call produces the ray summaries; the visualiser walks them to draw the overview map, detail panel, and waterfall breakdown.

Visualiser: geometries/get_powered_overlap.py:473 – 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.
find_closest_computation_index – maps a click to a computation (used by the visualiser).
PoweredOverlapVisualizer – interactive matplotlib panel.

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