Source code: GPS data analysis in Python with MovingPandas: trajectories, speed, stops, generalisation and flow maps

"""GPS trajectory analysis in Python with MovingPandas, on one day of real AIS vessel tracks in Puget Sound.

Steps: read the GPS fixes, build trajectories, add speed and direction, detect stops, split into trips,
generalise (Douglas-Peucker against time-based), aggregate into a flow map, and round-trip a track through GPX.

Data: NOAA / BOEM MarineCadastre.gov AIS, 1 June 2023 (US public domain), and US Census Bureau cartographic
boundary county subdivisions for Washington (public domain) as the land layer. No account or key is needed.
Run: python gps_trajectories.py   (downloads about 345 MB once into data/)
"""
import time
import warnings
import zipfile
from datetime import timedelta
from pathlib import Path

warnings.filterwarnings("ignore")

import geopandas as gpd
import gpxpy
import gpxpy.gpx
import matplotlib.pyplot as plt
import movingpandas as mpd
import numpy as np
import pandas as pd
import requests
import shapely
from matplotlib.collections import LineCollection

HERE = Path(__file__).parent
DATA, OUT = HERE / "data", HERE / "output"
DATA.mkdir(exist_ok=True)
OUT.mkdir(exist_ok=True)
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})

AIS_URL = "https://coast.noaa.gov/htdata/CMSP/AISDataHandler/2023/AIS_2023_06_01.zip"
LAND_URL = "https://www2.census.gov/geo/tiger/GENZ2023/shp/cb_2023_53_cousub_500k.zip"
BBOX = (-122.65, 47.47, -122.25, 47.75)          # central Puget Sound: Seattle, Bainbridge, Bremerton, Vashon
UTM = 32610                                       # WGS 84 / UTM zone 10N, metres
FOCUS = "WSF WALLA WALLA"                         # the ferry used for the single-trajectory steps
LAND, SEA = "#e9e5dc", "#f4f8fb"


def download(url, path):
    if not path.exists():
        print(f"downloading {url}")
        with requests.get(url, stream=True, timeout=600) as r:
            r.raise_for_status()
            with open(path, "wb") as f:
                for chunk in r.iter_content(1 << 20):
                    f.write(chunk)
    return path


def basemap(ax, land, extent=None):
    ax.set_facecolor(SEA)
    land.plot(ax=ax, color=LAND, edgecolor="#b8b2a6", linewidth=0.4, zorder=0)
    x0, y0, x1, y1 = extent or AREA
    ax.set_xlim(x0, x1)
    ax.set_ylim(y0, y1)
    ax.set_xticks([])
    ax.set_yticks([])
    ax.set_xlabel("")
    ax.set_ylabel("")
    ax.set_aspect("equal")


def save(fig, name):
    for ax in fig.axes:
        if "[metre]" in ax.get_xlabel():
            ax.set_xlabel("")
            ax.set_ylabel("")
    fig.savefig(OUT / name)
    plt.close(fig)

def km_bar(ax, length=5000):
    x0, x1 = ax.get_xlim()
    y0, y1 = ax.get_ylim()
    x, y = x0 + 0.04 * (x1 - x0), y0 + 0.05 * (y1 - y0)
    ax.plot([x, x + length], [y, y], color="#222", lw=2.5, solid_capstyle="butt")
    ax.text(x + length / 2, y + 0.015 * (y1 - y0), f"{length / 1000:g} km", ha="center", va="bottom", fontsize=8)


# ---------------------------------------------------------------- 1. read the GPS fixes
t0 = time.perf_counter()
zip_path = download(AIS_URL, DATA / "AIS_2023_06_01.zip")
cache = DATA / "puget_sound_ais_2023-06-01.csv"
if not cache.exists():
    with zipfile.ZipFile(zip_path) as z:
        parts = []
        for chunk in pd.read_csv(z.open(z.namelist()[0]), chunksize=1_000_000):
            keep = chunk.LON.between(BBOX[0], BBOX[2]) & chunk.LAT.between(BBOX[1], BBOX[3])
            parts.append(chunk[keep])
    pd.concat(parts).to_csv(cache, index=False)
