shuishen
2 days ago 4092854b0bea9e1fc02f29222c6f7cd876e565a5
capabilities/15-trajectory-analysis/run_trajectory_analysis.py
@@ -25,6 +25,7 @@
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib.patches import Patch
DEFAULT_THRESHOLDS = {
@@ -105,7 +106,7 @@
    return CRS.from_epsg(epsg)
def load_case(manifest_path: Path) -> tuple[dict[str, Any], pd.DataFrame, gpd.GeoDataFrame, gpd.GeoDataFrame, int, CRS]:
def load_case(manifest_path: Path) -> tuple[dict[str, Any], pd.DataFrame, gpd.GeoDataFrame, gpd.GeoDataFrame, gpd.GeoDataFrame | None, int, CRS]:
    manifest = load_json(manifest_path)
    for field in ("case_id", "crs", "observations", "reference_routes", "zones"):
        if not manifest.get(field):
@@ -142,6 +143,12 @@
        raise ValueError("Every track must contain at least two distinct timestamps.")
    routes = gpd.read_file(routes_path)
    zones = gpd.read_file(zones_path)
    flyable_zones: gpd.GeoDataFrame | None = None
    if manifest.get("flyable_zones"):
        flyable_path = resolve_input_path(manifest_path, str(manifest["flyable_zones"]))
        flyable_zones = gpd.read_file(flyable_path)
        if flyable_zones.crs is None:
            raise ValueError("Flyable-zone GeoJSON must declare a CRS.")
    if "track_id" not in routes.columns:
        raise ValueError("Reference route GeoJSON must contain a track_id property.")
    if "zone_id" not in zones.columns or "zone_type" not in zones.columns:
@@ -153,7 +160,15 @@
    if missing_routes:
        raise ValueError(f"Missing reference routes for: {', '.join(sorted(missing_routes))}")
    metric_crs = choose_metric_crs(float(frame["longitude"].mean()), float(frame["latitude"].mean()))
    return manifest, frame, routes.to_crs(metric_crs), zones.to_crs(metric_crs), duplicate_count, metric_crs
    return (
        manifest,
        frame,
        routes.to_crs(metric_crs),
        zones.to_crs(metric_crs),
        flyable_zones.to_crs(metric_crs) if flyable_zones is not None else None,
        duplicate_count,
        metric_crs,
    )
def true_runs(flags: list[bool], times: list[pd.Timestamp], max_gap: float) -> Iterable[tuple[int, int]]:
@@ -410,6 +425,7 @@
    tracks: dict[str, gpd.GeoDataFrame],
    routes_wgs84: gpd.GeoDataFrame,
    zones_wgs84: gpd.GeoDataFrame,
    flyable_zones_wgs84: gpd.GeoDataFrame | None,
    metadata: dict[str, Any],
) -> None:
    output_dir.mkdir(parents=True, exist_ok=True)
@@ -435,7 +451,22 @@
        for event in events
    ]
    write_geojson(output_dir / "events.geojson", event_features)
    render_map(output_dir / "analysis.png", tracks, routes_wgs84, zones_wgs84, events, manifest["case_id"])
    route_features = json.loads(routes_wgs84.to_json())["features"]
    write_geojson(output_dir / "reference_routes.geojson", route_features)
    zone_features = json.loads(zones_wgs84.to_json())["features"]
    write_geojson(output_dir / "zones.geojson", zone_features)
    if flyable_zones_wgs84 is not None:
        flyable_features = json.loads(flyable_zones_wgs84.to_json())["features"]
        write_geojson(output_dir / "flyable_zones.geojson", flyable_features)
    render_map(
        output_dir / "analysis.png",
        tracks,
        routes_wgs84,
        zones_wgs84,
        flyable_zones_wgs84,
        events,
        manifest["case_id"],
    )
    (output_dir / "run_metadata.json").write_text(
        json.dumps(metadata, ensure_ascii=False, indent=2), encoding="utf-8"
    )
