"""Prepare one DJI flight log, WPMZ route and no-fly dataset for analysis."""
|
|
from __future__ import annotations
|
|
import argparse
|
import gzip
|
import hashlib
|
import json
|
import sys
|
import warnings
|
import xml.etree.ElementTree as element_tree
|
from datetime import UTC
|
from pathlib import Path
|
from zipfile import ZipFile
|
|
import geopandas as gpd
|
import numpy as np
|
import pandas as pd
|
from shapely import make_valid
|
from shapely.geometry import LineString, Polygon, mapping
|
from shapely.ops import unary_union
|
|
|
WGS84 = "EPSG:4326"
|
LOCAL_TIMEZONE = "Asia/Shanghai"
|
KML_NAMESPACE = "http://www.opengis.net/kml/2.2"
|
WPML_NAMESPACE = "http://www.dji.com/wpmz/1.0.5"
|
EXCEL_COLUMNS = {
|
"mission_id": "\u98de\u884c\u4efb\u52a1ID",
|
"latitude": "\u7eac\u5ea6",
|
"longitude": "\u7ecf\u5ea6",
|
"absolute_height_m": "\u7edd\u5bf9\u9ad8\u5ea6(m)",
|
"relative_height_m": "\u5b9e\u65f6\u771f\u9ad8(m)",
|
"timestamp": "\u521b\u5efa\u65f6\u95f4",
|
"flight_phase": "\u98de\u884c\u7c7b\u578b",
|
}
|
|
|
def parse_args() -> argparse.Namespace:
|
root = Path(__file__).resolve().parents[2]
|
parser = argparse.ArgumentParser(description="Prepare one real DJI flight for trajectory analysis.")
|
parser.add_argument(
|
"--raw-dir",
|
type=Path,
|
default=root / "shared" / "data" / "raw" / "15-trajectory-analysis" / "tian-dun-demo-20260814",
|
)
|
parser.add_argument(
|
"--output",
|
type=Path,
|
default=root / "shared" / "data" / "processed" / "15-trajectory-analysis" / "tian-dun-flight-19578",
|
)
|
parser.add_argument("--case-id", default="tian-dun-flight-19578")
|
parser.add_argument("--overwrite", action="store_true", help="Update derived files in an existing output directory.")
|
parser.add_argument(
|
"--zone-search-radius-m",
|
type=float,
|
default=5000.0,
|
help="Keep source no-fly areas intersecting this buffer around the flight and route.",
|
)
|
return parser.parse_args()
|
|
|
def discover_single(directory: Path, suffix: str) -> Path:
|
matches = sorted(directory.glob(f"*{suffix}"))
|
if len(matches) != 1:
|
raise ValueError(f"Expected exactly one {suffix} file in {directory}, found {len(matches)}.")
|
return matches[0]
|
|
|
def sha256(path: Path) -> str:
|
digest = hashlib.sha256()
|
with path.open("rb") as stream:
|
for chunk in iter(lambda: stream.read(1024 * 1024), b""):
|
digest.update(chunk)
|
return digest.hexdigest()
|
|
|
def choose_metric_crs(longitude: float, latitude: float) -> str:
|
zone = int((longitude + 180) // 6) + 1
|
return f"EPSG:{(32600 if latitude >= 0 else 32700) + zone}"
|
|
|
def load_observations(path: Path) -> tuple[pd.DataFrame, dict[str, object]]:
|
warnings.filterwarnings("ignore", message="Workbook contains no default style")
|
source = pd.read_excel(path)
|
missing = [column for column in EXCEL_COLUMNS.values() if column not in source.columns]
|
if missing:
|
raise ValueError(f"Flight workbook is missing columns: {', '.join(missing)}")
|
missions = source[EXCEL_COLUMNS["mission_id"]].dropna().unique().tolist()
|
if len(missions) != 1:
|
raise ValueError(f"This single-flight preparer requires one mission ID, found: {missions}")
|
mission_id = str(missions[0]).strip()
|
timestamps = pd.to_datetime(source[EXCEL_COLUMNS["timestamp"]], errors="raise")
|
if getattr(timestamps.dt, "tz", None) is None:
|
timestamps = timestamps.dt.tz_localize(LOCAL_TIMEZONE, ambiguous="raise", nonexistent="raise")
|
timestamps = timestamps.dt.tz_convert(UTC)
|
longitude = pd.to_numeric(source[EXCEL_COLUMNS["longitude"]], errors="raise")
|
latitude = pd.to_numeric(source[EXCEL_COLUMNS["latitude"]], errors="raise")
|
valid = longitude.between(-180, 180) & latitude.between(-90, 90)
|
if not valid.all():
|
raise ValueError("Flight workbook contains longitude/latitude values outside WGS84 bounds.")
|
track_id = f"drone-{mission_id}"
|
observations = pd.DataFrame(
|
{
|
"track_id": track_id,
|
"entity_type": "drone",
|
"timestamp": timestamps.map(lambda value: value.strftime("%Y-%m-%dT%H:%M:%SZ")),
|
"longitude": longitude.astype(float),
|
"latitude": latitude.astype(float),
|
"flight_phase": source[EXCEL_COLUMNS["flight_phase"]].astype(str),
|
"absolute_height_m": pd.to_numeric(source[EXCEL_COLUMNS["absolute_height_m"]], errors="coerce"),
|
"relative_height_m": pd.to_numeric(source[EXCEL_COLUMNS["relative_height_m"]], errors="coerce"),
|
"source_row": range(2, len(source) + 2),
|
}
|
).sort_values("timestamp", kind="stable")
|
duplicates = int(observations.duplicated(["track_id", "timestamp"], keep="last").sum())
|
return observations, {
|
"source_row_count": int(len(source)),
|
"prepared_row_count": int(len(observations)),
|
"duplicate_timestamps_for_analyzer": duplicates,
|
"mission_id": mission_id,
|
"track_id": track_id,
|
"input_time_interpretation": f"Naive timestamps interpreted as {LOCAL_TIMEZONE}, then converted to UTC.",
|
"start_time_utc": observations["timestamp"].iloc[0],
|
"end_time_utc": observations["timestamp"].iloc[-1],
|
}
|
|
|
def extract_route(path: Path, track_id: str) -> tuple[dict[str, object], LineString]:
|
with ZipFile(path) as archive:
|
try:
|
payload = archive.read("wpmz/waylines.wpml")
|
except KeyError as exc:
|
raise ValueError("KMZ does not contain wpmz/waylines.wpml.") from exc
|
root = element_tree.fromstring(payload)
|
namespaces = {"k": KML_NAMESPACE, "w": WPML_NAMESPACE}
|
waypoints: list[tuple[int, list[float]]] = []
|
for placemark in root.findall(".//k:Placemark", namespaces):
|
coordinates = placemark.findtext("k:Point/k:coordinates", namespaces=namespaces)
|
index = placemark.findtext("w:index", namespaces=namespaces)
|
if coordinates is None or index is None:
|
continue
|
values = coordinates.strip().split(",")
|
longitude, latitude = float(values[0]), float(values[1])
|
if not (-180 <= longitude <= 180 and -90 <= latitude <= 90):
|
raise ValueError("KMZ waypoint is outside WGS84 bounds.")
|
waypoints.append((int(index), [longitude, latitude]))
|
waypoints.sort(key=lambda item: item[0])
|
coordinates = [item[1] for item in waypoints]
|
if len(coordinates) < 2:
|
raise ValueError("KMZ must contain at least two ordered waypoints.")
|
line = LineString(coordinates)
|
feature = {
|
"type": "Feature",
|
"properties": {
|
"track_id": track_id,
|
"source_format": "DJI WPMZ waylines.wpml",
|
"waypoint_count": len(coordinates),
|
},
|
"geometry": mapping(line),
|
}
|
return feature, line
|
|
|
def build_zones(
|
path: Path,
|
context_wgs84: LineString,
|
metric_crs: str,
|
search_radius_m: float,
|
) -> tuple[list[dict[str, object]], dict[str, object]]:
|
source = gpd.read_file(path)
|
if source.crs is None:
|
raise ValueError("No-fly GeoJSON has no declared CRS; it cannot be prepared safely.")
|
metric = source.to_crs(metric_crs)
|
context = gpd.GeoSeries([context_wgs84], crs=WGS84).to_crs(metric_crs).iloc[0].buffer(search_radius_m)
|
grouped: dict[str, dict[str, object]] = {}
|
skipped = 0
|
selected_source_features = 0
|
for row in metric.itertuples():
|
geometry = row.geometry
|
if geometry.geom_type == "Point":
|
radius_value = pd.to_numeric(getattr(row, "radius", None), errors="coerce")
|
radius = 0.0 if pd.isna(radius_value) else float(radius_value)
|
if radius <= 0:
|
skipped += 1
|
continue
|
effective_geometry = geometry.buffer(radius)
|
elif geometry.geom_type in {"Polygon", "MultiPolygon"}:
|
effective_geometry = geometry
|
else:
|
skipped += 1
|
continue
|
if not effective_geometry.intersects(context):
|
continue
|
selected_source_features += 1
|
area_id = str(getattr(row, "area_id"))
|
entry = grouped.setdefault(
|
area_id,
|
{
|
"geometry": [],
|
"name": str(getattr(row, "name", "")),
|
"city": str(getattr(row, "city", "")),
|
"level": getattr(row, "level", None),
|
"height": getattr(row, "height", None),
|
"source_feature_count": 0,
|
},
|
)
|
entry["geometry"].append(effective_geometry)
|
entry["source_feature_count"] = int(entry["source_feature_count"]) + 1
|
features: list[dict[str, object]] = []
|
for area_id, entry in sorted(grouped.items()):
|
geometry = unary_union(entry["geometry"])
|
geometry_wgs84 = gpd.GeoSeries([geometry], crs=metric_crs).to_crs(WGS84).iloc[0]
|
features.append(
|
{
|
"type": "Feature",
|
"properties": {
|
"zone_id": f"no-fly-{area_id}",
|
"zone_type": "restricted",
|
"source_area_id": area_id,
|
"source_name": entry["name"],
|
"source_city": entry["city"],
|
"source_level": entry["level"],
|
"source_height_m": entry["height"],
|
"source_feature_count": entry["source_feature_count"],
|
"mapping_note": "Source no-fly area mapped to restricted for this technical rule evaluation.",
|
},
|
"geometry": mapping(geometry_wgs84),
|
}
|
)
|
return features, {
|
"source_feature_count": int(len(source)),
|
"selected_source_feature_count": selected_source_features,
|
"selected_zone_count": len(features),
|
"skipped_source_feature_count": skipped,
|
"source_crs": str(source.crs),
|
"zone_search_radius_m": search_radius_m,
|
"mapping": "Polygon areas are retained; point areas are converted to metric-radius buffers; no-fly is mapped to restricted only for this technical rule evaluation.",
|
}
|
|
|
def build_flyable_zones(
|
path: Path,
|
context_wgs84: LineString,
|
metric_crs: str,
|
search_radius_m: float,
|
) -> tuple[list[dict[str, object]], dict[str, object]]:
|
"""Decode the supplied closed-ring UOM coordinate stream near one flight."""
|
values = np.frombuffer(gzip.decompress(path.read_bytes()), dtype="<i4")
|
if values.size < 3:
|
raise ValueError("Flyable-area Gzip is too small to contain a coordinate stream.")
|
ring_count = int(values[-1])
|
coordinate_values = values[:-1]
|
if ring_count <= 0 or coordinate_values.size % 2:
|
raise ValueError("Flyable-area Gzip has an invalid ring count or coordinate alignment.")
|
coordinates = coordinate_values.reshape(-1, 2)
|
context_metric = gpd.GeoSeries([context_wgs84], crs=WGS84).to_crs(metric_crs).iloc[0].buffer(search_radius_m)
|
context = gpd.GeoSeries([context_metric], crs=metric_crs).to_crs(WGS84).iloc[0]
|
minx, miny, maxx, maxy = context.bounds
|
features: list[dict[str, object]] = []
|
start = 0
|
skipped = 0
|
for ring_index in range(ring_count):
|
end = start + 4
|
while end <= len(coordinates) and not np.array_equal(coordinates[end - 1], coordinates[start]):
|
end += 1
|
if end > len(coordinates):
|
raise ValueError(f"Flyable-area ring {ring_index} is not closed within the coordinate stream.")
|
ring = coordinates[start:end].astype(np.float64) / 1e7
|
start = end
|
ring_minx, ring_miny = ring.min(axis=0)
|
ring_maxx, ring_maxy = ring.max(axis=0)
|
if ring_maxx < minx or ring_minx > maxx or ring_maxy < miny or ring_miny > maxy:
|
continue
|
polygon = Polygon(ring)
|
if not polygon.is_valid:
|
polygon = make_valid(polygon)
|
if polygon.is_empty or polygon.geom_type not in {"Polygon", "MultiPolygon"} or not polygon.intersects(context):
|
skipped += 1
|
continue
|
features.append(
|
{
|
"type": "Feature",
|
"properties": {
|
"zone_id": f"flyable-bin-{ring_index:05d}",
|
"zone_type": "flyable",
|
"source_format": "UOM closed-ring coordinate stream (int32 / 1e7)",
|
"source_ring_index": ring_index,
|
"vertex_count": len(ring),
|
},
|
"geometry": mapping(polygon),
|
}
|
)
|
if start > len(coordinates):
|
raise ValueError("Flyable-area coordinate parser exceeded the available stream.")
|
return features, {
|
"source_ring_count": ring_count,
|
"consumed_coordinate_pair_count": start,
|
"selected_zone_count": len(features),
|
"skipped_candidate_count": skipped,
|
"zone_search_radius_m": search_radius_m,
|
"decoder": "Closed rings parsed from little-endian int32 longitude/latitude values scaled by 1e7.",
|
"limitation": "The supplied file contains geometry only; its publication time, authority and operational semantics remain unverified.",
|
}
|
|
|
def write_json(path: Path, payload: dict[str, object]) -> None:
|
path.write_text(json.dumps(payload, ensure_ascii=False, indent=2), encoding="utf-8")
|
|
|
def main() -> int:
|
args = parse_args()
|
raw_dir = args.raw_dir.resolve()
|
output_dir = args.output.resolve()
|
if args.zone_search_radius_m <= 0:
|
print("Preparation failed: --zone-search-radius-m must be greater than zero.", file=sys.stderr)
|
return 2
|
if output_dir.exists() and not args.overwrite:
|
print(f"Preparation failed: output already exists: {output_dir}", file=sys.stderr)
|
return 2
|
try:
|
flight_path = discover_single(raw_dir / "tracks", ".xlsx")
|
route_path = discover_single(raw_dir / "routes", ".kmz")
|
zone_path = discover_single(raw_dir / "areas", ".geojson")
|
flyable_matches = sorted((raw_dir / "areas").glob("*.gzip"))
|
if len(flyable_matches) > 1:
|
raise ValueError(f"Expected zero or one .gzip file in {raw_dir / 'areas'}, found {len(flyable_matches)}.")
|
flyable_path = flyable_matches[0] if flyable_matches else None
|
observations, observation_metadata = load_observations(flight_path)
|
metric_crs = choose_metric_crs(float(observations["longitude"].mean()), float(observations["latitude"].mean()))
|
route_feature, route = extract_route(route_path, str(observation_metadata["track_id"]))
|
track_line = LineString(observations[["longitude", "latitude"]].to_numpy().tolist())
|
context = unary_union([track_line, route])
|
zones, zone_metadata = build_zones(zone_path, context, metric_crs, args.zone_search_radius_m)
|
flyable_zones: list[dict[str, object]] = []
|
flyable_metadata: dict[str, object] | None = None
|
if flyable_path is not None:
|
flyable_zones, flyable_metadata = build_flyable_zones(
|
flyable_path, context, metric_crs, args.zone_search_radius_m
|
)
|
if not zones:
|
raise ValueError("No source no-fly areas intersect the configured flight context buffer.")
|
observations_metric = gpd.GeoSeries(
|
gpd.points_from_xy(observations["longitude"], observations["latitude"]), crs=WGS84
|
).to_crs(metric_crs)
|
route_metric = gpd.GeoSeries([route], crs=WGS84).to_crs(metric_crs).iloc[0]
|
output_dir.mkdir(parents=True, exist_ok=args.overwrite)
|
observations.to_csv(output_dir / "observations.csv", index=False, encoding="utf-8-sig")
|
write_json(
|
output_dir / "reference_routes.geojson",
|
{
|
"type": "FeatureCollection",
|
"crs": {"type": "name", "properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}},
|
"features": [route_feature],
|
},
|
)
|
write_json(
|
output_dir / "zones.geojson",
|
{
|
"type": "FeatureCollection",
|
"crs": {"type": "name", "properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}},
|
"features": zones,
|
},
|
)
|
if flyable_path is not None:
|
write_json(
|
output_dir / "flyable_zones.geojson",
|
{
|
"type": "FeatureCollection",
|
"crs": {"type": "name", "properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}},
|
"features": flyable_zones,
|
},
|
)
|
manifest = {
|
"case_id": args.case_id,
|
"description": "One user-provided DJI flight, prepared from an XLSX log, WPMZ route and local no-fly-zone subset.",
|
"crs": WGS84,
|
"observations": "observations.csv",
|
"reference_routes": "reference_routes.geojson",
|
"zones": "zones.geojson",
|
}
|
if flyable_path is not None:
|
manifest["flyable_zones"] = "flyable_zones.geojson"
|
write_json(output_dir / f"{args.case_id}.case.json", manifest)
|
metadata = {
|
"case_id": args.case_id,
|
"prepared_at": pd.Timestamp.now(tz=UTC).isoformat(),
|
"raw_inputs": [
|
{"path": str(path), "sha256": sha256(path)}
|
for path in (flight_path, route_path, zone_path)
|
if path is not None
|
] + ([{"path": str(flyable_path), "sha256": sha256(flyable_path)}] if flyable_path is not None else []),
|
"observation_preparation": observation_metadata,
|
"route_preparation": {
|
"source_format": "DJI WPMZ waylines.wpml",
|
"waypoint_count": route_feature["properties"]["waypoint_count"],
|
"track_to_route_min_m": round(float(observations_metric.distance(route_metric).min()), 3),
|
"track_to_route_mean_m": round(float(observations_metric.distance(route_metric).mean()), 3),
|
"track_to_route_max_m": round(float(observations_metric.distance(route_metric).max()), 3),
|
},
|
"zone_preparation": zone_metadata,
|
"flyable_zone_preparation": flyable_metadata,
|
"crs": {"output": WGS84, "metric": metric_crs},
|
"limitations": [
|
"No-fly source semantics, publication time and operational authority have not been independently verified.",
|
"Spatial intersection with a source no-fly area is a technical result, not a legal or operational violation conclusion.",
|
"The reference route does not encode flight-phase semantics such as takeoff, return-to-home or landing.",
|
],
|
}
|
write_json(output_dir / "preparation_metadata.json", metadata)
|
except (OSError, ValueError, KeyError, element_tree.ParseError, pd.errors.ParserError) as exc:
|
print(f"Preparation failed: {exc}", file=sys.stderr)
|
return 2
|
print(f"Prepared 1 real flight case in {output_dir}")
|
return 0
|
|
|
if __name__ == "__main__":
|
raise SystemExit(main())
|