ais = pd.read_csv(cache, parse_dates=["BaseDateTime"])
print(f"AIS fixes in the Puget Sound box on 1 June 2023: {len(ais):,} from {ais.MMSI.nunique()} vessels "
      f"(read in {time.perf_counter() - t0:.1f} s)")
print(ais[["MMSI", "BaseDateTime", "LAT", "LON", "SOG", "COG", "VesselName", "VesselType"]].head(3).to_string(index=False))

n0 = len(ais)
ais = ais.drop_duplicates(["MMSI", "BaseDateTime"]).sort_values(["MMSI", "BaseDateTime"])
print(f"duplicate timestamps dropped: {n0 - len(ais):,}")
gaps = ais.groupby("MMSI").BaseDateTime.diff().dt.total_seconds().dropna()
print(f"time between fixes: median {gaps.median():.0f} s, 90th pct {gaps.quantile(0.9):.0f} s, max {gaps.max() / 3600:.1f} h")

land_zip = download(LAND_URL, DATA / "cb_2023_53_cousub_500k.zip")
land = gpd.read_file(land_zip).to_crs(UTM)
pts = gpd.GeoDataFrame(ais, geometry=gpd.points_from_xy(ais.LON, ais.LAT), crs=4326).to_crs(UTM)
AREA = tuple(gpd.GeoSeries([shapely.box(*BBOX)], crs=4326).to_crs(UTM).total_bounds)
land = land.clip(shapely.box(*AREA).buffer(2000))

# ---------------------------------------------------------------- 2. build trajectories
pts = pts.rename(columns={"BaseDateTime": "t"})
t0 = time.perf_counter()
tc = mpd.TrajectoryCollection(pts, "MMSI", t="t", min_length=1000)
print(f"\nTrajectoryCollection: {len(tc)} trajectories longer than 1 km, built in {time.perf_counter() - t0:.1f} s "
      f"({ais.MMSI.nunique() - len(tc)} vessels moved less than 1 km and were left out)")
# a vessel that leaves the box and comes back later would be joined by a straight line across land
tc = mpd.ObservationGapSplitter(tc).split(gap=timedelta(minutes=15), min_length=1000)
print(f"ObservationGapSplitter(gap=15 min): {len(tc)} trajectories after splitting at gaps in reception")
focus = next(tr for tr in tc if tr.df.VesselName.iloc[0] == FOCUS)
print(f"{FOCUS}: {sum(tr.df.VesselName.iloc[0] == FOCUS for tr in tc)} trajectory, {len(focus.df):,} fixes, "
      f"{focus.get_length() / 1000:.1f} km, {focus.get_start_time():%H:%M} to {focus.get_end_time():%H:%M} UTC")

fig, axs = plt.subplots(1, 2, figsize=(11, 5.2))
basemap(axs[0], land)
pts.plot(ax=axs[0], markersize=0.2, color="#3288bd", alpha=0.35)
axs[0].set_title(f"(a) {len(pts):,} AIS fixes: unordered points")
basemap(axs[1], land)
lines = tc.to_traj_gdf()
lines.plot(ax=axs[1], linewidth=0.5, color="#5e4fa2", alpha=0.6)
focus.to_line_gdf().plot(ax=axs[1], linewidth=1.6, color="#d53e4f")
axs[1].set_title(f"(b) {len(tc)} trajectories (red: {FOCUS.title()})")
km_bar(axs[0])
km_bar(axs[1])
save(fig, "01-points-to-trajectories.png")

# ---------------------------------------------------------------- 3. speed and direction
focus.add_speed(overwrite=True, units=("km", "h"))
focus.add_direction(overwrite=True)
f = focus.df
sog_kmh = f.SOG * 1.852
moving = sog_kmh > 5
diff = (f.speed - sog_kmh)[moving]
ang = ((f.direction - f.COG + 180) % 360 - 180)[moving]
print(f"\n{FOCUS} speed from positions vs reported SOG (moving fixes, n={moving.sum()}): "
      f"median |difference| {diff.abs().median():.2f} km/h, r = {np.corrcoef(f.speed[moving], sog_kmh[moving])[0, 1]:.3f}")
print(f"direction from positions vs reported COG: median |difference| {ang.abs().median():.1f} degrees")
print(f"max computed speed {f.speed.max():.1f} km/h ({f.speed.max() / 1.852:.1f} kn); "
      f"max reported SOG {f.SOG.max():.1f} kn")