@@ -457,17 +488,20 @@
    tracks: dict[str, gpd.GeoDataFrame],
    routes: gpd.GeoDataFrame,
    zones: gpd.GeoDataFrame,
    flyable_zones: gpd.GeoDataFrame | None,
    events: list[dict[str, Any]],
    case_id: str,
) -> None:
    figure, axis = plt.subplots(figsize=(10, 7), dpi=140)
    zones.plot(ax=axis, facecolor="#f0e442", edgecolor="#8c6d1f", alpha=0.28, linewidth=1.5)
    routes.plot(ax=axis, color="#666666", linestyle="--", linewidth=1.2)
    if flyable_zones is not None and not flyable_zones.empty:
        flyable_zones.plot(ax=axis, facecolor="#36a269", edgecolor="#167344", alpha=0.22, linewidth=1.3)
    palette = ["#0072b2", "#009e73", "#d55e00", "#cc79a7", "#56b4e9"]
    for index, (track_id, track) in enumerate(sorted(tracks.items())):
        color = palette[index % len(palette)]
        axis.plot(track.geometry.x, track.geometry.y, color=color, linewidth=2.2, marker="o", markersize=2.8, label=track_id)
        axis.annotate(track_id, (track.geometry.x.iloc[-1], track.geometry.y.iloc[-1]), fontsize=7, color=color)
    routes.plot(ax=axis, color="#444444", linestyle="--", linewidth=1.5, zorder=4, label="reference route")
    for event in events:
        event_type = event["event_type"]
        axis.scatter(
@@ -481,6 +515,11 @@
            zorder=5,
        )
    handles, labels = axis.get_legend_handles_labels()
    handles.append(Patch(facecolor="#f0e442", edgecolor="#8c6d1f", alpha=0.45, label="restricted area"))
    labels.append("restricted area")
    if flyable_zones is not None and not flyable_zones.empty:
        handles.append(Patch(facecolor="#36a269", edgecolor="#167344", alpha=0.45, label="flyable area"))
        labels.append("flyable area")
    present_event_types = {event["event_type"] for event in events}
    handles.extend(
        Line2D(
@@ -503,6 +542,13 @@
    axis.set_ylabel("Latitude (WGS84)")
    axis.grid(alpha=0.18)
    axis.ticklabel_format(useOffset=False)
    visible_geometries = [geometry for track in tracks.values() for geometry in track.geometry]
    visible_geometries.extend(geometry for geometry in routes.geometry if geometry is not None)
    bounds = gpd.GeoSeries(visible_geometries, crs="EPSG:4326").total_bounds
    width = max(bounds[2] - bounds[0], 0.0005)
    height = max(bounds[3] - bounds[1], 0.0005)
    axis.set_xlim(bounds[0] - width * 0.08, bounds[2] + width * 0.08)
    axis.set_ylim(bounds[1] - height * 0.08, bounds[3] + height * 0.08)
    figure.tight_layout()
    figure.savefig(path)
    plt.close(figure)
@@ -510,7 +556,7 @@
def process_manifest(manifest_path: Path, output_dir: Path) -> dict[str, Any]:
    started = time.perf_counter()
    manifest, frame, routes_metric, zones_metric, duplicates, metric_crs = load_case(manifest_path)
    manifest, frame, routes_metric, zones_metric, flyable_zones_metric, duplicates, metric_crs = load_case(manifest_path)
    thresholds = dict(DEFAULT_THRESHOLDS)
    for key, value in manifest.get("thresholds", {}).items():
        if key not in thresholds:
@@ -566,6 +612,7 @@
        tracks,
        routes_metric.to_crs("EPSG:4326"),
        zones_metric.to_crs("EPSG:4326"),
        flyable_zones_metric.to_crs("EPSG:4326") if flyable_zones_metric is not None else None,
        metadata,
    )
    return metadata