Files
venus/dbus-anchor-alarm/analysis/live_track_debug.py
2026-03-26 14:15:02 +00:00

1209 lines
46 KiB
Python

#!/usr/bin/env python3
"""
Live anchor-estimation debug tool.
Connects to a Venus OS device (cerbo) via MQTT, pulls current anchor alarm
state (GPS track, wind, chain, estimations), and generates an HTML report
showing which track segments are reliable for anchor position estimation.
Reliability is based on two criteria:
1. Wind speed >= threshold (default 7 kts) — chain is under tension
2. Vessel distance from marked anchor >= ratio * chain_length (default 0.60)
— the rode is taut enough to point toward the anchor
Usage:
python3 live_track_debug.py --host cerbo
python3 live_track_debug.py --host 192.168.1.100 --wind-threshold 8 --tension-ratio 0.55
"""
import argparse
import json
import math
import os
import sys
import threading
import time
sys.path.insert(0, os.path.join(os.path.dirname(__file__), ".."))
from catenary import catenary_distance, wind_force_lbs
try:
import paho.mqtt.client as mqtt
except ImportError:
sys.exit("paho-mqtt is required: pip install paho-mqtt")
EARTH_RADIUS_FT = 20_902_231.0
NM_TO_FT = 6076.12
PATHS = [
"/Anchor/Set",
"/Anchor/Marked/Latitude",
"/Anchor/Marked/Longitude",
"/Anchor/Estimated/Latitude",
"/Anchor/Estimated/Longitude",
"/Anchor/UncertaintyRadius",
"/Anchor/Drift",
"/Anchor/EstimatedHistory/Json",
"/Anchor/ChainLength",
"/Anchor/EstimatedDistance",
"/Vessel/Latitude",
"/Vessel/Longitude",
"/Vessel/Speed",
"/Vessel/Course",
"/Vessel/Heading",
"/Vessel/Depth",
"/Vessel/WindSpeed",
"/Vessel/WindDirection",
"/Track/Json",
"/Track/PointCount",
"/Settings/ChainLength",
"/Settings/FreeboardHeight",
"/Settings/WindageArea",
"/Settings/DragCoefficient",
"/Settings/ChainWeight",
"/Settings/Anchor/Latitude",
"/Settings/Anchor/Longitude",
]
# ---------------------------------------------------------------------------
# Geo helpers (same as anchor_tracker.py)
# ---------------------------------------------------------------------------
def haversine_ft(lat1, lon1, lat2, lon2):
rlat1, rlat2 = math.radians(lat1), math.radians(lat2)
dlat = math.radians(lat2 - lat1)
dlon = math.radians(lon2 - lon1)
a = (math.sin(dlat / 2) ** 2
+ math.cos(rlat1) * math.cos(rlat2) * math.sin(dlon / 2) ** 2)
return EARTH_RADIUS_FT * 2.0 * math.atan2(math.sqrt(a), math.sqrt(1.0 - a))
def project_position(lat, lon, distance_ft, bearing_deg):
d = distance_ft / EARTH_RADIUS_FT
brng = math.radians(bearing_deg)
rlat = math.radians(lat)
rlon = math.radians(lon)
new_lat = math.asin(
math.sin(rlat) * math.cos(d)
+ math.cos(rlat) * math.sin(d) * math.cos(brng)
)
new_lon = rlon + math.atan2(
math.sin(brng) * math.sin(d) * math.cos(rlat),
math.cos(d) - math.sin(rlat) * math.sin(new_lat),
)
return math.degrees(new_lat), math.degrees(new_lon)
def to_local_xy(lat, lon, ref_lat, ref_lon):
lat_rad = math.radians(ref_lat)
dx = (lon - ref_lon) * 60.0 * NM_TO_FT * math.cos(lat_rad)
dy = (lat - ref_lat) * 60.0 * NM_TO_FT
return dx, dy
def from_local_xy(x_ft, y_ft, ref_lat, ref_lon):
lat_rad = math.radians(ref_lat)
cos_lat = math.cos(lat_rad)
if abs(cos_lat) < 1e-12:
cos_lat = 1e-12
dlon = x_ft / (60.0 * NM_TO_FT * cos_lat)
dlat = y_ft / (60.0 * NM_TO_FT)
return ref_lat + dlat, ref_lon + dlon
def circle_fit(points_xy):
"""Algebraic circle fit (Kasa method). Returns (cx, cy, radius) or None."""
n = len(points_xy)
if n < 3:
return None
sx = sy = sxx = syy = sxy = 0.0
sxxx = syyy = sxyy = syxx = 0.0
for x, y in points_xy:
x2, y2 = x * x, y * y
sx += x; sy += y
sxx += x2; syy += y2; sxy += x * y
sxxx += x2 * x; syyy += y2 * y
sxyy += x * y2; syxx += y * x2
a11 = 2.0 * (sx * sx - n * sxx)
a12 = 2.0 * (sx * sy - n * sxy)
a22 = 2.0 * (sy * sy - n * syy)
b1 = (sx * sxx - n * sxxx) + (sx * syy - n * sxyy)
b2 = (sy * sxx - n * syxx) + (sy * syy - n * syyy)
det = a11 * a22 - a12 * a12
if abs(det) < 1e-12:
return None
cx = (b1 * a22 - b2 * a12) / det
cy = (a11 * b2 - a12 * b1) / det
radius_sq = (sxx + syy - 2.0 * cx * sx - 2.0 * cy * sy) / n + cx * cx + cy * cy
if radius_sq < 0:
return None
return cx, cy, math.sqrt(radius_sq)
# ---------------------------------------------------------------------------
# Estimation algorithms
# ---------------------------------------------------------------------------
def run_estimators(track, chain_length_ft, tension_ratio, ref_lat, ref_lon):
"""Run competing anchor estimation algorithms on track data.
Each algorithm produces an (x, y) estimate in local coordinates relative
to ref_lat/ref_lon, converted back to lat/lon. Returns a list of dicts.
The track has no per-point wind data, so algorithms use distance-from-
centroid as a proxy for chain tension (farther = more wind = more taut).
"""
if len(track) < 3:
return []
pts_ll = [(p["lat"], p["lon"]) for p in track]
pts_xy = [to_local_xy(la, lo, ref_lat, ref_lon) for la, lo in pts_ll]
results = []
# ---- A. Simple centroid ----
cx = sum(x for x, _ in pts_xy) / len(pts_xy)
cy = sum(y for _, y in pts_xy) / len(pts_xy)
results.append(_make_result(
"centroid", "Simple Centroid",
"Unweighted average of all track positions.",
cx, cy, ref_lat, ref_lon,
))
# ---- B. Distance-weighted centroid (iterative) ----
wcx, wcy = cx, cy
for _ in range(5):
total_w = 0.0
wx = wy = 0.0
for x, y in pts_xy:
d = math.sqrt((x - wcx) ** 2 + (y - wcy) ** 2) + 1.0
wx += x * d
wy += y * d
total_w += d
if total_w > 0:
wcx = wx / total_w
wcy = wy / total_w
results.append(_make_result(
"dist_weighted", "Distance-Weighted Centroid",
"Points farther from center get higher weight (proxy for chain "
"tension / wind). Iterates 5x.",
wcx, wcy, ref_lat, ref_lon,
))
# ---- C. Distance-squared-weighted centroid ----
wcx2, wcy2 = cx, cy
for _ in range(5):
total_w = 0.0
wx = wy = 0.0
for x, y in pts_xy:
d2 = (x - wcx2) ** 2 + (y - wcy2) ** 2 + 1.0
wx += x * d2
wy += y * d2
total_w += d2
if total_w > 0:
wcx2 = wx / total_w
wcy2 = wy / total_w
results.append(_make_result(
"dist2_weighted", "Distance²-Weighted Centroid",
"Squared distance weighting — strongly favors outermost points "
"where chain is taut.",
wcx2, wcy2, ref_lat, ref_lon,
))
# ---- D. Circle fit (all points) ----
fit = circle_fit(pts_xy)
if fit:
fcx, fcy, frad = fit
r = _make_result(
"circle_fit", "Circle Fit (all points)",
f"Kasa algebraic circle fit on all {len(pts_xy)} points. "
f"Fitted radius: {frad:.0f} ft.",
fcx, fcy, ref_lat, ref_lon,
)
r["fit_radius"] = round(frad, 1)
results.append(r)
# ---- E. Tension-gated centroid ----
dists_from_centroid = [
math.sqrt((x - cx) ** 2 + (y - cy) ** 2) for x, y in pts_xy
]
max_dist = max(dists_from_centroid) if dists_from_centroid else 1.0
gate = tension_ratio * max_dist
gated = [(x, y) for (x, y), d in zip(pts_xy, dists_from_centroid) if d >= gate]
if len(gated) >= 3:
gcx = sum(x for x, _ in gated) / len(gated)
gcy = sum(y for _, y in gated) / len(gated)
results.append(_make_result(
"tension_gated", f"Tension-Gated Centroid (>{tension_ratio:.0%})",
f"Only uses {len(gated)}/{len(pts_xy)} points beyond "
f"{tension_ratio:.0%} of max swing radius ({gate:.0f} ft). "
f"These points had taut chain.",
gcx, gcy, ref_lat, ref_lon,
))
# ---- F. Boundary circle fit ----
if len(pts_xy) >= 10:
sorted_by_dist = sorted(
zip(pts_xy, dists_from_centroid), key=lambda t: t[1], reverse=True,
)
boundary_n = max(5, len(pts_xy) // 4)
boundary_pts = [p for p, _ in sorted_by_dist[:boundary_n]]
bfit = circle_fit(boundary_pts)
if bfit:
bcx, bcy, brad = bfit
r = _make_result(
"boundary_fit", f"Boundary Circle Fit (top 25%)",
f"Circle fit on outermost {boundary_n} points only. "
f"These had the highest chain tension. "
f"Fitted radius: {brad:.0f} ft.",
bcx, bcy, ref_lat, ref_lon,
)
r["fit_radius"] = round(brad, 1)
results.append(r)
# ---- G. Tension-gated circle fit ----
if len(gated) >= 5:
tfit = circle_fit(gated)
if tfit:
tcx, tcy, trad = tfit
r = _make_result(
"tension_circle", f"Tension-Gated Circle Fit (>{tension_ratio:.0%})",
f"Circle fit on {len(gated)} points beyond tension threshold. "
f"Fitted radius: {trad:.0f} ft.",
tcx, tcy, ref_lat, ref_lon,
)
r["fit_radius"] = round(trad, 1)
results.append(r)
return results
def _make_result(algo_id, name, description, x, y, ref_lat, ref_lon):
lat, lon = from_local_xy(x, y, ref_lat, ref_lon)
error_ft = round(haversine_ft(ref_lat, ref_lon, lat, lon), 1)
return {
"id": algo_id,
"name": name,
"description": description,
"x": round(x, 1),
"y": round(y, 1),
"lat": round(lat, 7),
"lon": round(lon, 7),
"error_ft": error_ft,
}
# ---------------------------------------------------------------------------
# MQTT data collector
# ---------------------------------------------------------------------------
class AnchorDataCollector:
"""Connects to Venus MQTT and collects anchor alarm state."""
def __init__(self, host, port=1883):
self.host = host
self.port = port
self.portal_id = None
self.values = {}
self._ready = threading.Event()
self._got_portal = threading.Event()
self._received_count = 0
def collect(self, timeout=10):
client = mqtt.Client(client_id="anchor-debug", protocol=mqtt.MQTTv311)
client.on_connect = self._on_connect
client.on_message = self._on_message
print(f"Connecting to {self.host}:{self.port} ...")
client.connect(self.host, self.port, keepalive=60)
client.loop_start()
if not self._got_portal.wait(timeout=timeout):
client.loop_stop()
client.disconnect()
sys.exit(f"Timed out waiting for portal ID from {self.host}")
print(f"Portal ID: {self.portal_id}")
prefix = f"N/{self.portal_id}/anchoralarm/0"
client.subscribe(f"{prefix}/#")
time.sleep(0.5)
for path in PATHS:
topic = f"R/{self.portal_id}/anchoralarm/0{path}"
client.publish(topic, "")
deadline = time.time() + timeout
while time.time() < deadline:
time.sleep(0.5)
if self._received_count >= len(PATHS) - 2:
time.sleep(1.0)
break
client.loop_stop()
client.disconnect()
print(f"Collected {self._received_count} values")
return self.values
def _on_connect(self, client, _userdata, _flags, rc):
if rc != 0:
sys.exit(f"MQTT connection failed: rc={rc}")
client.subscribe("N/+/anchoralarm/0/#")
def _on_message(self, _client, _userdata, msg):
parts = msg.topic.split("/")
if len(parts) < 4:
return
if self.portal_id is None:
self.portal_id = parts[1]
self._got_portal.set()
aa_idx = None
for i, p in enumerate(parts):
if p == "anchoralarm":
aa_idx = i
break
if aa_idx is None or aa_idx + 1 >= len(parts):
return
path = "/" + "/".join(parts[aa_idx + 2:])
try:
payload = json.loads(msg.payload.decode())
value = payload.get("value", payload)
except (json.JSONDecodeError, UnicodeDecodeError):
return
self.values[path] = value
self._received_count += 1
# ---------------------------------------------------------------------------
# Analysis
# ---------------------------------------------------------------------------
def analyze(values, wind_threshold, tension_ratio,
known_anchor_lat=None, known_anchor_lon=None):
"""Run reliability analysis on collected data.
If known_anchor_lat/lon are provided, they are used as the true anchor
position for all distance and XY calculations. The service's marked
and estimated positions are still reported for comparison.
"""
track_raw = values.get("/Track/Json", "[]")
if isinstance(track_raw, str):
try:
track = json.loads(track_raw)
except json.JSONDecodeError:
track = []
else:
track = track_raw if isinstance(track_raw, list) else []
marked_lat = _float(values.get("/Anchor/Marked/Latitude"))
marked_lon = _float(values.get("/Anchor/Marked/Longitude"))
est_lat = _float(values.get("/Anchor/Estimated/Latitude"))
est_lon = _float(values.get("/Anchor/Estimated/Longitude"))
chain_length = _float(values.get("/Anchor/ChainLength")) or _float(values.get("/Settings/ChainLength")) or 150.0
est_distance = _float(values.get("/Anchor/EstimatedDistance")) or 0.0
vessel_wind = _float(values.get("/Vessel/WindSpeed")) or 0.0
vessel_wind_dir = _float(values.get("/Vessel/WindDirection"))
vessel_heading = _float(values.get("/Vessel/Heading"))
vessel_depth = _float(values.get("/Vessel/Depth")) or 0.0
freeboard = _float(values.get("/Settings/FreeboardHeight")) or 4.0
windage_area = _float(values.get("/Settings/WindageArea")) or 200.0
drag_coeff = _float(values.get("/Settings/DragCoefficient")) or 1.1
chain_weight = _float(values.get("/Settings/ChainWeight")) or 2.25
uncertainty = _float(values.get("/Anchor/UncertaintyRadius")) or 0.0
drift = _float(values.get("/Anchor/Drift")) or 0.0
anchor_set = values.get("/Anchor/Set", 0)
est_history_raw = values.get("/Anchor/EstimatedHistory/Json", "[]")
if isinstance(est_history_raw, str):
try:
est_history = json.loads(est_history_raw)
except json.JSONDecodeError:
est_history = []
else:
est_history = est_history_raw if isinstance(est_history_raw, list) else []
if not track:
print("WARNING: No track data received. Is the anchor set?")
using_known = known_anchor_lat is not None and known_anchor_lon is not None
if using_known:
ref_lat = known_anchor_lat
ref_lon = known_anchor_lon
print(f"Using known anchor position: {ref_lat:.6f}, {ref_lon:.6f}")
elif marked_lat is not None and marked_lon is not None:
ref_lat = marked_lat
ref_lon = marked_lon
else:
print("WARNING: No marked anchor position. Using estimated position as reference.")
if est_lat is not None and est_lon is not None:
ref_lat, ref_lon = est_lat, est_lon
elif track:
ref_lat, ref_lon = track[0]["lat"], track[0]["lon"]
else:
ref_lat, ref_lon = 0.0, 0.0
distance_threshold = tension_ratio * chain_length
classified = []
for pt in track:
lat, lon = pt["lat"], pt["lon"]
ts = pt.get("ts", 0)
dist = haversine_ft(ref_lat, ref_lon, lat, lon)
ratio = dist / chain_length if chain_length > 0 else 0.0
wind_ok = vessel_wind >= wind_threshold
tension_ok = dist >= distance_threshold
if wind_ok and tension_ok:
reliability = "reliable"
elif wind_ok or tension_ok:
reliability = "marginal"
else:
reliability = "unreliable"
x, y = to_local_xy(lat, lon, ref_lat, ref_lon)
classified.append({
"ts": ts, "lat": lat, "lon": lon,
"x": round(x, 1), "y": round(y, 1),
"dist": round(dist, 1), "ratio": round(ratio, 3),
"reliability": reliability,
})
est_history_xy = []
for eh in est_history:
x, y = to_local_xy(eh["lat"], eh["lon"], ref_lat, ref_lon)
est_history_xy.append({
"x": round(x, 1), "y": round(y, 1),
"ts": eh.get("ts", 0),
"uncertainty_ft": eh.get("uncertainty_ft", 0),
})
n_reliable = sum(1 for c in classified if c["reliability"] == "reliable")
n_marginal = sum(1 for c in classified if c["reliability"] == "marginal")
n_unreliable = sum(1 for c in classified if c["reliability"] == "unreliable")
n_total = len(classified)
pct_reliable = round(100.0 * n_reliable / n_total, 1) if n_total else 0.0
est_x, est_y = (0.0, 0.0)
if est_lat is not None and est_lon is not None:
est_x, est_y = to_local_xy(est_lat, est_lon, ref_lat, ref_lon)
force = wind_force_lbs(vessel_wind, windage_area, drag_coeff)
cat = catenary_distance(chain_length, vessel_depth, freeboard, force, chain_weight)
estimators = run_estimators(track, chain_length, tension_ratio, ref_lat, ref_lon)
marked_x, marked_y = (0.0, 0.0)
if marked_lat is not None and marked_lon is not None:
marked_x, marked_y = to_local_xy(marked_lat, marked_lon, ref_lat, ref_lon)
est_error_ft = None
marked_error_ft = None
if using_known:
if est_lat is not None and est_lon is not None:
est_error_ft = round(haversine_ft(ref_lat, ref_lon, est_lat, est_lon), 1)
if marked_lat is not None and marked_lon is not None:
marked_error_ft = round(haversine_ft(ref_lat, ref_lon, marked_lat, marked_lon), 1)
return {
"track": classified,
"estimators": estimators,
"est_history": est_history_xy,
"ref": {"lat": ref_lat, "lon": ref_lon, "using_known": using_known},
"marked": {
"lat": marked_lat, "lon": marked_lon,
"x": round(marked_x, 1), "y": round(marked_y, 1),
"error_ft": marked_error_ft,
},
"estimated": {
"lat": est_lat, "lon": est_lon,
"x": round(est_x, 1), "y": round(est_y, 1),
"error_ft": est_error_ft,
},
"params": {
"chain_length": chain_length,
"wind_threshold": wind_threshold,
"tension_ratio": tension_ratio,
"distance_threshold": round(distance_threshold, 1),
"vessel_wind": round(vessel_wind, 1),
"vessel_wind_dir": vessel_wind_dir,
"vessel_heading": vessel_heading,
"vessel_depth": round(vessel_depth, 1),
"freeboard": freeboard,
"windage_area": windage_area,
"drag_coeff": drag_coeff,
"chain_weight": chain_weight,
"est_distance": round(est_distance, 1),
"catenary_distance": round(cat.total_distance_ft, 1),
"chain_on_bottom": round(cat.chain_on_bottom_ft, 1),
"uncertainty": round(uncertainty, 1),
"drift": round(drift, 1),
"anchor_set": bool(anchor_set),
},
"stats": {
"total_points": n_total,
"reliable": n_reliable,
"marginal": n_marginal,
"unreliable": n_unreliable,
"pct_reliable": pct_reliable,
},
}
def _float(v):
if v is None:
return None
try:
return float(v)
except (TypeError, ValueError):
return None
# ---------------------------------------------------------------------------
# HTML report
# ---------------------------------------------------------------------------
HTML_TEMPLATE = r"""<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Anchor Estimation Debug</title>
<script src="https://cdn.jsdelivr.net/npm/chart.js@4"></script>
<style>
:root {
--bg: #0f1117; --surface: #1a1d27; --border: #2a2d3a;
--text: #e0e0e8; --muted: #8a8d9a; --accent: #4fc3f7;
--green: #66bb6a; --yellow: #fdd835; --red: #ef5350;
--accent2: #f4845f; --accent3: #7ec8a0; --accent4: #c39af4;
}
* { box-sizing: border-box; margin: 0; padding: 0; }
body {
background: var(--bg); color: var(--text);
font-family: 'SF Mono', 'Fira Code', monospace;
padding: 2rem; line-height: 1.6;
}
h1 { color: var(--accent); font-size: 1.6rem; margin-bottom: 0.3rem; }
.subtitle { color: var(--muted); margin-bottom: 1.5rem; font-size: 0.85rem; }
h2 { color: var(--accent); font-size: 1.1rem; margin: 2rem 0 0.75rem;
border-bottom: 1px solid var(--border); padding-bottom: 0.3rem; }
.stats-grid {
display: grid; grid-template-columns: repeat(auto-fit, minmax(180px, 1fr));
gap: 0.75rem; margin: 1rem 0;
}
.stat-card {
background: var(--surface); border: 1px solid var(--border);
border-radius: 6px; padding: 0.75rem;
}
.stat-card .label { color: var(--muted); font-size: 0.75rem; text-transform: uppercase; }
.stat-card .value { font-size: 1.3rem; font-weight: 700; }
.val-accent { color: var(--accent); }
.val-green { color: var(--green); }
.val-yellow { color: var(--yellow); }
.val-red { color: var(--red); }
.chart-row { display: grid; grid-template-columns: 1fr 1fr; gap: 1.5rem; margin: 1rem 0; }
.chart-box {
background: var(--surface); border: 1px solid var(--border);
border-radius: 8px; padding: 1rem;
}
.chart-box.full { grid-column: 1 / -1; }
.chart-box canvas { width: 100% !important; }
.insight {
color: var(--muted); font-size: 0.85rem; margin-top: 0.5rem; font-style: italic;
background: var(--surface); border-left: 3px solid var(--accent);
padding: 0.75rem 1rem; border-radius: 0 6px 6px 0;
}
.param-table { width: 100%; border-collapse: collapse; margin: 1rem 0; }
.param-table td { padding: 0.4rem 0.75rem; border-bottom: 1px solid var(--border); }
.param-table td:first-child { color: var(--muted); font-size: 0.85rem; width: 45%; }
.param-table td:last-child { color: var(--accent); font-weight: 600; }
@media (max-width: 900px) { .chart-row { grid-template-columns: 1fr; } }
</style>
</head>
<body>
<h1>Anchor Estimation Debug Report</h1>
<p class="subtitle">Generated TIMESTAMP &bull; Host: HOST &bull; Wind threshold: WIND_THR kts &bull; Tension ratio: TENSION_RATIO</p>
<h2>Current State</h2>
<div class="stats-grid">
<div class="stat-card"><div class="label">Anchor Set</div><div class="value val-accent">ANCHOR_SET</div></div>
<div class="stat-card"><div class="label">Chain Length</div><div class="value val-accent">CHAIN_LEN ft</div></div>
<div class="stat-card"><div class="label">Current Wind</div><div class="value WIND_CLASS">WIND_NOW kts</div></div>
<div class="stat-card"><div class="label">Est. Distance</div><div class="value val-accent">EST_DIST ft</div></div>
<div class="stat-card"><div class="label">Catenary Distance</div><div class="value val-accent">CAT_DIST ft</div></div>
<div class="stat-card"><div class="label">Chain on Bottom</div><div class="value val-accent">CHAIN_BOTTOM ft</div></div>
<div class="stat-card"><div class="label">Uncertainty</div><div class="value val-accent">UNCERTAINTY ft</div></div>
<div class="stat-card"><div class="label">Drift from Marked</div><div class="value val-accent">DRIFT ft</div></div>
</div>
KNOWN_ANCHOR_SECTION
<h2>Reliability Analysis</h2>
<div class="stats-grid">
<div class="stat-card"><div class="label">Track Points</div><div class="value val-accent">N_TOTAL</div></div>
<div class="stat-card"><div class="label">Reliable</div><div class="value val-green">N_RELIABLE (PCT_REL%)</div></div>
<div class="stat-card"><div class="label">Marginal</div><div class="value val-yellow">N_MARGINAL</div></div>
<div class="stat-card"><div class="label">Unreliable</div><div class="value val-red">N_UNRELIABLE</div></div>
<div class="stat-card"><div class="label">Distance Threshold</div><div class="value val-accent">DIST_THR ft</div></div>
<div class="stat-card"><div class="label">Depth</div><div class="value val-accent">DEPTH ft</div></div>
</div>
<p class="insight">
Reliable = wind &ge; WIND_THR kts AND distance from marked anchor &ge; TENSION_RATIO &times; chain length (DIST_THR ft).
In light winds the chain coils on the seabed in a spiral; the vessel drifts randomly and heading does not point toward the anchor.
Only when the rode is taut (&asymp;60% of chain length) under steady wind can the vessel heading reliably indicate anchor bearing.
</p>
<!-- 1. Track map -->
<h2>1. Track Map &mdash; Reliability Overlay</h2>
<div class="chart-row">
<div class="chart-box full"><canvas id="trackChart" height="450"></canvas></div>
</div>
<p class="insight">Vessel positions in local X/Y (ft) relative to ANCHOR_REF_LABEL. Green = reliable, yellow = marginal, red = unreliable. Dashed circles show the distance threshold and chain length radius.</p>
<!-- 2. Distance + wind timeline -->
<h2>2. Distance &amp; Wind vs Time</h2>
<div class="chart-row">
<div class="chart-box full"><canvas id="distWindChart" height="300"></canvas></div>
</div>
<p class="insight">Distance from marked anchor (blue) and chain-length distance threshold (dashed). Only points above the threshold with sufficient wind are considered reliable.</p>
<!-- 3. Reliability breakdown bar -->
<h2>3. Distance Ratio Distribution</h2>
<div class="chart-row">
<div class="chart-box"><canvas id="ratioHist" height="250"></canvas></div>
<div class="chart-box"><canvas id="reliabilityPie" height="250"></canvas></div>
</div>
<p class="insight">Left: histogram of distance-from-anchor / chain-length ratios. The vertical line marks the tension threshold. Right: overall reliability breakdown.</p>
<!-- 4. Estimation history -->
<h2>4. Estimated Anchor Position History</h2>
<div class="chart-row">
<div class="chart-box full"><canvas id="estHistoryChart" height="350"></canvas></div>
</div>
<p class="insight">History of estimated anchor positions (X/Y). Spread shows estimation uncertainty over time. Tight clusters indicate confidence; scattered points indicate the estimator is chasing noise.</p>
<!-- 5. Algorithm comparison -->
<h2>5. Estimation Algorithm Comparison</h2>
<div class="chart-row">
<div class="chart-box full"><canvas id="algoChart" height="400"></canvas></div>
</div>
ALGO_TABLE
<p class="insight">
Each algorithm estimates the anchor position using only the GPS track geometry.
Distance from center acts as a proxy for chain tension (farther = more wind = taut rode).
Algorithms that gate or weight by distance effectively filter for wind conditions, even without
per-point wind data. Lower error = closer to ALGO_REF_LABEL.
</p>
<!-- 6. Parameters -->
<h2>6. System Parameters</h2>
<div class="chart-row">
<div class="chart-box">
<table class="param-table">
<tr><td>Chain length</td><td>CHAIN_LEN ft</td></tr>
<tr><td>Vessel depth</td><td>DEPTH ft</td></tr>
<tr><td>Freeboard</td><td>FREEBOARD ft</td></tr>
<tr><td>Windage area</td><td>WINDAGE sqft</td></tr>
<tr><td>Drag coefficient</td><td>DRAG_COEFF</td></tr>
<tr><td>Chain weight</td><td>CHAIN_WT lb/ft</td></tr>
</table>
</div>
<div class="chart-box">
<table class="param-table">
<tr><td>Wind threshold</td><td>WIND_THR kts</td></tr>
<tr><td>Tension ratio</td><td>TENSION_RATIO</td></tr>
<tr><td>Distance threshold</td><td>DIST_THR ft</td></tr>
<tr><td>Current wind</td><td>WIND_NOW kts</td></tr>
<tr><td>Wind direction</td><td>WIND_DIR&deg;</td></tr>
<tr><td>Vessel heading</td><td>HEADING&deg;</td></tr>
</table>
</div>
</div>
<script>
const D = DATA_JSON;
const commonOpts = {
responsive: true,
plugins: { legend: { labels: { color: '#8a8d9a' } } },
scales: {
x: { ticks: { color: '#8a8d9a' }, grid: { color: '#2a2d3a' } },
y: { ticks: { color: '#8a8d9a' }, grid: { color: '#2a2d3a' } },
},
};
const COLORS = { reliable: '#66bb6a', marginal: '#fdd835', unreliable: '#ef5350' };
// --- 1. Track scatter with reliability colors ---
(function() {
const pts = D.track;
const data = pts.map(p => ({ x: p.x, y: p.y }));
const colors = pts.map(p => COLORS[p.reliability]);
const datasets = [{
label: 'Vessel track',
data, backgroundColor: colors,
pointRadius: 2.5, showLine: false,
}];
// Distance threshold circle
const thrRadius = D.params.distance_threshold;
const circlePts = (r, label, color, dash) => {
const pts = [];
for (let a = 0; a <= 360; a += 2) {
const rad = a * Math.PI / 180;
pts.push({ x: r * Math.cos(rad), y: r * Math.sin(rad) });
}
return {
label, data: pts, borderColor: color, borderDash: dash,
pointRadius: 0, showLine: true, fill: false, borderWidth: 1.5,
};
};
datasets.push(circlePts(thrRadius, `Tension threshold (${thrRadius} ft)`, '#fdd835', [6, 3]));
datasets.push(circlePts(D.params.chain_length, `Chain length (${D.params.chain_length} ft)`, '#8a8d9a', [3, 3]));
// Markers
if (D.ref.using_known) {
// Known anchor at origin — large filled marker
datasets.push({
label: 'Known anchor (0,0)', data: [{ x: 0, y: 0 }],
backgroundColor: '#00e676', borderColor: '#00e676',
pointRadius: 12, pointStyle: 'triangle', showLine: false,
});
// Crosshairs at origin for visibility
datasets.push({
label: '_knownH', data: [{ x: -20, y: 0 }, { x: 20, y: 0 }],
borderColor: '#00e676', pointRadius: 0, showLine: true, fill: false,
borderWidth: 2, borderDash: [2, 2],
});
datasets.push({
label: '_knownV', data: [{ x: 0, y: -20 }, { x: 0, y: 20 }],
borderColor: '#00e676', pointRadius: 0, showLine: true, fill: false,
borderWidth: 2, borderDash: [2, 2],
});
if (D.marked.x !== undefined && D.marked.x !== null) {
datasets.push({
label: `Marked anchor (${D.marked.error_ft || '?'} ft error)`,
data: [{ x: D.marked.x, y: D.marked.y }],
backgroundColor: '#ff5252', borderColor: '#ff5252',
pointRadius: 8, pointStyle: 'crossRot', showLine: false,
});
}
} else {
datasets.push({
label: 'Marked anchor (0,0)', data: [{ x: 0, y: 0 }],
backgroundColor: '#ff5252', borderColor: '#ff5252',
pointRadius: 10, pointStyle: 'triangle', showLine: false,
});
}
if (D.estimated.x !== null) {
const errLabel = D.estimated.error_ft !== null ? ` (${D.estimated.error_ft} ft error)` : '';
datasets.push({
label: 'Service estimate' + errLabel, data: [{ x: D.estimated.x, y: D.estimated.y }],
backgroundColor: '#4fc3f7', borderColor: '#4fc3f7',
pointRadius: 9, pointStyle: 'star', showLine: false,
});
// Uncertainty circle around service estimate
const unc = D.params.uncertainty;
if (unc > 0) {
const uncPts = [];
for (let a = 0; a <= 360; a += 3) {
const rad = a * Math.PI / 180;
uncPts.push({ x: D.estimated.x + unc * Math.cos(rad), y: D.estimated.y + unc * Math.sin(rad) });
}
datasets.push({
label: `Uncertainty (${unc} ft)`,
data: uncPts, borderColor: 'rgba(79,195,247,0.5)', borderDash: [4, 4],
pointRadius: 0, showLine: true, fill: false, borderWidth: 1.5,
});
}
}
new Chart('trackChart', {
type: 'scatter', data: { datasets },
options: {
...commonOpts, aspectRatio: 1.1,
plugins: {
legend: {
labels: {
color: '#8a8d9a',
filter: (item) => !item.text.startsWith('_'),
},
},
},
scales: {
x: { title: { display: true, text: 'East-West (ft)', color: '#8a8d9a' }, ticks: { color: '#8a8d9a' }, grid: { color: '#2a2d3a' } },
y: { title: { display: true, text: 'North-South (ft)', color: '#8a8d9a' }, ticks: { color: '#8a8d9a' }, grid: { color: '#2a2d3a' } },
},
},
});
})();
// --- 2. Distance + wind timeline ---
(function() {
const pts = D.track;
if (!pts.length) return;
const t0 = pts[0].ts;
const hours = pts.map(p => ((p.ts - t0) / 3600).toFixed(2));
new Chart('distWindChart', {
type: 'scatter',
data: {
datasets: [
{
label: 'Distance from marked (ft)',
data: pts.map((p, i) => ({ x: parseFloat(hours[i]), y: p.dist })),
backgroundColor: pts.map(p => COLORS[p.reliability]),
pointRadius: 2, showLine: false, yAxisID: 'y',
},
{
label: `Threshold (${D.params.distance_threshold} ft)`,
data: [{ x: 0, y: D.params.distance_threshold }, { x: parseFloat(hours[hours.length-1]) || 1, y: D.params.distance_threshold }],
borderColor: '#fdd835', borderDash: [6, 3],
pointRadius: 0, showLine: true, fill: false, borderWidth: 1.5, yAxisID: 'y',
},
],
},
options: {
...commonOpts,
scales: {
x: { title: { display: true, text: 'Hours', color: '#8a8d9a' }, ticks: { color: '#8a8d9a' }, grid: { color: '#2a2d3a' } },
y: { position: 'left', title: { display: true, text: 'Distance (ft)', color: '#8a8d9a' }, ticks: { color: '#8a8d9a' }, grid: { color: '#2a2d3a' } },
},
},
});
})();
// --- 3a. Ratio histogram ---
(function() {
const ratios = D.track.map(p => p.ratio);
if (!ratios.length) return;
const bins = 20;
const maxR = Math.min(Math.max(...ratios) * 1.1, 2.0);
const step = maxR / bins;
const counts = new Array(bins).fill(0);
const labels = [];
for (let i = 0; i < bins; i++) {
labels.push((i * step).toFixed(2));
ratios.forEach(r => { if (r >= i * step && r < (i+1) * step) counts[i]++; });
}
const thrBin = Math.floor(D.params.tension_ratio / step);
const bgColors = counts.map((_, i) => i >= thrBin ? 'rgba(102,187,106,0.6)' : 'rgba(239,83,80,0.6)');
new Chart('ratioHist', {
type: 'bar',
data: { labels, datasets: [{ label: 'Distance / Chain Length', data: counts, backgroundColor: bgColors, borderWidth: 0 }] },
options: {
...commonOpts,
scales: {
x: { title: { display: true, text: 'Distance Ratio', color: '#8a8d9a' }, ticks: { color: '#8a8d9a' }, grid: { color: '#2a2d3a' } },
y: { title: { display: true, text: 'Count', color: '#8a8d9a' }, ticks: { color: '#8a8d9a' }, grid: { color: '#2a2d3a' } },
},
},
});
})();
// --- 3b. Reliability pie ---
(function() {
new Chart('reliabilityPie', {
type: 'doughnut',
data: {
labels: ['Reliable', 'Marginal', 'Unreliable'],
datasets: [{
data: [D.stats.reliable, D.stats.marginal, D.stats.unreliable],
backgroundColor: ['#66bb6a', '#fdd835', '#ef5350'],
borderWidth: 0,
}],
},
options: {
responsive: true,
plugins: {
legend: { labels: { color: '#8a8d9a' } },
},
},
});
})();
// --- 4. Estimation history scatter ---
(function() {
const eh = D.est_history;
if (!eh.length) return;
const data = eh.map(p => ({ x: p.x, y: p.y }));
const unc = eh.map(p => Math.max(2, Math.min(p.uncertainty_ft / 10, 10)));
const datasets = [{
label: 'Estimated anchor',
data, backgroundColor: 'rgba(195,154,244,0.5)',
pointRadius: unc, showLine: false,
}];
const refLabel = D.ref.using_known ? 'Known anchor (0,0)' : 'Marked anchor (0,0)';
const refColor = D.ref.using_known ? '#00e676' : '#ff5252';
datasets.push({
label: refLabel, data: [{ x: 0, y: 0 }],
backgroundColor: refColor, pointRadius: 8, pointStyle: 'crossRot', showLine: false,
});
new Chart('estHistoryChart', {
type: 'scatter', data: { datasets },
options: {
...commonOpts, aspectRatio: 1.2,
scales: {
x: { title: { display: true, text: 'East-West (ft)', color: '#8a8d9a' }, ticks: { color: '#8a8d9a' }, grid: { color: '#2a2d3a' } },
y: { title: { display: true, text: 'North-South (ft)', color: '#8a8d9a' }, ticks: { color: '#8a8d9a' }, grid: { color: '#2a2d3a' } },
},
},
});
})();
// --- 5. Algorithm comparison scatter ---
(function() {
const algos = D.estimators;
if (!algos || !algos.length) return;
const ALGO_COLORS = [
'#4fc3f7', '#f4845f', '#c39af4', '#66bb6a', '#fdd835', '#e040fb', '#26c6da'
];
const trackPts = D.track;
const datasets = [];
if (trackPts.length) {
datasets.push({
label: 'Vessel track',
data: trackPts.map(p => ({ x: p.x, y: p.y })),
backgroundColor: 'rgba(255,255,255,0.08)',
pointRadius: 1.2, showLine: false,
});
}
const refLabel2 = D.ref.using_known ? 'Known anchor' : 'Marked anchor';
const refColor2 = D.ref.using_known ? '#00e676' : '#ff5252';
datasets.push({
label: refLabel2 + ' (0,0)', data: [{ x: 0, y: 0 }],
backgroundColor: refColor2, pointRadius: 10, pointStyle: 'crossRot', showLine: false,
});
if (D.estimated.x !== null) {
datasets.push({
label: `Service estimate (${D.estimated.error_ft ?? '?'} ft)`,
data: [{ x: D.estimated.x, y: D.estimated.y }],
backgroundColor: '#ffffff', pointRadius: 7, pointStyle: 'star', showLine: false,
});
}
algos.forEach((a, i) => {
const c = ALGO_COLORS[i % ALGO_COLORS.length];
datasets.push({
label: `${a.name} (${a.error_ft} ft)`,
data: [{ x: a.x, y: a.y }],
backgroundColor: c, pointRadius: 8,
pointStyle: 'rectRot', showLine: false,
});
});
new Chart('algoChart', {
type: 'scatter', data: { datasets },
options: {
...commonOpts, aspectRatio: 1.1,
scales: {
x: { title: { display: true, text: 'East-West (ft)', color: '#8a8d9a' }, ticks: { color: '#8a8d9a' }, grid: { color: '#2a2d3a' } },
y: { title: { display: true, text: 'North-South (ft)', color: '#8a8d9a' }, ticks: { color: '#8a8d9a' }, grid: { color: '#2a2d3a' } },
},
},
});
})();
</script>
</body>
</html>"""
def render_report(results, host, wind_threshold, tension_ratio):
html = HTML_TEMPLATE
p = results["params"]
s = results["stats"]
html = html.replace("DATA_JSON", json.dumps(results))
html = html.replace("TIMESTAMP", time.strftime("%Y-%m-%d %H:%M:%S"))
html = html.replace("HOST", host)
ref = results["ref"]
mk = results["marked"]
est = results["estimated"]
if ref["using_known"]:
known_section = (
'<h2>Known Anchor vs Service</h2>'
'<div class="stats-grid">'
f' <div class="stat-card"><div class="label">Known Anchor</div>'
f' <div class="value val-green" style="font-size:0.95rem">{ref["lat"]:.6f}, {ref["lon"]:.6f}</div></div>'
)
if mk["lat"] is not None:
known_section += (
f' <div class="stat-card"><div class="label">Marked Anchor</div>'
f' <div class="value val-accent" style="font-size:0.95rem">{mk["lat"]:.6f}, {mk["lon"]:.6f}</div></div>'
f' <div class="stat-card"><div class="label">Marked Error</div>'
f' <div class="value val-red">{mk["error_ft"]} ft</div></div>'
)
if est["lat"] is not None:
known_section += (
f' <div class="stat-card"><div class="label">Service Estimate</div>'
f' <div class="value val-accent" style="font-size:0.95rem">{est["lat"]:.6f}, {est["lon"]:.6f}</div></div>'
f' <div class="stat-card"><div class="label">Estimate Error</div>'
f' <div class="value val-red">{est["error_ft"]} ft</div></div>'
)
known_section += '</div>'
else:
known_section = ''
html = html.replace("KNOWN_ANCHOR_SECTION", known_section)
html = html.replace("ANCHOR_SET", "Yes" if p["anchor_set"] else "No")
html = html.replace("CHAIN_LEN", str(p["chain_length"]))
html = html.replace("CHAIN_BOTTOM", str(p["chain_on_bottom"]))
html = html.replace("CAT_DIST", str(p["catenary_distance"]))
html = html.replace("EST_DIST", str(p["est_distance"]))
html = html.replace("UNCERTAINTY", str(p["uncertainty"]))
html = html.replace("DRIFT", str(p["drift"]))
html = html.replace("DEPTH", str(p["vessel_depth"]))
html = html.replace("FREEBOARD", str(p["freeboard"]))
html = html.replace("WINDAGE", str(p["windage_area"] if p.get("windage_area") else ""))
html = html.replace("DRAG_COEFF", str(p.get("drag_coeff", "")))
html = html.replace("CHAIN_WT", str(p.get("chain_weight", "")))
wind_class = "val-green" if p["vessel_wind"] >= wind_threshold else "val-red"
html = html.replace("WIND_CLASS", wind_class)
html = html.replace("WIND_NOW", str(p["vessel_wind"]))
html = html.replace("WIND_DIR", str(p["vessel_wind_dir"] if p["vessel_wind_dir"] is not None else ""))
html = html.replace("HEADING", str(p["vessel_heading"] if p["vessel_heading"] is not None else ""))
html = html.replace("WIND_THR", str(wind_threshold))
html = html.replace("TENSION_RATIO", str(tension_ratio))
html = html.replace("DIST_THR", str(p["distance_threshold"]))
anchor_ref_label = "known anchor position" if ref["using_known"] else "marked anchor"
html = html.replace("ANCHOR_REF_LABEL", anchor_ref_label)
html = html.replace("ALGO_REF_LABEL", anchor_ref_label)
algos = results.get("estimators", [])
if algos:
algo_table = (
'<div class="chart-row"><div class="chart-box full">'
'<table class="param-table">'
'<tr><td style="font-weight:700;color:var(--accent)">Algorithm</td>'
'<td style="font-weight:700;color:var(--accent)">Error (ft)</td></tr>'
)
if est["lat"] is not None and est["error_ft"] is not None:
algo_table += (
f'<tr><td>Service estimate (current)</td>'
f'<td>{est["error_ft"]}</td></tr>'
)
for a in sorted(algos, key=lambda a: a["error_ft"]):
algo_table += (
f'<tr><td>{a["name"]}</td>'
f'<td>{a["error_ft"]}</td></tr>'
)
algo_table += '</table></div></div>'
else:
algo_table = '<p style="color:var(--muted)">Not enough track points for estimation algorithms.</p>'
html = html.replace("ALGO_TABLE", algo_table)
html = html.replace("N_TOTAL", str(s["total_points"]))
html = html.replace("N_RELIABLE", str(s["reliable"]))
html = html.replace("N_MARGINAL", str(s["marginal"]))
html = html.replace("N_UNRELIABLE", str(s["unreliable"]))
html = html.replace("PCT_REL", str(s["pct_reliable"]))
return html
# ---------------------------------------------------------------------------
# CLI
# ---------------------------------------------------------------------------
def main():
parser = argparse.ArgumentParser(
description="Pull anchor alarm data from Venus MQTT and analyze estimation reliability.",
)
parser.add_argument("--host", default="cerbo",
help="MQTT broker hostname (default: cerbo)")
parser.add_argument("--port", type=int, default=1883,
help="MQTT broker port (default: 1883)")
parser.add_argument("--wind-threshold", type=float, default=7.0,
help="Min wind speed (kts) for reliable classification (default: 7)")
parser.add_argument("--tension-ratio", type=float, default=0.60,
help="Min distance/chain_length ratio for reliable (default: 0.60)")
parser.add_argument("--anchor-lat", type=float, default=None,
help="Known anchor latitude (overrides service marked position)")
parser.add_argument("--anchor-lon", type=float, default=None,
help="Known anchor longitude (overrides service marked position)")
parser.add_argument("--chain-length", type=float, default=None,
help="Override chain length in feet (default: read from service)")
parser.add_argument("--timeout", type=float, default=10.0,
help="Seconds to wait for MQTT data (default: 10)")
parser.add_argument("--output", default=None,
help="Output HTML path (default: anchor_debug.html in script dir)")
args = parser.parse_args()
out_path = args.output or os.path.join(os.path.dirname(__file__), "anchor_debug.html")
collector = AnchorDataCollector(args.host, args.port)
values = collector.collect(timeout=args.timeout)
if not values:
sys.exit("No data received from MQTT. Check host and that anchor alarm service is running.")
if args.chain_length is not None:
values["/Anchor/ChainLength"] = args.chain_length
print("Analyzing ...")
results = analyze(values, args.wind_threshold, args.tension_ratio,
args.anchor_lat, args.anchor_lon)
print("Rendering HTML report ...")
html = render_report(results, args.host, args.wind_threshold, args.tension_ratio)
with open(out_path, "w") as f:
f.write(html)
s = results["stats"]
p = results["params"]
est = results["estimated"]
mk = results["marked"]
print(f"\nReport: {out_path}")
print(f" Track points: {s['total_points']}")
print(f" Reliable: {s['reliable']} ({s['pct_reliable']}%)")
print(f" Marginal: {s['marginal']}")
print(f" Unreliable: {s['unreliable']}")
print(f" Chain length: {p['chain_length']} ft")
print(f" Current wind: {p['vessel_wind']} kts")
print(f" Threshold: wind >= {args.wind_threshold} kts, dist >= {p['distance_threshold']} ft")
if args.anchor_lat is not None:
print(f" Known anchor: {args.anchor_lat:.6f}, {args.anchor_lon:.6f}")
if est["error_ft"] is not None:
print(f" Estimate err: {est['error_ft']} ft from known position")
if mk["error_ft"] is not None:
print(f" Marked error: {mk['error_ft']} ft from known position")
algos = results.get("estimators", [])
if algos:
print("\n Algorithm comparison (sorted by error):")
for a in sorted(algos, key=lambda a: a["error_ft"]):
print(f" {a['error_ft']:7.1f} ft {a['name']}")
if __name__ == "__main__":
main()