tc.add_speed(overwrite=True, units=("km", "h"))
allp = tc.to_point_gdf()
fast = allp[allp.speed > 100]
print(f"fixes implying more than 100 km/h across all vessels: {len(fast)} in {fast.MMSI.nunique()} trajectories")
cleaned = mpd.OutlierCleaner(tc).clean(v_max=100, units=("km", "h"))
print(f"OutlierCleaner(v_max=100 km/h) removed {len(allp) - len(cleaned.to_point_gdf()):,} fixes")
tc = cleaned
focus_clean = mpd.OutlierCleaner(focus).clean(v_max=45, units=("km", "h"))
focus_clean.add_speed(overwrite=True, units=("km", "h"))
print(f"OutlierCleaner(v_max=45 km/h) on {FOCUS}: removed {len(focus.df) - len(focus_clean.df)} fixes; "
      f"max speed now {focus_clean.df.speed.max():.1f} km/h")

fig = plt.figure(figsize=(11, 7))
gs = fig.add_gridspec(2, 2, height_ratios=[0.9, 1], hspace=0.28)
ax = fig.add_subplot(gs[0, :])
basemap(ax, land)
ln = focus.to_line_gdf()
ln["speed"] = focus.df.speed.iloc[1:].values
seg = LineCollection([np.asarray(g.coords) for g in ln.geometry], cmap="Spectral_r", linewidths=2)
seg.set_array(ln.speed.values)
seg.set_clim(0, 40)
ax.add_collection(seg)
x0, y0, x1, y1 = focus.to_traj_gdf().total_bounds
ax.set_xlim(x0 - 1500, x1 + 1500)
ax.set_ylim(y0 - 1500, y1 + 1500)
fig.colorbar(seg, ax=ax, shrink=0.8, pad=0.01, label="speed from positions (km/h)")
km_bar(ax, 2000)
ax.set_title(f"(a) {FOCUS.title()}, Bainbridge Island (left) to Seattle (right): speed from add_speed()")
ax = fig.add_subplot(gs[1, 0])
ax.plot(f.index, sog_kmh, color="#3288bd", lw=1.6, label="reported SOG")
ax.plot(f.index, f.speed, color="#d53e4f", lw=0.8, label="add_speed()")
ax.set_ylabel("km/h")
ax.legend(loc="upper center", fontsize=7, ncol=2)
ax.set_ylim(0, 70)
ax.set_title("(b) speed: from positions vs reported")
ax.set_xlabel("1 June 2023 (UTC)")
ax.xaxis.set_major_formatter(plt.matplotlib.dates.DateFormatter("%H:%M"))
ax = fig.add_subplot(gs[1, 1], sharex=ax)
ax.scatter(f.index[moving], f.COG[moving], s=4, color="#3288bd", label="reported COG")
ax.scatter(f.index[moving], f.direction[moving], s=1.5, color="#d53e4f", label="add_direction()")
ax.set_ylabel("degrees from north")
ax.set_yticks([0, 90, 180, 270, 360])
ax.set_ylim(0, 400)
ax.legend(loc="upper center", fontsize=7, ncol=2, markerscale=3)
ax.set_title("(c) direction while moving (> 5 km/h)")
ax.set_xlabel("1 June 2023 (UTC)")
ax.xaxis.set_major_formatter(plt.matplotlib.dates.DateFormatter("%H:%M"))
save(fig, "02-speed-direction.png")
focus = focus_clean

# ---------------------------------------------------------------- 4. stop detection
ferries = mpd.TrajectoryCollection([tr for tr in tc if str(tr.df.VesselName.iloc[0]).startswith("WSF ")])
print(f"\nWashington State Ferries in the box: {len(ferries)} trajectories from "
      + ", ".join(sorted({tr.df.VesselName.iloc[0].replace('WSF ', '').title() for tr in ferries})))
t0 = time.perf_counter()
detector = mpd.TrajectoryStopDetector(ferries)
stops = detector.get_stop_points(max_diameter=150, min_duration=timedelta(minutes=5))
print(f"TrajectoryStopDetector(max_diameter=150 m, min_duration=5 min): {len(stops)} stops in "
      f"{time.perf_counter() - t0:.1f} s")

