.. _drifting: ============================ Drifting Risk Calculations ============================ This chapter explains, step by step, how OMRAT calculates the risk that a ship which loses propulsion (blackout) drifts into a shallow area (grounding), a structure (allision) or succeeds in anchoring before it does so. The chapter is built around the five worked examples in ``drifting/debug/level_1 ... level_5``. Each example runs OMRAT's actual calculation functions, prints every intermediate value, and produces a figure showing the leg, the obstacles and the drift corridor. Re-running the scripts reproduces the numbers quoted below exactly. .. contents:: In this chapter :local: :depth: 2 .. _drifting-overview: Overview of the Drifting Model =============================== When a ship suffers a blackout it becomes powerless and drifts under wind and current until one of the following happens: 1. The crew repairs the engine, 2. The crew successfully drops the anchor (in water shallow enough for anchoring), or 3. The ship reaches an obstacle (grounds on a shallow, collides with a structure or reaches the reach-distance limit without incident). The probability that a drifting ship causes an accident on an obstacle :math:`X` from leg :math:`i` in compass direction :math:`d` is .. math:: :label: eq-drift-per-direction P_{i,d,X} \;=\; \mathrm{base}_{i,k} \cdot r_{p,d} \cdot h_{X,\mathrm{eff}} \cdot P_{NR}(d_X) and for anchoring (which saves ships rather than causing an accident) .. math:: :label: eq-drift-anchor P_{\mathrm{anchor},i,d,X} \;=\; \mathrm{base}_{i,k}\cdot r_{p,d}\cdot a_p\cdot h_{\mathrm{anchor}} The total drifting accident rate is the sum over all wind-rose directions, all obstacles, all legs and all ship categories: .. math:: :label: eq-drift-total R_{\mathrm{drift}} \;=\; \sum_{i}\ \sum_{k}\ \sum_{X\in \{\mathrm{gr,al}\}}\ \sum_{d=1}^{8} r_{p,d}\cdot \mathrm{base}_{i,k}\cdot h_{X,\mathrm{eff},d}\cdot P_{NR,d}(d_X) with :math:`\sum_{d=1}^{8} r_{p,d} = 1`. For a uniform wind rose (used in the worked examples) :math:`r_{p,d} = 0.125`. Anchoring events are summed separately since they represent saved — not lost — ships. The base exposure is .. math:: :label: eq-base \mathrm{base}_{i,k} \;=\; \frac{L_i}{V_k\cdot 1852}\cdot f_{i,k}\cdot \lambda_{bo,\mathrm{hour}} \qquad \lambda_{bo,\mathrm{hour}} \;=\; \frac{\lambda_{bo}}{365.25\cdot 24} The symbols are listed in the table below. Sections :ref:`drifting-base`, :ref:`drifting-repair`, :ref:`drifting-hole`, :ref:`drifting-distance` and :ref:`drifting-run-vs-analysis` define each term in detail, and the five worked examples in Sections :ref:`drifting-level1` through :ref:`drifting-level5` show how they combine. .. list-table:: Symbols used throughout this chapter :header-rows: 1 :widths: 18 62 20 * - Symbol - Meaning - Unit / type * - :math:`L_i` - Length of leg :math:`i` - metres * - :math:`V_k` - Service speed of ship category :math:`k` - knots * - :math:`1852` - Metres per nautical mile - constant * - :math:`f_{i,k}` - Transits per year of category :math:`k` on leg :math:`i` - 1/year * - :math:`\lambda_{bo}` - Blackout rate while at sea - 1/year * - :math:`\lambda_{bo,\mathrm{hour}}` - Blackouts per hour = :math:`\lambda_{bo}/(365.25\cdot 24)` - 1/hour * - :math:`r_{p,d}` - Wind-rose probability for compass direction :math:`d` (8 directions, sum = 1) - dimensionless * - :math:`h_X` - Probability hole of obstacle :math:`X` (the fraction of the lateral traffic distribution whose drift rays intersect :math:`X`) - dimensionless * - :math:`h_{X,\mathrm{eff},d}` - Effective hole of obstacle :math:`X` in direction :math:`d` after upstream shadowing / anchoring - dimensionless * - :math:`d_X` - Along-drift distance from the leg reference to obstacle :math:`X` - metres * - :math:`V_{\mathrm{drift}}` - Drift speed - m/s * - :math:`P_{NR}(d_X)` - Probability the blackout is *not* repaired by the time the ship has drifted :math:`d_X` - dimensionless * - :math:`a_p` - Anchoring success probability - dimensionless * - :math:`a_d` - Anchoring depth multiplier (anchoring applies if depth is less than :math:`a_d\cdot T_{\mathrm{ship}}`) - dimensionless * - :math:`T_{\mathrm{ship}}` - Ship draught - metres * - :math:`S_X` - Quad-sweep shadow polygon of obstacle :math:`X` along the drift direction - polygon .. _drifting-run-vs-analysis: ``run analysis`` vs ``run model`` ================================== Two buttons in the OMRAT QGIS plugin trigger different code paths. .. list-table:: :header-rows: 1 :widths: 22 20 58 * - Button - Entry point - What it does * - **Run analysis** - :code:`run_drift_analysis()` in ``omrat.py:884`` → ``DriftCorridorTask`` - Builds the **drift corridor polygons** (per leg, per direction) and draws them as QGIS vector layers. Uses the quad-sweep shadow algorithm and the corridor-clipping code in ``geometries/drift/``. **No probabilities are computed.** * - **Run model** - :code:`run_calculation()` in ``omrat.py:599`` → ``CalculationTask`` → :code:`run_drifting_model()` in ``compute/drifting_model.py:2096`` - Computes the actual drifting risk — :math:`P(\mathrm{ground})`, :math:`P(\mathrm{allision})`, :math:`P(\mathrm{anchor})` — for every (leg, direction, obstacle, ship-category). This is the path described in the rest of this chapter. The quad-sweep algorithm described in :ref:`drifting-shadow-algorithm` is **only used by "run analysis"** (for drawing the corridors); the probability calculation uses the shadow-coverage formulation described per level below, which operates directly on the obstacle polygons and their quad-sweep shadows as a geometric subtraction rather than a corridor-clipping step. .. _drifting-parameters: Common Parameters Used in the Examples ======================================= All five worked examples use Leg 3 from the ``proj_3_3`` test scenario and an 8-vertex 12 m depth polygon (``BD5A4C46``) as the target. The values below are the exact inputs the scripts pass to the OMRAT functions — re-running the scripts reproduces every number shown in this chapter. .. list-table:: :header-rows: 1 :widths: 30 20 50 * - Parameter - Value - Symbol / note * - Leg 3 length - 34 113.2 m - :math:`L` * - Leg 3 start (lon, lat) - (14.24187°, 55.16728°) - EPSG:4326 * - Leg 3 end (lon, lat) - (14.59271°, 55.39937°) - EPSG:4326 * - Ship speed (service) - 12.5 kts - :math:`V` * - Ship frequency (Leg 3) - 610 transits/year - :math:`f` * - Ship draught (oil tanker) - 14.27 m - :math:`T_{\mathrm{ship}}` * - Drift speed - 1.94 kts (0.998 m/s) - :math:`V_{\mathrm{drift}}` * - Blackout rate - 1.0 /year - :math:`\lambda_{bo}` * - Blackouts per hour - :math:`1.1408\times 10^{-4}` - :math:`\lambda_{bo,\mathrm{hour}}` * - Wind rose (uniform) - 1/8 = 0.125 per direction - :math:`r_p` * - Lateral sigma - 500 m - :math:`\sigma` * - Reach distance - 50 000 m - max drift distance * - Repair time distribution - lognormal(std=1, loc=0, scale=1) - :math:`F_{\mathrm{lognorm}}` * - Anchoring success - :math:`a_p = 0.70` - * - Anchoring depth multiplier - :math:`a_d = 7.0` - threshold = :math:`a_d\cdot T_{\mathrm{ship}}` * - Drift direction (examples) - 315° (NW) - :math:`\theta_{\mathrm{compass}}` Compass-to-math conversion used throughout OMRAT: .. math:: \theta_{\mathrm{math}} = (90\degree - \theta_{\mathrm{compass}}) \bmod 360\degree For :math:`\theta_{\mathrm{compass}}=315\degree` (NW) this gives :math:`\theta_{\mathrm{math}}=135\degree` and a unit drift vector :math:`\hat v = (-0.7071, 0.7071)`. .. _drifting-base: Exposure Base ============== The base exposure counts how many blackouts per year are expected while ships are on the leg (equation :eq:`eq-base`). With the shared parameters: .. math:: \mathrm{hours\_present} = \frac{34113.2}{12.5\cdot 1852}\cdot 610 = 898.878\ \mathrm{h/yr} \qquad \mathrm{base} = 898.878\cdot 1.1408\times10^{-4} = 1.02541\times 10^{-1} The unit of :math:`\mathrm{base}` is *blackouts per year on this leg from ships of this category*. Every grounding, allision and anchoring probability below is multiplied by this factor and the wind-rose probability. .. container:: source-code-ref ``compute/basic_equations.py:9`` -- `get_drifting_prob() `__ .. _drifting-repair: Repair Time and :math:`P_{NR}` ============================== The time required to repair a blackout is modelled by a lognormal distribution with parameters :math:`(\sigma, \mu_{\mathrm{loc}}, s)`. The probability that the ship is *not* repaired by the time it has drifted a distance :math:`d` is .. math:: P_{NR}(d) \;=\; 1 - F_{\mathrm{lognorm}}\! \left( \frac{d}{V_{\mathrm{drift}}\cdot 3600};\ \sigma,\ \mu_{\mathrm{loc}},\ s \right) In the examples we use :math:`\sigma=1.0`, :math:`\mu_{\mathrm{loc}}=0`, :math:`s=1.0`. For the average target edge at :math:`d=11\,775\,\mathrm{m}` this evaluates to .. math:: \frac{d}{V_{\mathrm{drift}}\cdot 3600} = \frac{11775}{0.998\cdot 3600} = 3.28\ \mathrm{h} \qquad P_{NR}(11775) \approx 0.1176 i.e. about 11.8% of blackouts are still unresolved three hours after the event. .. container:: source-code-ref ``compute/basic_equations.py:30`` -- `get_not_repaired() `__ .. _drifting-hole: The Probability "Hole" :math:`h_X` =================================== The **hole** of an obstacle :math:`X` is the fraction of the lateral traffic distribution whose straight-line drift rays intersect :math:`X`. OMRAT computes it analytically by slicing the leg into N cross-sections and integrating the lateral PDF over the y-intervals where the drift ray hits the polygon (``geometries/analytical_probability.py``). For a lateral normal :math:`\mathcal N(0, \sigma^2)` the result is deterministic (no Monte Carlo noise): .. math:: h_X \;=\; \int_0^1 \!\!\!\!\int_{R(s)\subset[-W/2,+W/2]} \phi\!\left(\frac{y}{\sigma}\right) dy\, ds where :math:`R(s)` is the union of y-intervals at slice position :math:`s` where the drift ray hits :math:`X`, and :math:`W=5\sigma` is the lateral integration range. Numerical value in the examples: for the 8-vertex 12 m target polygon with drift direction NW, :math:`h_{\mathrm{target}} = 2.4915\times 10^{-2}`. .. container:: source-code-ref ``geometries/analytical_probability.py`` -- `compute_probability_analytical() `__. The analytical path is the default (``data['use_analytical']=True`` in ``run_drifting_model``); a Monte Carlo fallback exists in ``geometries/calculate_probability_holes.py`` but is only used when the user explicitly disables the analytical integrator. .. _drifting-distance: Directional Distance to an Obstacle ==================================== The distance :math:`d_X` that goes into :math:`P_{NR}(d_X)` is the *along-drift* distance from a leg reference line to the obstacle. OMRAT measures it by casting a **reverse ray** from each vertex of the obstacle back along the drift direction and intersecting with the reference line. The per-edge distance is the average of its two vertex distances (``edge_average_distance_m`` in ``drifting/engine.py``). Two reference lines are supported: - **Leg centreline** (``use_leg_offset = False``, default): distance is measured from the charted leg line itself. - **Mean-offset line** (``use_leg_offset = True``): distance is measured from a parallel line offset by :math:`\mathrm{mean\_offset\_m}` from the leg. This is the correct choice when the traffic in a direction is systematically biased to one side of the charted leg. For the target polygon (NW drift, leg centreline), the three front-facing edges give =========== ================= ================= =============== Edge vertex 0 (m) vertex 1 (m) average (m) =========== ================= ================= =============== 5 11 662.9 11 578.4 11 620.6 6 11 578.4 11 461.5 11 519.9 7 11 461.5 12 066.2 11 763.8 =========== ================= ================= =============== .. container:: source-code-ref ``drifting/engine.py:215`` -- `directional_distance_to_point_from_offset_leg() `__ | ``drifting/engine.py:288`` -- `edge_average_distance_m() `__ .. _drifting-level1: Level 1 — Single Polygon, Single Ship, Single Leg ================================================== **Script:** ``drifting/debug/level_1_single_polygon.py`` .. image:: ../../drifting/debug/level_1_single_polygon.png :width: 100% :alt: Level 1 drifting calculation The simplest case. Leg 3 drifts NW toward the 8-vertex 12 m polygon ``BD5A4C46``. One ship category (oil tanker 225–250 m, draught 14.27 m, 610 transits/year) contributes because the target depth (12 m) is less than the draught. Step 1 — base exposure (section :ref:`drifting-base`): .. math:: \mathrm{base} \;=\; 1.02541\times 10^{-1} Step 2 — analytical target hole (section :ref:`drifting-hole`): .. math:: h_{\mathrm{target}} \;=\; 2.4979\times 10^{-2} Step 3 — distribute the hole over the three front-facing edges in proportion to their length (total = 1215.0 m): =========== ============ ================= ===================== Edge idx length (m) distance (m) :math:`h_{\mathrm{edge}}` =========== ============ ================= ===================== 5 168.2 11 620.6 :math:`3.458\times 10^{-3}` 6 119.7 11 519.9 :math:`2.462\times 10^{-3}` 7 927.0 11 763.8 :math:`1.906\times 10^{-2}` =========== ============ ================= ===================== Step 4 — :math:`P_{NR}` per edge (section :ref:`drifting-repair`): .. math:: \begin{aligned} P_{NR}(11620.6)\text{ (Edge 5)} &= 0.12023 \\ P_{NR}(11519.9)\text{ (Edge 6)} &= 0.12198 \\ P_{NR}(11763.8)\text{ (Edge 7)} &= 0.11780 \end{aligned} Step 5 — per-edge grounding contribution: .. math:: P_{\mathrm{edge}} \;=\; \mathrm{base}\cdot r_p\cdot h_{\mathrm{edge}}\cdot P_{NR}(d_{\mathrm{edge}}) ================== ================================================ Edge Contribution (events/year) ================== ================================================ 5 :math:`5.329\times 10^{-6}` 6 :math:`3.849\times 10^{-6}` 7 :math:`2.878\times 10^{-5}` **Total** :math:`3.7955\times 10^{-5}` ================== ================================================ So **P(grounding) = 3.7955 × 10⁻⁵ events/year** for this single (leg, ship, direction) combination. .. _drifting-level1b: Level 1b — Distance from the Leg Centreline vs the Distribution Centre ======================================================================= **Script:** ``drifting/debug/level_1_single_polygon.py`` (second figure) .. image:: ../../drifting/debug/level_1b_distribution_center.png :width: 100% :alt: Level 1b — two mean-offset reference lines (one per traffic direction) A *charted* leg is a single line, but the **traffic** on that leg is typically split into two directions (ships going one way and ships going the other way) and those two populations may be centred on slightly different parallel lines. OMRAT supports this via two reference-line modes for the :math:`d_X` measurement (:ref:`drifting-distance`): - **Mode (a) — ``leg_center`` (default):** distance is measured from the charted leg line. Used everywhere in Levels 1–5 unless stated. - **Mode (b) — ``distribution_center``:** distance is measured from a mean-offset line, one per traffic direction. This shifts the distance values by the mean lateral offset of the traffic for that direction. The same front-facing target edges computed with both modes (:math:`\mathrm{mean\_offset}_A = +500\,\mathrm{m}`, :math:`\mathrm{mean\_offset}_B = -500\,\mathrm{m}`): .. list-table:: Edge distances per reference line (metres) :header-rows: 1 :widths: 10 22 22 22 22 * - Edge - Leg centreline (v0 / v1 / avg) - Offset A = +500 (v0 / v1 / avg) - Offset B = -500 (v0 / v1 / avg) - Shift * - 5 - 11 662.9 / 11 578.4 / 11 620.6 - 11 161.8 / 11 077.4 / 11 119.6 - 12 163.9 / 12 079.5 / 12 121.7 - ±500 m * - 6 - 11 578.4 / 11 461.5 / 11 519.9 - 11 077.4 / 10 960.4 / 11 018.9 - 12 079.5 / 12 562.5 / 12 320.9 - ±500 m * - 7 - 11 461.5 / 12 066.2 / 11 763.8 - 10 960.4 / 11 565.1 / 11 262.8 - 11 961.5 / 12 566.8 / 12 264.0 - ±500 m Using the mean-offset distance changes :math:`P_{NR}(d_X)` (shorter distance → ship more likely to still be drifting). Combined over the two traffic directions this produces a more realistic risk than treating both directions as centred on the charted leg. To enable this mode set ``use_leg_offset_for_distance = True`` in the :code:`DriftConfig` and provide a non-zero :code:`LegState.mean_offset_m` per direction. .. container:: source-code-ref ``drifting/engine.py:215`` -- `directional_distance_to_point_from_offset_leg(use_leg_offset=True) `__ | ``drifting/engine.py:46`` -- `DriftConfig.use_leg_offset_for_distance `__ .. _drifting-level2: Level 2 — Two Legs, Multiple Ship Categories ============================================= **Script:** ``drifting/debug/level_2_two_legs_multi_ships.py`` .. image:: ../../drifting/debug/level_2_two_legs_multi_ships.png :width: 100% :alt: Level 2 — two legs, multiple ship categories This example introduces: 1. A second traffic leg (Leg 6, 7 489 m, bearing ~41° NE) with its own drift direction (N = 0°). 2. Five ship categories with different draughts and frequencies. Only ship categories with :math:`T_{\mathrm{ship}} > d_{\mathrm{polygon}} = 12\,\mathrm{m}` can ground on this polygon. Three of the five qualify: ========================== ============= ========= ========== ========== ======= Ship category Draught (m) Grounds? Freq L3 Freq L6 Speed ========================== ============= ========= ========== ========== ======= Oil tanker 225-250 m 14.27 yes 610 50 12.5 General cargo 225-250 m 11.82 no 450 40 13.0 Bulk carrier 250-275 m 16.53 yes 180 15 13.5 Container 275-300 m 13.50 yes 95 8 18.0 Passenger 100-125 m 5.80 no 320 280 16.0 ========================== ============= ========= ========== ========== ======= Per-leg analytical holes (NW for Leg 3, N for Leg 6): ========== ===================================================== Leg :math:`h_X` (NW for L3, N for L6) ========== ===================================================== LEG_3 :math:`2.4979\times 10^{-2}` LEG_6 :math:`8.1012\times 10^{-2}` ========== ===================================================== The per-(leg, ship) exposure base scales with that ship's speed and frequency: .. math:: \mathrm{base}_{i,k} = \frac{L_i}{V_k\cdot 1852}\cdot f_{i,k}\cdot \lambda_{bo,\mathrm{hour}} Summing over all valid (leg, ship, edge) combinations and the single NW direction yields .. list-table:: :header-rows: 1 :widths: 40 20 20 20 * - Leg / ship - Base - P(ground) - Share * - LEG_3 / Oil tanker 225-250 m - 0.10254 - :math:`3.7955\times 10^{-5}` - 70.8% * - LEG_3 / Bulk carrier 250-275 m - 0.02802 - :math:`1.0370\times 10^{-5}` - 19.4% * - LEG_3 / Container 275-300 m - 0.01109 - :math:`4.1049\times 10^{-6}` - 7.7% * - LEG_6 / Oil tanker 225-250 m - 0.001845 - :math:`8.281\times 10^{-7}` - 1.5% * - LEG_6 / Bulk carrier 250-275 m - 0.000513 - :math:`2.300\times 10^{-7}` - 0.4% * - LEG_6 / Container 275-300 m - 0.000205 - :math:`9.201\times 10^{-8}` - 0.2% * - **Total (NW)** - - :math:`\mathbf{5.358\times 10^{-5}}` - Multiplying by the rose probability and summing over all 8 wind directions gives the full per-leg risk used in equation :eq:`eq-drift-total`. .. _drifting-level3: Level 3 — Blocking Polygon (Shadow Coverage) ============================================= **Script:** ``drifting/debug/level_3_blocking_polygon.py`` .. image:: ../../drifting/debug/level_3_blocking_polygon.png :width: 100% :alt: Level 3 — three shadow-coverage scenarios When an obstacle sits upstream of another obstacle on the drift axis, drifting ships that would otherwise reach the downstream obstacle may be blocked by the upstream one. In OMRAT the blocking effect is **purely geometric**: a blocker removes exactly the drift rays it physically intercepts, nothing more. Shadow-Coverage Model --------------------- Define the blocker's shadow polygon :math:`S_B` as the quad-sweep of the blocker along the drift direction (see :ref:`drifting-shadow-algorithm`). The effective hole of the downstream target becomes .. math:: :label: eq-shadow-cover h_{\mathrm{target,eff}} \;=\; h\!\big(\mathrm{target} - S_B\big) which is computed by the same analytical integrator used for any other hole, but evaluated on the *unshadowed part* of the target polygon. Equivalently, .. math:: h_{\mathrm{target,eff}} \;=\; h_{\mathrm{target}}\cdot\big(1 - \mathrm{cov}_B\big) \qquad \mathrm{cov}_B \;=\; \frac{h_{\mathrm{target}} - h_{\mathrm{target,eff}}} {h_{\mathrm{target}}} where :math:`\mathrm{cov}_B` is the *probability-hole coverage* (not the geometric area fraction — the lateral PDF weights the target unevenly). The grounding probability is .. math:: P(\mathrm{target}) \;=\; \mathrm{base}\cdot r_p\cdot h_{\mathrm{target,eff}} \cdot P_{NR}(d_{\mathrm{target}}) .. important:: There is **no** multiplicative "remaining" factor in the blocker cascade. The blocker only changes the target's effective hole via equation :eq:`eq-shadow-cover`. Three Scenarios --------------- Using the same target polygon as Level 1/2 and a rectangular blocker placed 4 km upstream of the target, the script generates three coverage scenarios by varying the blocker's cross-drift width. For each, it reports the effective hole, the blocker's own hole and the target probability. .. list-table:: Level 3 — three shadow coverages :header-rows: 1 :widths: 30 15 25 30 * - Scenario - cov(hole) - :math:`h_{\mathrm{target,eff}}` - P(target) * - A — blocker laterally offset - 0.0% - :math:`2.4915\times 10^{-2}` - :math:`3.7559\times 10^{-5}` * - B — 30 % coverage - 30.0% - :math:`1.7441\times 10^{-2}` - :math:`2.6291\times 10^{-5}` * - C — 100 % coverage - 100.0% - 0 - 0 * - Baseline (no blocker) - - :math:`2.4915\times 10^{-2}` - :math:`3.7559\times 10^{-5}` Detailed Calculation (Scenario B, 30% coverage) ----------------------------------------------- .. math:: \begin{aligned} h_{\mathrm{target}} &= 2.4915\times 10^{-2} \\ \mathrm{cov}_B\text{ (hole)} &= 30.00\% \\ h_{\mathrm{target,eff}} &= h_{\mathrm{target}}\cdot(1-0.30) = 1.7441\times 10^{-2} \\ \mathrm{base} &= 1.02541\times 10^{-1} \\ r_p &= 0.125 \\ P_{NR}(d_{\mathrm{target}}=11\,775) &= 0.11761 \\ P(\mathrm{target}) &= 1.02541\times 10^{-1}\cdot 0.125 \cdot 1.7441\times 10^{-2}\cdot 0.11761 \\ &= 2.6291\times 10^{-5}\ \text{events/year} \end{aligned} The blocker itself may also be a grounding or allision hazard; the script reports its contribution separately using its own hole, distance and (for grounding/allision) :math:`P_{NR}`. .. _drifting-level4: Level 4 — Anchoring ==================== **Script:** ``drifting/debug/level_4_anchoring.py`` .. image:: ../../drifting/debug/level_4_anchoring.png :width: 100% :alt: Level 4 — three anchor coverage scenarios An *anchoring polygon* is an area where ships can successfully drop anchor because the water depth is shallower than :math:`a_d\cdot T_{\mathrm{ship}}` but deeper than the ship's draught (so it is not a grounding hazard). Anchoring is modelled *probabilistically*: a ship drifting through the anchor zone anchors with probability :math:`a_p` and keeps drifting with probability :math:`(1-a_p)`. The anchoring formulation uses the same shadow-coverage geometry as Level 3, combined with the probabilistic success factor. If :math:`S_A` is the anchor polygon's quad-sweep shadow along the drift direction, the effective target hole is .. math:: h_{\mathrm{target,eff}} \;=\; h_{\mathrm{target}} - a_p\cdot h\!\big(\mathrm{target}\cap S_A\big) The grounding and anchoring contributions are .. math:: P(\mathrm{target}) = \mathrm{base}\cdot r_p \cdot h_{\mathrm{target,eff}} \cdot P_{NR}(d_{\mathrm{target}}) \qquad P(\mathrm{anchor}) = \mathrm{base}\cdot r_p\cdot a_p\cdot h_{\mathrm{anchor}} When :math:`a_p \to 1`, the anchor behaves like a pure blocker (the Level 3 full-coverage case). When :math:`a_p = 0` it has no effect. Parameters specific to this example: - :math:`a_p = 0.70`, :math:`a_d = 7.0`, :math:`T_{\mathrm{ship}} = 14.27\,\mathrm{m}` - Anchoring threshold: :math:`a_d\cdot T_{\mathrm{ship}} = 99.9\,\mathrm{m}` - Anchor polygon depth: 50 m (< 99.9 m, so anchoring applies) Three Scenarios --------------- The script varies the anchor polygon's cross-drift width to produce three coverage cases. .. list-table:: Level 4 — anchor shadow coverages :header-rows: 1 :widths: 30 15 20 20 20 * - Scenario - cov(hole) - :math:`h_{\mathrm{target,eff}}` - P(anchor) - P(target) * - A — no anchor polygon - 0.0% - :math:`2.4915\times 10^{-2}` - 0 - :math:`3.7559\times 10^{-5}` * - B — 30% coverage - 30.0% - :math:`1.9479\times 10^{-2}` - :math:`9.353\times 10^{-5}` - :math:`2.9672\times 10^{-5}` * - C — 100% coverage - 100.0% - :math:`7.475\times 10^{-3}` - :math:`3.283\times 10^{-4}` - :math:`1.1268\times 10^{-5}` In Scenario C the anchor reduces the target hole by exactly :math:`a_p = 70\%`, giving :math:`P(\mathrm{target})` equal to 30% of the baseline. In Scenario B it reduces it by :math:`0.3\cdot 0.7 = 21\%`, consistent with the reported -21%. Detailed Calculation (Scenario C, 100% coverage) ------------------------------------------------ .. math:: \begin{aligned} h_{\mathrm{anchor}} &= 3.6585\times 10^{-2} \\ h(\mathrm{target}\cap S_A) &= h_{\mathrm{target}} = 2.4915\times 10^{-2} \\ h_{\mathrm{target,eff}} &= h_{\mathrm{target}} - a_p\cdot h(\mathrm{target}\cap S_A) \\ &= 2.4915\times 10^{-2} - 0.7\cdot 2.4915\times 10^{-2} \\ &= 7.475\times 10^{-3} \\ P(\mathrm{anchor}) &= \mathrm{base}\cdot r_p\cdot a_p \cdot h_{\mathrm{anchor}} \\ &= 1.02541\times 10^{-1}\cdot 0.125 \cdot 0.7\cdot 3.6585\times 10^{-2} \\ &= 3.283\times 10^{-4}\ \text{/yr} \\ P(\mathrm{target}) &= \mathrm{base}\cdot r_p \cdot h_{\mathrm{target,eff}} \cdot P_{NR}(d_{\mathrm{target}}) \\ &= 1.02541\times 10^{-1}\cdot 0.125 \cdot 7.475\times 10^{-3} \cdot 0.11761 \\ &= 1.1268\times 10^{-5}\ \text{/yr} \end{aligned} .. _drifting-level5: Level 5 — Complete Cascade =========================== **Script:** ``drifting/debug/level_5_complete_cascade.py`` .. image:: ../../drifting/debug/level_5_complete_cascade.png :width: 100% :alt: Level 5 — anchoring + allision + grounding The full cascade combines all three obstacle types in the same drift path: 1. **Anchoring polygon** (50 m depth, 8 km upstream, wider than target): probabilistic, reduces subsequent rays by :math:`a_p = 0.70` over the cross-drift intersection. 2. **Allision structure** (platform / turbine cluster, 600 m × 200 m, 4 km upstream of target, directly in front on the cross-drift axis): pure geometric blocker. 3. **Grounding target** (same 8-vertex 12 m polygon). Combined Shadow-Coverage Formulation ------------------------------------- Let :math:`S_s` be the structure's quad-sweep shadow and :math:`S_a` be the anchor's. Only the part of the target *not* shadowed by the structure is reachable: .. math:: \mathrm{target}_{\mathrm{reach}} \;=\; \mathrm{target} - S_s Of that reachable part, some rays also pass through the anchor zone and are reduced by :math:`a_p`: .. math:: h_{\mathrm{target,eff}} \;=\; h\!\big(\mathrm{target}_{\mathrm{reach}}\big) - a_p\cdot h\!\big(\mathrm{target}_{\mathrm{reach}} \cap S_a\big) The structure itself is reduced by the anchor for its own rays that pass through the anchor zone first: .. math:: h_{\mathrm{struct,eff}} \;=\; h_{\mathrm{struct}} - a_p\cdot h\!\big(\mathrm{struct}\cap S_a\big) Anchoring is unconditional — every ray passing through the anchor zone is counted, regardless of whether a downstream structure would later have blocked it: .. math:: P(\mathrm{anchor}) \;=\; \mathrm{base}\cdot r_p\cdot a_p\cdot h_{\mathrm{anchor}} Probability contributions: .. math:: \begin{aligned} P(\mathrm{allision}) &= \mathrm{base}\cdot r_p \cdot h_{\mathrm{struct,eff}} \cdot P_{NR}(d_{\mathrm{struct}}) \\ P(\mathrm{ground}) &= \mathrm{base}\cdot r_p \cdot h_{\mathrm{target,eff}} \cdot P_{NR}(d_{\mathrm{target}}) \end{aligned} The total accident rate is :math:`P(\mathrm{allision}) + P(\mathrm{ground})`. :math:`P(\mathrm{anchor})` is shown separately since those ships are saved. Numerical Walk-Through ----------------------- Values computed by the script (see ``drifting/debug/level_5_complete_cascade.py``): .. list-table:: Obstacle holes and distances :header-rows: 1 :widths: 40 30 30 * - Obstacle - :math:`h_X` - :math:`d_X` (m) * - Anchor (50 m, 8 km upstream) - :math:`4.2542\times 10^{-2}` - 3 734 * - Structure (600 m × 200 m, 4 km upstream of target) - :math:`1.7640\times 10^{-2}` - 7 734 * - Target (12 m) - :math:`2.4915\times 10^{-2}` - 11 775 Shadow coverages (computed by the script): .. list-table:: Shadow coverages of downstream obstacles :header-rows: 1 :widths: 50 50 * - Coverage - Value * - Structure shadow on target hole - 70.8% * - Anchor shadow on REACHABLE target hole - 100.0% * - Anchor shadow on structure hole - 100.0% Effective holes: .. math:: \begin{aligned} h_{\mathrm{target,eff}} &= h(\mathrm{target}-S_s) - a_p\cdot h((\mathrm{target}-S_s)\cap S_a) \\ &= 7.277\times 10^{-3} - 0.7\cdot 7.277\times 10^{-3} \\ &= 2.183\times 10^{-3} \\ h_{\mathrm{struct,eff}} &= h_{\mathrm{struct}} - a_p\cdot h(\mathrm{struct}\cap S_a) \\ &= 1.7640\times 10^{-2} - 0.7\cdot 1.7640\times 10^{-2} \\ &= 5.292\times 10^{-3} \end{aligned} Probability contributions (using :math:`P_{NR}(7734)=0.22164` and :math:`P_{NR}(11775)=0.11761`): .. math:: \begin{aligned} P(\mathrm{anchor}) &= 1.02541\times 10^{-1}\cdot 0.125\cdot 0.70\cdot 4.2542\times 10^{-2} = 3.817\times 10^{-4} \\ P(\mathrm{allision}) &= 1.02541\times 10^{-1}\cdot 0.125\cdot 5.292\times 10^{-3}\cdot 0.22164 \\ &= 1.503\times 10^{-5} \\ P(\mathrm{ground}) &= 1.02541\times 10^{-1}\cdot 0.125\cdot 2.183\times 10^{-3}\cdot 0.11761 \\ &= 3.290\times 10^{-6} \\ \text{Accident rate} &= P(\mathrm{allision}) + P(\mathrm{ground}) = 1.832\times 10^{-5}\ \text{/yr} \end{aligned} Scenario Comparison ------------------- Dropping one obstacle at a time shows how each component affects the total: .. list-table:: :header-rows: 1 :widths: 35 20 20 20 * - Scenario - P(anchor) - P(allision) - P(ground) * - target only - 0 - 0 - :math:`3.756\times 10^{-5}` * - target + anchor - :math:`3.817\times 10^{-4}` - 0 - :math:`1.127\times 10^{-5}` * - target + structure - 0 - :math:`5.011\times 10^{-5}` - :math:`1.097\times 10^{-5}` * - target + anchor + structure - :math:`3.817\times 10^{-4}` - :math:`1.503\times 10^{-5}` - :math:`3.290\times 10^{-6}` The structure blocks 70.8% of the target's rays (the part directly behind it); the anchor reduces the remaining target rays (and the structure's rays) by :math:`a_p = 0.7`. .. _drifting-shadow-algorithm: Shadow Algorithm (Quad-Sweep) ============================== .. note:: The quad-sweep algorithm described in this section is the same geometric operation used everywhere else in the chapter (for example, :math:`S_B` in :eq:`eq-shadow-cover` and :math:`S_s`, :math:`S_a` in Level 5). It is also the only place a shadow polygon is explicitly **drawn**: in the ``run analysis`` code path that produces the QGIS drift-corridor layers. During the probability calculation (``run model``) the shadow polygons are built on the fly per obstacle and used directly as a :code:`shapely.Polygon.difference` / :code:`Polygon.intersection` argument — the result is the effective-hole equation already shown per level. There is no corridor-clipping or MultiPolygon reachability filter on the probability path. OMRAT computes shadow polygons by sweeping each obstacle along the drift direction and joining the original polygon to its translated copy via a quadrilateral per edge. .. image:: _static/images/shadow_algorithm.svg :width: 100% :alt: Quad-sweep shadow algorithm Given polygon :math:`P = (v_0, v_1, \ldots, v_{n-1})` and drift vector :math:`\hat v \cdot L`: 1. Translate :math:`P` by :math:`\hat v \cdot L` to produce :math:`P'`. 2. For each edge :math:`(v_i, v_{i+1})` create a quadrilateral :math:`Q_i = (v_i, v_{i+1}, v'_{i+1}, v'_i)`. 3. The shadow is :math:`S = P \cup P' \cup \bigcup_i Q_i`. This preserves the original contour of the obstacle even for **concave** shapes (e.g. a reef with a navigable gap) — the quads fill only the region physically swept by each edge. A convex-hull shortcut would incorrectly fill navigable channels. Reproducing this operation manually in the worked examples: .. code-block:: python def build_quad_shadow(poly, extrude_length): far = translate(poly, xoff=drift_ux * extrude_length, yoff=drift_uy * extrude_length) orig = list(poly.exterior.coords)[:-1] far_c = list(far.exterior.coords)[:-1] quads = [] for i in range(len(orig)): j = (i + 1) % len(orig) q = Polygon([orig[i], orig[j], far_c[j], far_c[i]]) if q.is_valid and q.area > 0: quads.append(q) return unary_union([poly, far] + quads) .. container:: source-code-ref ``geometries/drift/shadow.py:18`` -- `create_obstacle_shadow() `__ *(visualisation path only; used by ``run analysis`` to produce the corridor layers)* .. _drifting-corridor-viz: Drift Corridor Visualisation (``run analysis`` only) ====================================================== The **drift corridor** layers drawn on the map come from a separate code path (``run analysis``) and are not part of the probability calculation. They are a convenience for the user to see, for each leg and each of the 8 compass directions, the reachable area after all obstacles have been clipped out. Base surface width: .. math:: \text{width} = 2\cdot Z_{0.995}\cdot \sigma \approx 2\cdot 2.576\cdot \sigma Maximum drift distance: .. math:: d_{\mathrm{max}} = F^{-1}_{\mathrm{lognorm}}(0.999) \cdot 3600\cdot V_{\mathrm{drift}} clamped to :math:`[10\,000, 50\,000]\,\mathrm{m}` for numerical stability. The visualisation corridor is the convex hull of the base surface and its translation by :math:`d_{\mathrm{max}}\hat v`, clipped by every obstacle's quad-sweep shadow (:ref:`drifting-shadow-algorithm`) and filtered for reachability. **This clipping happens only on the visualisation path** — during ``run model`` the probability calculation uses the shadow polygons directly (via shadow-coverage subtraction, per-level equations above). .. container:: source-code-ref ``omrat.py:884`` -- `run_drift_analysis() `__ | ``geometries/drift/corridor.py:16`` -- `create_base_surface() `__ | ``geometries/drift/clipping.py:16`` -- `clip_corridor_at_obstacles() `__ .. _drifting-pipeline: Pipeline and Source Code Pointers ================================== .. container:: source-code-ref pipeline **Probability pipeline (``run model``):** ``omrat.py:599 run_calculation()`` → ``CalculationTask`` → ``compute/drifting_model.py:2096 run_drifting_model()`` → `analytical hole `__ + `edge distance `__ + `P_NR `__ + shadow-coverage subtraction per Levels 3–5. **Visualisation pipeline (``run analysis``):** ``omrat.py:884 run_drift_analysis()`` → ``DriftCorridorTask`` → `DriftCorridorGenerator `__ → `quad-sweep shadow `__ + `corridor clipping `__. The worked-example scripts in ``drifting/debug/`` import the same functions OMRAT uses during ``run model``, so the numbers printed above match the results the QGIS plugin produces for the same inputs. Summary ======== - For every (leg, direction, obstacle, ship) tuple the drifting probability is :math:`P = \mathrm{base}\cdot r_p\cdot h_{\mathrm{eff}}\cdot P_{NR}(d)`, and anchoring is :math:`P_{\mathrm{anchor}} = \mathrm{base}\cdot r_p\cdot a_p\cdot h_{\mathrm{anchor}}`. - :math:`h_{\mathrm{eff}}` is the analytical hole of the obstacle *minus* everything that is geometrically or probabilistically blocked upstream. Blocking polygons (Level 3) remove the rays they intercept; anchoring zones (Levels 4–5) reduce those rays by the factor :math:`a_p`. No multiplicative "remaining" factor is applied across obstacles — the correct formulation is the shadow-coverage model, identical in form across blockers (:math:`h - h(\cap S)`) and anchoring (:math:`h - a_p\cdot h(\cap S_A)`). - The final drifting rate is the sum over all 8 wind-rose directions, all obstacles, all legs and all ship categories (equation :eq:`eq-drift-total`). - Levels 1–5 of the worked examples in ``drifting/debug/`` reproduce every intermediate value in this chapter. - The quad-sweep shadow algorithm is used in both the probability calculation (via polygon subtraction) and the drift-corridor visualisation; the corridor-clipping step is **visualisation-only**.