# group the stop points into places: stops within 500 m of each other belong to the same terminal
blobs = gpd.GeoDataFrame(geometry=list(stops.buffer(250).union_all().geoms), crs=UTM)
blobs["place"] = range(len(blobs))
stops = gpd.sjoin(stops, blobs, predicate="within").drop(columns="index_right")
TERMINALS = {"Seattle (Colman Dock)": (-122.3385, 47.6025), "Bainbridge Island": (-122.5097, 47.6229),
             "Bremerton": (-122.6237, 47.5623), "Fauntleroy": (-122.3963, 47.5232),
             "Vashon": (-122.4637, 47.5108), "Southworth": (-122.4957, 47.5130)}
term = gpd.GeoDataFrame({"name": list(TERMINALS)}, geometry=gpd.points_from_xy(*zip(*TERMINALS.values())),
                        crs=4326).to_crs(UTM)
places = stops.dissolve("place", aggfunc={"duration_s": ["count", "median", "sum"]})
places.columns = ["geometry", "n_stops", "median_dwell_s", "total_s"]
places["geometry"] = places.centroid
places = places.set_crs(UTM, allow_override=True)
places = gpd.sjoin_nearest(places, term, max_distance=600, how="left").drop(columns="index_right")
places["name"] = places.name.fillna("unnamed")
places = places.sort_values("n_stops", ascending=False)
print("stop places (stops within 500 m grouped):")
print(places.assign(median_dwell_min=(places.median_dwell_s / 60).round(1), hours_stopped=(places.total_s / 3600).round(1),
                    lon=places.to_crs(4326).geometry.x.round(4), lat=places.to_crs(4326).geometry.y.round(4))
      [["name", "n_stops", "median_dwell_min", "hours_stopped", "lon", "lat"]].to_string(index=False))

fig, ax = plt.subplots(figsize=(8.5, 6.4))
basemap(ax, land)
ferries.to_traj_gdf().plot(ax=ax, color="#5e4fa2", lw=0.4, alpha=0.5)
sc = ax.scatter(stops.geometry.x, stops.geometry.y, s=8 + stops.duration_s.clip(upper=3600) / 60 * 3,
                c=stops.duration_s / 60, cmap="Spectral_r", vmin=5, vmax=60, edgecolor="#333", linewidth=0.3,
                alpha=0.85, zorder=3)
OFFSETS = {"Seattle (Colman Dock)": (8, 10), "Bainbridge Island": (-10, 14), "Bremerton": (6, 12),
           "Fauntleroy": (10, 4), "Vashon": (14, -26), "Southworth": (-70, 14)}
for _, p in places.iterrows():
    label = p["name"] if p["name"] != "unnamed" else "berthed all day"
    dwell = (f"{p.n_stops} stops, median {p.median_dwell_s / 60:.0f} min" if p.n_stops > 1
             else f"1 stop of {p.median_dwell_s / 3600:.0f} h")
    ax.annotate(f"{label}\n{dwell}",
                (p.geometry.x, p.geometry.y), xytext=OFFSETS.get(p["name"], (8, -22)), textcoords="offset points",
                fontsize=7, bbox=dict(boxstyle="square,pad=0.2", fc="white", ec="none", alpha=0.85), zorder=4)
fig.colorbar(sc, ax=ax, shrink=0.6, label="stop duration (minutes; marker size capped at 60)", extend="max")
km_bar(ax)
ax.set_title(f"{len(stops)} ferry stops detected (max_diameter 150 m, min_duration 5 min)")
save(fig, "03-stop-detection.png")

# ---------------------------------------------------------------- 5. splitting into trips
trips = mpd.StopSplitter(ferries).split(max_diameter=150, min_duration=timedelta(minutes=5), min_length=1000)
print(f"\nStopSplitter(max_diameter=150 m, min_duration=5 min, min_length=1 km): {len(ferries)} ferry "
      f"trajectories -> {len(trips)} trips")


def nearest_place(pt):
    d = places.distance(pt)
    return places.name.iloc[int(np.argmin(d.values))] if d.min() < 1500 else "edge of box"


rows = []
for tr in trips:
    rows.append({"ferry": tr.df.VesselName.iloc[0].replace("WSF ", "").title(),
                 "from": nearest_place(tr.get_start_location()), "to": nearest_place(tr.get_end_location()),
                 "minutes": tr.get_duration().total_seconds() / 60, "km": tr.get_length() / 1000,
                 "traj": tr})
trip_df = pd.DataFrame(rows)
trip_df["route"] = trip_df.apply(lambda r: " - ".join(sorted([r["from"], r["to"]])), axis=1)
summary = (trip_df.groupby("route").agg(trips=("minutes", "size"), median_min=("minutes", "median"),
                                        median_km=("km", "median")).sort_values("trips", ascending=False))
print(summary.round(1).to_string())
ft = trip_df[trip_df.ferry == FOCUS.replace("WSF ", "").title()]
print(f"{FOCUS} trips: {len(ft)}, median crossing {ft.minutes.median():.0f} min over {ft.km.median():.1f} km")

fig = plt.figure(figsize=(10, 9))
gs = fig.add_gridspec(2, 1, height_ratios=[1.35, 1], hspace=0.22)
ax0 = fig.add_subplot(gs[0])
basemap(ax0, land)
cmap = plt.get_cmap("Spectral")
for i, (_, r) in enumerate(ft.iterrows()):
    r.traj.to_traj_gdf().plot(ax=ax0, color=cmap(i / max(len(ft) - 1, 1)), lw=1.3)
x0, y0, x1, y1 = focus.to_traj_gdf().total_bounds
ax0.set_xlim(x0 - 1500, x1 + 1500)
ax0.set_ylim(y0 - 1500, y1 + 1500)
km_bar(ax0, 2000)
ax0.set_title(f"(a) {FOCUS.title()} split into {len(ft)} trips, one colour each")
ax1 = fig.add_subplot(gs[1])
top = summary[summary.trips >= 4].index
ROUTE_COLOURS = ["#d53e4f", "#f46d43", "#66c2a5", "#3288bd", "#5e4fa2", "#999999"]
for i, route in enumerate(top):
    v = trip_df[trip_df.route == route].minutes
    ax1.hist(v, bins=np.arange(0, 80, 2.5), alpha=0.8, color=ROUTE_COLOURS[i % len(ROUTE_COLOURS)],
             label=f"{route} ({len(v)} trips)")
ax1.set_xlabel("trip duration (minutes)")
ax1.set_ylabel("trips")
ax1.legend(fontsize=7.5)
ax1.set_title("(b) duration of every ferry trip by route")
save(fig, "04-split-trips.png")

# ---------------------------------------------------------------- 6. generalisation: Douglas-Peucker vs time-based
fdf = focus.df
tt = (fdf.index - fdf.index[0]).total_seconds().values
xy = np.c_[fdf.geometry.x, fdf.geometry.y]


def errors(gen):
    """Spatial error (distance to the generalised line) and synchronous error (distance to where the
    generalised trajectory says the vessel was at the same moment), for every original fix."""
    g = gen.df
    gt = (g.index - fdf.index[0]).total_seconds().values
    gx, gy = np.interp(tt, gt, g.geometry.x), np.interp(tt, gt, g.geometry.y)
    sync = np.hypot(xy[:, 0] - gx, xy[:, 1] - gy)
    line = gen.to_linestring()
    spatial = shapely.distance(shapely.points(xy), line)
    errors.last_sync_xy = np.c_[gx, gy]
    return len(g), spatial, sync


sweep = []
methods = {"Douglas-Peucker": (mpd.DouglasPeuckerGeneralizer, [2, 5, 10, 20, 50, 100, 200, 400], "m"),
           "Top-down time ratio": (mpd.TopDownTimeRatioGeneralizer, [2, 5, 10, 20, 50, 100, 200, 400], "m"),
           "Min. time delta": (mpd.MinTimeDeltaGeneralizer,
                               [timedelta(minutes=m) for m in (2, 3, 5, 10, 15, 20, 30, 60)], "min")}
for name, (cls, tols, unit) in methods.items():
    for tol in tols:
        t0 = time.perf_counter()
        gen = cls(focus).generalize(tolerance=tol)
        secs = time.perf_counter() - t0
        n, sp, sy = errors(gen)
        sweep.append({"method": name, "tolerance": tol if unit == "m" else tol.total_seconds() / 60,
                      "unit": unit, "points": n, "kept_pct": 100 * n / len(fdf), "spatial_max_m": sp.max(),
                      "spatial_mean_m": sp.mean(), "sync_mean_m": sy.mean(), "sync_max_m": sy.max(), "secs": secs})
sweep = pd.DataFrame(sweep)
print(f"\ngeneralising {FOCUS} ({len(fdf):,} fixes):")
print(sweep[["method", "tolerance", "unit", "points", "kept_pct", "spatial_max_m", "spatial_mean_m",
             "sync_mean_m", "sync_max_m", "secs"]].round(2).to_string(index=False))

picks = {"Douglas-Peucker": 50, "Top-down time ratio": 50, "Min. time delta": 10}
crossing = ft.sort_values("minutes").iloc[len(ft) // 2].traj
c0, c1 = crossing.get_start_time(), crossing.get_end_time()
win = (fdf.index >= c0) & (fdf.index <= c1)
appr = fdf[(fdf.index >= c0) & (fdf.index <= c0 + timedelta(minutes=9))]
ax_, ay_ = appr.geometry.x.mean(), appr.geometry.y.mean()
half = 0.6 * max(np.ptp(appr.geometry.x), np.ptp(appr.geometry.y)) + 300
zoom = (ax_ - half, ay_ - half, ax_ + half, ay_ + half)
print(f"figure 5 zooms on the Bainbridge approach of one crossing: {c0:%H:%M}-{c1:%H:%M} UTC, "
      f"{win.sum()} original fixes")
fig = plt.figure(figsize=(11, 8.4))
gs = fig.add_gridspec(2, 6, height_ratios=[1, 0.95], hspace=0.3, wspace=0.6)
for i, (name, tol) in enumerate(picks.items()):
    cls = methods[name][0]
    gen = cls(focus).generalize(tolerance=tol if name != "Min. time delta" else timedelta(minutes=tol))
    n, sp, sy = errors(gen)
    sxy = errors.last_sync_xy
    ax = fig.add_subplot(gs[0, 2 * i:2 * i + 2])
    basemap(ax, land, extent=zoom)
    for (ox, oy), (sx, sy_) in zip(xy[win], sxy[win]):
        ax.plot([ox, sx], [oy, sy_], color="#3288bd", lw=0.9, zorder=2)
    ax.plot(xy[win, 0], xy[win, 1], color="#bbb", lw=2.2, zorder=1)
    ax.scatter(xy[win, 0], xy[win, 1], s=6, color="#555", zorder=2)
    g = gen.df
    gw = (g.index >= c0) & (g.index <= c1)
    ax.plot(g.geometry.x[gw], g.geometry.y[gw], color="#d53e4f", lw=1.1, marker="o", ms=3.5, zorder=3)
    km_bar(ax, 500)
    unit = "min" if name == "Min. time delta" else "m"
    ax.set_title(f"{name}, {tol} {unit}\n{n} of {len(fdf)} fixes kept\nmean sync. error, this crossing: "
                 f"{sy[win].mean():.0f} m", fontsize=8)
for j, (metric, label) in enumerate([("spatial_max_m", "max. distance to generalised line (m)"),
                                     ("sync_mean_m", "mean synchronous error (m)")]):
    ax = fig.add_subplot(gs[1, 3 * j:3 * j + 3])
    for k, (name, grp) in enumerate(sweep.groupby("method", sort=False)):
        ax.plot(grp.kept_pct, grp[metric], marker="o", ms=4, color=["#d53e4f", "#3288bd", "#66a61e"][k], label=name)
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.set_xlabel("fixes kept (%)")
    ax.set_ylabel(label)
    ax.grid(alpha=0.3, which="both")
    ax.legend(fontsize=7)
fig.suptitle(f"Generalising {FOCUS.title()}: grey = original fixes, red = kept fixes, "
             f"blue = offset to the synchronous position", y=0.97, fontsize=9.5)
save(fig, "05-generalisation.png")

# ---------------------------------------------------------------- 7. aggregation into a flow map
movers = mpd.TrajectoryCollection([tr for tr in tc if tr.df.SOG.max() > 3])
gen_tc = mpd.MinDistanceGeneralizer(movers).generalize(tolerance=100)
t0 = time.perf_counter()
agg = mpd.TrajectoryCollectionAggregator(gen_tc, max_distance=2000, min_distance=200,
                                         min_stop_duration=timedelta(minutes=5))
sig = agg.get_significant_points_gdf()
clusters = agg.get_clusters_gdf()
flows = agg.get_flows_gdf()
print(f"\nflow map from {len(movers)} vessels that moved faster than 3 kn: {len(sig):,} significant points, "
      f"{len(clusters)} clusters, {len(flows)} flows in {time.perf_counter() - t0:.1f} s")
flows = flows.sort_values("weight")
print(f"flow weights: median {flows.weight.median():.0f}, max {flows.weight.max():.0f}; "
      f"flows with weight >= 5: {(flows.weight >= 5).sum()}")

fig, axs = plt.subplots(1, 2, figsize=(11, 5.2))
fig.subplots_adjust(right=0.9, wspace=0.04)
basemap(axs[0], land)
movers.to_traj_gdf().plot(ax=axs[0], lw=0.4, color="#5e4fa2", alpha=0.5)
sig.plot(ax=axs[0], markersize=1.5, color="#f46d43", alpha=0.6)
clusters.plot(ax=axs[0], markersize=12, color="#222", marker="x")
axs[0].set_title(f"(a) {len(sig):,} significant points, {len(clusters)} cluster centres")
basemap(axs[1], land)
w = flows.weight.clip(upper=60)
cax = fig.add_axes([0.915, 0.25, 0.012, 0.5])
flows.plot(ax=axs[1], column="weight", cmap="Spectral_r", linewidth=0.4 + w / 8, vmin=1, vmax=60, legend=True,
           cax=cax, legend_kwds={"label": "trajectories per flow"})
axs[1].set_title(f"(b) {len(flows)} aggregated flows between clusters")
km_bar(axs[0])
km_bar(axs[1])
save(fig, "06-flow-map.png")

# ---------------------------------------------------------------- 8. GPX in and out
trip = ft.sort_values("minutes").iloc[len(ft) // 2].traj.to_crs(4326)
gpx = gpxpy.gpx.GPX()
track = gpxpy.gpx.GPXTrack(name=f"{FOCUS.title()} {trip.get_start_time():%Y-%m-%d %H:%M} UTC")
gpx.tracks.append(track)
seg = gpxpy.gpx.GPXTrackSegment()
track.segments.append(seg)
for t, p in zip(trip.df.index, trip.df.geometry):
    seg.points.append(gpxpy.gpx.GPXTrackPoint(p.y, p.x, time=t.to_pydatetime()))
gpx_path = OUT / "ferry_trip.gpx"
gpx_path.write_text(gpx.to_xml(), encoding="utf-8")
print(f"\nwrote {gpx_path.name}: {len(seg.points)} track points, {gpx_path.stat().st_size / 1024:.1f} KB")

gpx_pts = gpd.read_file(gpx_path, layer="track_points")
print(f"read back with geopandas (layer='track_points'): {len(gpx_pts)} rows, columns include "
      f"{[c for c in ('track_fid', 'track_seg_point_id', 'time', 'geometry') if c in gpx_pts.columns]}")
gpx_pts["time"] = pd.to_datetime(gpx_pts.time).dt.tz_localize(None)
gpx_traj = mpd.Trajectory(gpx_pts.to_crs(UTM), traj_id="gpx", t="time")
gpx_traj.add_speed(units=("km", "h"))
print(f"GPX trajectory: {gpx_traj.get_length() / 1000:.2f} km in {gpx_traj.get_duration()}, "
      f"max speed {gpx_traj.df.speed.max():.1f} km/h (original trip {trip.to_crs(UTM).get_length() / 1000:.2f} km)")

with open(gpx_path, encoding="utf-8") as fh:
    parsed = gpxpy.parse(fh)
first = parsed.tracks[0].segments[0].points[0]
print(f"gpxpy: {parsed.tracks[0].name}; first point {first.latitude:.5f}, {first.longitude:.5f} at {first.time}; "
      f"gpxpy length_2d {parsed.length_2d() / 1000:.2f} km")