commit 0891cf12a5aad987fcc525444eaa16a38d8ba2c0 Author: Olaf Date: Tue Oct 20 14:18:35 2026 +0200 Initial Vibe coded project diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..971e588 --- /dev/null +++ b/.gitignore @@ -0,0 +1,11 @@ +.venv/ +__pycache__/ +.pytest_cache/ +.mypy_cache/ +.ruff_cache/ +*.egg-info/ +dist/ +*.sqlite +*.sqlite-journal +*.svg +.flight-map-cache/ diff --git a/README.md b/README.md new file mode 100644 index 0000000..5695cdd --- /dev/null +++ b/README.md @@ -0,0 +1,284 @@ +# Flight Map + +A Python command-line application that **records OpenSky aircraft positions to +SQLite**, then **independently generates a pen-plotter SVG** from a saved session. +Collection never downloads maps or terrain; plotting never contacts OpenSky. + +## Install + +Python 3.11 or newer: + +```sh +python3 -m venv .venv +source .venv/bin/activate +python -m pip install -e . +flight-map --help +``` + +Alternatively, use `uv venv .venv` and `uv pip install -e .`. + +## 1. Collect flights + +The default map rectangle is in Uithoorn, Netherlands: + +| Corner | Latitude | Longitude | +| --- | --- | --- | +| Southwest | 52.232025 | 4.784498 | +| Northeast | 52.253985 | 4.863525 | + +OpenSky supports OAuth client credentials. Obtain your own credentials through +OpenSky and set them in your shell, **not in source code**: + +```sh +export OPENSKY_CLIENT_ID='your-client-id' +export OPENSKY_CLIENT_SECRET='your-client-secret' +flight-map collect --db flights.sqlite --duration 3600 +``` + +An existing token can instead be supplied through `OPENSKY_ACCESS_TOKEN`. +Anonymous access is attempted when no credentials are supplied, with a warning; +availability and daily quotas depend on OpenSky's current policy. Client +credentials automatically obtain and refresh tokens; a manually provided token +does not refresh. + +Collection prints the session ID immediately. Omit `--duration` to run until +Ctrl-C. Observations and progress are committed after each successful poll. +Ctrl-C preserves the observations and marks the session interrupted. Service +errors stop collection with a nonzero exit status and mark the session failed; +HTTP rate-limit responses include `Retry-After` when available. + +Custom bounds are **longitude, latitude, longitude, latitude**, not latitude first: + +```sh +flight-map collect --db flights.sqlite \ + --bounds 4.784498 52.232025 4.863525 52.253985 \ + --interval 30 --collection-margin 2000 --duration 3600 +flight-map sessions --db flights.sqlite +flight-map collect --db flights.sqlite --resume SESSION_ID --duration 1800 +``` + +Reversed corners are normalized. The requested map rectangle, expanded collection +rectangle, interval, source, session times, status, and observation associations +are saved in the database. Resume uses those saved settings; it does not change +the map area or interval. After an abrupt process termination, the last committed +checkpoint remains available even if the status still says `recording`. +Do not run two collectors against the same session simultaneously. + +The default poll interval is 30 seconds, with a minimum of 10 seconds. **A valid +interval does not guarantee enough daily API credits**: select an interval and +duration appropriate to your account's current quota. No historical observations +from before collection starts are fetched. + +## Purge a recording + +To permanently delete one session, including its saved map bounds and all flight +observations: + +```sh +flight-map sessions --db flights.sqlite +flight-map purge --db flights.sqlite --session SESSION_ID --yes +``` + +Stop any collector writing to that session first. The exact session ID and +`--yes` confirmation flag are required. Deletion is transactional: either both +observations and session metadata are removed, or neither is. Other sessions, +map cache entries, and previously exported SVGs are unchanged. Missing sessions +or databases produce an error; a mistyped database path does not create a file. +Deletion cannot be undone without a backup. This is logical deletion, not secure +erasure or database-file compaction. + +## 2. Plot a recorded session + +```sh +flight-map plot --db flights.sqlite --session SESSION_ID --output map.svg +``` + +The map area is read from the session, not from current application defaults. +Flight data is opened read-only. Supply a shorter recording interval if needed: + +```sh +flight-map plot --db flights.sqlite --session SESSION_ID --output morning.svg \ + --from 2026-10-20T08:00:00Z --to 2026-10-20T08:30:00Z +``` + +Timestamps require an explicit time zone; labels are UTC. Both endpoints must be +within the saved recording. Only observations within the selected interval are +used. At least two usable consecutive airborne observations are required per +track. Observations are separated into tracks at callsign changes, ground +observations, and gaps over `--max-gap` seconds (default 120); unknown missing +positions are not bridged beyond that threshold. Callsigns are not verified +passenger-facing flight numbers. Missing callsigns are labeled explicitly, with +the aircraft ICAO identifier retained. + +Flight paths use smooth **centripetal Catmull-Rom interpolation**, exported as +cubic Bezier SVG paths passing through every recorded position. Two-point tracks +remain straight; duplicate positions are handled without invalid coordinates. +Smoothing never joins separate tracks or fills long recording gaps. Curves are +interpolations, not evidence of the aircraft's precise motion between samples. +The actual curves are clipped geometrically to the output rectangle, and their +in-map portions determine the time-weighted height calculation. + +### Map cache and offline plotting + +The first plot downloads OpenStreetMap **vector features**, not raster tiles, +using Overpass. Roads, waterways, buildings, and water outlines are rendered as +unfilled lines. Minor pedestrian/cycle paths, steps, proposed roads, and roads +under construction are omitted from the drawing. Building outlines whose shorter +dimension would be smaller than three pen widths are omitted to prevent dense +housing from turning into ink blobs; larger paper or a finer pen retains more +buildings without another download. Multipolygon relations are not assembled in this initial version; +only returned way geometry is drawn. + +Downloaded map data is cached in `.flight-map-cache/`. The cache key includes +the service URL, exact map bounds, and feature query. Different time ranges, +page sizes, and flight selections for the same map reuse that data, including +when offline. Cache entries do not expire automatically: + +```sh +flight-map plot --db flights.sqlite --session SESSION_ID --output map.svg \ + --cache /path/to/shared-map-cache +flight-map plot --db flights.sqlite --session SESSION_ID --output updated.svg \ + --cache /path/to/shared-map-cache --refresh-map +``` + +Refresh validates the response before atomically replacing the cached entry. +A failed refresh leaves the previous entry intact and reports an error. A damaged +entry is reported rather than silently ignored; use `--refresh-map` to replace it. +Cached responses have no embedded credentials. + +For entirely offline plotting, pass a local Overpass JSON export containing +`elements` with way `geometry`, produced by an `out geom;` query: + +```sh +flight-map plot --db flights.sqlite --session SESSION_ID --output offline.svg \ + --map-data local-map.json +``` + +### Height above the map's average ground level + +**Without suitable elevation inputs the map still plots, but heights are `N/A` +and the program emits warnings.** It never substitutes pressure altitude, zero +ground elevation, or sea-level altitude. + +The calculation is: + +```text +height = time-weighted mean aircraft altitude inside the saved map bounds + minus area-weighted mean ground elevation inside those same bounds +``` + +Outside-map extensions do not contribute. Heights are in metres, not terrain +clearance under the aircraft. + +Supply a terrain GeoTIFF whose first band contains ground elevation in metres, +has a CRS, and covers the entire map rectangle. Its declared vertical reference +must be known. Sampling uses an equal-area 64 by 64 geographic grid and honors +nodata, raster band scale, and offset. Any missing sampled ground value makes +the mean unavailable. Choose terrain resolution appropriate to the map; the +grid is an approximation, not a survey-grade terrain integral. + +For terrain already referenced to the WGS84 ellipsoid: + +```sh +flight-map plot --db flights.sqlite --session SESSION_ID --output heights.svg \ + --terrain terrain-ellipsoid.tif --terrain-reference wgs84-ellipsoid \ + --confirm-geometric-ellipsoid +``` + +For orthometric terrain, also supply a geoid separation GeoTIFF in metres, +representing **N = h - H**, matching the terrain's vertical datum and covering +the track observations used for conversion: + +```sh +flight-map plot --db flights.sqlite --session SESSION_ID --output heights.svg \ + --terrain terrain.tif --terrain-reference orthometric \ + --geoid matching-geoid.tif --confirm-geometric-ellipsoid +``` + +The application computes `H = h - N` for each geometric aircraft altitude. +OpenSky documents its field as geometric altitude but does not establish every +observation's datum or accuracy. The confirmation flag is a deliberate assertion +by the operator that treating these heights as WGS84 ellipsoidal is appropriate. +Do not set it without verifying your source's interpretation. A missing geometric +altitude or geoid value on a contributing segment makes that track's height +unavailable. Pressure altitude is recorded for provenance, not used for this +calculation. For Dutch terrain in NAP, use a conversion consistent with NAP; +do not assume an arbitrary global geoid matches it. + +### Pen-plotter settings + +Defaults: A4 landscape, 12 mm page margin, 0.3 mm pen, 1.8 mm lettering, and up to +4 mm of recorded flight extension beyond the map border. + +```sh +flight-map plot --db flights.sqlite --session SESSION_ID --output large.svg \ + --width 420 --height 297 --pen-width 0.3 --font-height 2 --extension 5 +``` + +SVG dimensions are in millimetres. A local equirectangular projection preserves +the rectangular geographic bounds and uses a common scale for map and tracks; +it is intended for local maps, not large geographic extents. Geometry is +physically clipped, not hidden with masks. Groups `map`, `flights`, and +`lettering` can be selected separately. Text, including attribution, is drawn +with an original single-stroke alphabet: no installed fonts, raster images, +filled glyphs, or SVG clipping support are required. + +Each track has a numeric identifier and a legend entry with UTC dates/times, +callsign, ICAO identifier, and average height. For dense recordings, select +shorter periods or larger paper. The program fails explicitly if the legend +does not fit or a track identifier cannot be placed legibly. Flight identifiers +avoid each other; intersections with map geometry or flight paths are not fully +optimized. Map lines are simplified at 0.05 mm tolerance and identical paths +deduplicated; physical pen-path travel optimization is left to plotter software. + +Existing output is protected unless `--force` is supplied. Map bounds crossing +the antimeridian, widths above 180 degrees, and latitudes outside [-85, 85] are +rejected. + +## Offline synthetic demo + +This creates **invented** tracks, map features, terrain, and geoid data, for +checking the complete workflow without network access: + +```sh +python examples/make_demo.py demo +``` + +Run the plotting command printed by the script. Its heights use synthetic, +known-compatible references; these inputs are not a real Uithoorn dataset. + +## Development + +```sh +python -m pip install -e '.[dev]' +python -m pytest +ruff check flight_map tests examples +ruff format --check flight_map tests examples +mypy flight_map +``` + +Tests cover recording and resume settings, persistence, deduplication, bounds, +projection, track gaps, clipping, time-weighted heights, missing terrain, vertical +references, cache reuse/refresh/failure, and repeatable offline SVG generation +without mutating the database. + +## Data policies and attribution + +- Map: [OpenStreetMap contributors](https://www.openstreetmap.org/copyright), + licensed under ODbL. Respect the chosen Overpass service's + [usage guidance](https://dev.overpass-api.de/overpass-doc/en/preface/commons.html); + cache and avoid unnecessary requests. +- Flights: [The OpenSky Network](https://opensky-network.org). This initial + application is intended for non-commercial use, subject to its current + [terms](https://opensky-network.org/about/terms-of-use) and + [API documentation](https://openskynetwork.github.io/opensky-api/rest.html). + Coverage and fields vary; collection does not guarantee complete tracks. + Cite Matthias Schaefer, Martin Strohmeier, Vincent Lenders, Ivan Martinovic, + and Matthias Wilhelm, "Bringing Up OpenSky: A Large-scale ADS-B Sensor Network + for Research", IPSN 2014, pp. 83-94, when publishing derived work. +- Terrain and geoid: supplied by the operator; follow their separate licenses + and attribution requirements when distributing outputs. + +Map and flight attribution is included as plotter paths, with further provenance +in SVG metadata. This is an artistic recording/plotting tool, not a navigation +or aviation safety product. diff --git a/REQUIREMENTS.md b/REQUIREMENTS.md new file mode 100644 index 0000000..3d5ddbe --- /dev/null +++ b/REQUIREMENTS.md @@ -0,0 +1,188 @@ + +# Flight Map Requirements + +## Purpose + +Create a Python program that generates an SVG map of a rectangular geographic +area with recorded aircraft flight tracks overlaid, suitable for drawing with +a physical pen plotter. + +## Separate collection and plotting workflow + +- Provide separate command-line commands for flight-data collection and SVG + generation. Run collection first; generate maps later from the saved database. +- The collection command accepts optional opposite map corners and the SQLite + database destination. If corners are omitted, use the default map area + including Uithoorn in the Netherlands. It records flight data without + generating an SVG or requiring map or terrain downloads. +- Store a recording-session record containing a session identifier, the + normalized requested map bounds, the expanded collection bounds, the flight + data source, the collection interval, and the recording start and end times. + Preserve this metadata across interruption or restart. +- Associate every stored flight observation with its recording session so its + map area and collection context remain identifiable. +- The plotting command accepts the database, a recording-session selection, + an optional time interval within that recording, and output settings. + Use the selected session's stored map bounds rather than requiring the user + to enter the map area again. +- Plotting must not contact OpenSky or start flight-data collection. It may + obtain OpenStreetMap-derived vector data and terrain elevation data + separately as needed. +- Allow multiple maps with different time selections and output settings to be + generated from the same recording without modifying the recorded flight + observations or session metadata. +- Report missing sessions, invalid time selections, and absence of usable + observations explicitly. +- Provide a purge command requiring an exact session identifier and explicit + confirmation to permanently delete that session's observations and metadata + in one transaction. Preserve all other sessions and their observations. +- Purging a missing session or database must report an error without creating + a database. Stop collection for the selected session before purging it. + +## Map input and output + +- At collection time, accept the longitude and latitude of two opposite corners + of the map, validate coordinates, and normalize their order. +- Provide a predefined default rectangular map area that includes the city of + Uithoorn in the Netherlands. Explicitly supplied corners override this default. +- Default southwest corner: latitude 52.232025, longitude 4.784498. + Default northeast corner: latitude 52.253985, longitude 4.863525. + Command-line bounds use longitude first: 4.784498 52.232025 4.863525 52.253985. +- Store the actual bounds used, including when the default is selected, in the + recording session. Later plotting must use those stored bounds rather than + a potentially changed application default. +- At plotting time, load the map bounds from the selected recording session. +- Obtain OpenStreetMap-derived vector map data, either as vector tiles or + geographic feature data, following the source's usage policy and including + the required attribution. Standard raster map tiles are not suitable for + the pen-plotter output. +- Cache downloaded vector map data persistently by source, geographic bounds, + and feature query. Reuse it across plotting runs without downloading it + again. Provide a configurable cache directory and an explicit refresh option. +- Write cache entries atomically so interrupted downloads or writes do not + replace a previously usable entry. Report damaged cache entries explicitly. +- Generate an SVG using a consistent geographic projection for the map and + flight tracks. +- Handle invalid input and download failures with clear error messages. + +## Pen-plotter output + +- Render map features and flight tracks as vector paths or line geometry, + without embedded raster images or reliance on filled areas, shading, or + other screen-only effects. +- Render all labels and attribution as paths, using single-line stroke + lettering suitable for pen plotting rather than relying on installed fonts + or tracing filled glyph outlines. +- Resolve clipping into the exported geometry so map features stop at the map + boundary and flight-track extensions remain within the output page. Do not + rely on SVG clipping masks being supported by the plotter software. +- Specify physical page dimensions and a matching SVG viewBox so the intended + drawing size is preserved when imported into plotter software. +- Keep map geometry, flight tracks, and lettering in separate named SVG groups + for independent selection and plotting. +- Select and simplify map features to keep the drawing legible at its physical + output size. Omit minor pedestrian and cycle paths and suppress building + outlines whose shorter dimension is smaller than three pen widths. + +## Flight data and local recording + +- Use OpenSky Network as the initial flight-data source. +- Retrieve live aircraft state vectors for the selected geographic area. +- Query an area slightly larger than the map so recorded flight paths can + extend beyond the map boundary. +- Collect positions periodically, respecting OpenSky's current authentication + requirements, rate limits, and usage terms. +- Store observations locally in SQLite so maps can be generated later without + fetching those observations again. +- Store the aircraft's ICAO identifier, callsign when available, position + timestamp, longitude, latitude, and available altitude measurements. +- Support generating a map from locally recorded observations for a selected + time interval. +- Do not require historical data from before local recording began. Unrestricted + historical access is not assumed. +- Report recording failures explicitly and handle missing fields and gaps in + coverage without inventing observations. +- Initial use is non-commercial, subject to OpenSky's applicable terms and + attribution requirements. + +## Flight-track rendering + +- Overlay recorded flight paths that cross the requested map area. +- Fit smooth interpolating curves through recorded flight positions rather + than drawing angular polylines. Export flight paths as cubic Bezier geometry + suitable for a pen plotter. Two-point tracks may remain straight. +- Do not smooth across track gaps, on-ground interruptions, or flight changes. + Handle repeated positions without generating invalid curve coordinates. +- Clip the actual curves to the output bounds, not just their sample points, + and use the same curved in-map portions for time and average-height selection. +- Extend the visible tracks a little beyond the map boundary using actual + recorded positions from the surrounding area, not fabricated extrapolation. +- Keep the map background within its requested bounds and provide space in the + SVG for the extended tracks and their labels. +- Identify each flight track with its date, time, callsign or flight number when + available, and average height above the map's average ground elevation. +- Do not present an OpenSky callsign as a verified passenger-facing flight + number. Callsign-to-flight-number lookup is not required for the initial + version. +- Calculate the map's average ground elevation from terrain elevation data + covering the requested map bounds, excluding the surrounding collection area. + Use an area-weighted mean or an equivalent equal-area sampling method. +- Calculate each track's average height as its average aircraft altitude over + the portion inside the requested map bounds minus the map's average ground + elevation. Outside-map extensions must not affect this calculation. +- Use compatible units and a common vertical reference for aircraft altitude + and terrain elevation, applying conversions where needed. Do not subtract + terrain elevation directly from uncorrected barometric pressure altitude. +- Label the result explicitly as average height above the map's average ground + level, with units. This is not height above the terrain directly beneath each + aircraft position and must not be presented as terrain clearance. +- Indicate when the height cannot be calculated because compatible aircraft + altitude or sufficient terrain data is unavailable; do not silently substitute + altitude above sea level or assume zero ground elevation. +- Make labels readable and indicate unavailable metadata explicitly. + +## Implementation and validation + +- Implement the program in Python with a command-line interface. +- Document installation, the separate collection and plotting commands, + session selection, map generation, configuration, data-source limitations, + and attribution. +- Test coordinate handling, projection, local recording, track selection, + average ground elevation, vertical-reference compatibility, average-height + calculation, missing elevation data, and SVG generation. +- Test persistence of session bounds and observation associations, plotting + after collection has stopped, and repeated plotting without changing the + recorded flight data. Verify that collection does not require map or terrain + services and plotting makes no OpenSky requests. +- Test that omitting map corners selects the Uithoorn default, explicit corners + override it, and the selected bounds are persisted for later plotting. +- Validate that the exported SVG contains drawable vector geometry, path-based + lettering, correct physical dimensions, and geometry clipped to its intended + bounds, with no raster images or dependence on text rendering or clipping + masks. +- Include an end-to-end test using fixed sample data without requiring live + network access. +- Test cache reuse, separate entries for different map areas, explicit refresh, + and preservation of cached data when refresh fails. + +## Initial implementation defaults and limitations + +- Use Overpass to obtain OpenStreetMap roads, waterways, building outlines, and + water outlines, with local Overpass JSON also accepted for offline plotting. +- Default to A4 landscape (297 by 210 mm), 12 mm page margins, a 0.3 mm pen, + and track extensions of up to 4 mm beyond the map border. +- Collect every 30 seconds within a 2,000 m margin around the saved map area. +- Display timestamps in UTC using date and 24-hour time, and heights in metres. +- Accept local terrain GeoTIFFs with a declared vertical reference. For + orthometric terrain, accept a matching geoid-separation raster for converting + ellipsoidal aircraft heights. Require explicit confirmation of the geometric + altitude interpretation before calculating heights; otherwise label them + unavailable. Never silently substitute barometric altitude. +- Estimate average ground elevation using a 64 by 64 equal-area sampling grid. + Any missing terrain sample makes the average unavailable. +- Use time-weighted aircraft altitude with linear interpolation between recorded + timestamps, restricted to the smoothed path's portion inside the map. Use + centripetal Catmull-Rom interpolation converted to cubic Bezier curves for + flight paths. Split tracks at changes + in callsign, on-ground observations, or gaps exceeding 120 seconds by default. +- Reject antimeridian-crossing bounds and latitudes outside [-85, 85] explicitly. \ No newline at end of file diff --git a/examples/make_demo.py b/examples/make_demo.py new file mode 100644 index 0000000..d9b3b3a --- /dev/null +++ b/examples/make_demo.py @@ -0,0 +1,104 @@ +"""Create entirely synthetic offline demo inputs; no real map or flight observations.""" + +import json +import sys +from pathlib import Path + +import numpy as np +import rasterio +from rasterio.transform import from_bounds + +from flight_map.geo import DEFAULT_BOUNDS +from flight_map.storage import ( + Observation, + checkpoint, + create_session, + database, + initialize, + save_observations, +) + + +def main() -> None: + if len(sys.argv) != 2: + raise SystemExit("Usage: python examples/make_demo.py OUTPUT_DIRECTORY") + folder = Path(sys.argv[1]) + folder.mkdir(parents=True, exist_ok=True) + db = folder / "demo.sqlite" + if db.exists(): + raise SystemExit(f"Refusing to replace {db}") + started = 1_792_497_600.0 + with database(db) as connection: + initialize(connection) + session = create_session(connection, DEFAULT_BOUNDS, 2000, 30, started) + items = [] + northings = [52.235, 52.237, 52.244, 52.247, 52.246] + for index in range(5): + items.append( + Observation( + "abc123", + "DEMO123", + started + index * 30, + 4.77 + index * 0.027, + northings[index], + 900 + index * 30, + 1043 + index * 30, + ) + ) + items.append( + Observation( + "def456", + "DEMO456", + started + index * 30, + 4.85 - index * 0.013, + 52.26 - index * 0.008, + 700, + 843, + ) + ) + save_observations(connection, session, items, started + 150) + checkpoint(connection, session.id, started + 150, "stopped") + ways = [ + [(4.77, 52.235), (4.81, 52.24), (4.87, 52.25)], + [(4.80, 52.23), (4.82, 52.245), (4.81, 52.26)], + [(4.835, 52.23), (4.84, 52.26)], + [(4.82, 52.242), (4.823, 52.242), (4.823, 52.244), (4.82, 52.244), (4.82, 52.242)], + ] + (folder / "map.json").write_text( + json.dumps( + { + "elements": [ + { + "type": "way", + "id": index, + "geometry": [{"lon": lon, "lat": lat} for lon, lat in points], + } + for index, points in enumerate(ways) + ] + } + ) + ) + for name, value in (("terrain.tif", -2), ("geoid.tif", 43)): + with rasterio.open( + folder / name, + "w", + driver="GTiff", + height=64, + width=64, + count=1, + dtype="float32", + crs="EPSG:4326", + transform=from_bounds(*session.collection_bounds.values, 64, 64), + ) as dataset: + dataset.write(np.full((1, 64, 64), value, dtype="float32")) + print(f"Synthetic demo session: {session.id}") + print( + f"flight-map plot --db '{db}' --session {session.id} " + f"--map-data '{folder / 'map.json'}' --terrain '{folder / 'terrain.tif'}' " + f"--terrain-reference orthometric --geoid '{folder / 'geoid.tif'}' " + f"--confirm-geometric-ellipsoid --output '{folder / 'demo.svg'}'" + ) + + +if __name__ == "__main__": + main() diff --git a/flight_map/__init__.py b/flight_map/__init__.py new file mode 100644 index 0000000..abe7ce6 --- /dev/null +++ b/flight_map/__init__.py @@ -0,0 +1 @@ +"""OpenSky recording and independent pen-plotter map generation.""" diff --git a/flight_map/__main__.py b/flight_map/__main__.py new file mode 100644 index 0000000..45504aa --- /dev/null +++ b/flight_map/__main__.py @@ -0,0 +1,3 @@ +from flight_map.cli import main + +raise SystemExit(main()) diff --git a/flight_map/cli.py b/flight_map/cli.py new file mode 100644 index 0000000..95e8a76 --- /dev/null +++ b/flight_map/cli.py @@ -0,0 +1,194 @@ +import argparse +import logging +import math +import sqlite3 +import sys +import time +from datetime import UTC, datetime +from pathlib import Path + +from flight_map.geo import DEFAULT_BOUNDS, Bounds +from flight_map.storage import create_session, database, get_session, initialize, purge_session + + +def positive(value: str) -> float: + try: + result = float(value) + except ValueError as exc: + raise argparse.ArgumentTypeError("Expected a positive finite number.") from exc + if not math.isfinite(result) or result <= 0: + raise argparse.ArgumentTypeError("Expected a positive finite number.") + return result + + +def timestamp(value: str) -> float: + try: + parsed = datetime.fromisoformat(value) + except ValueError as exc: + raise argparse.ArgumentTypeError( + "Use an ISO 8601 timestamp, e.g. 2026-10-20T12:00:00Z." + ) from exc + if parsed.tzinfo is None: + raise argparse.ArgumentTypeError("Timestamps must include a time zone, such as Z for UTC.") + return parsed.timestamp() + + +def parser() -> argparse.ArgumentParser: + root = argparse.ArgumentParser(description="Record flights, then draw pen-plotter SVG maps.") + sub = root.add_subparsers(dest="command", required=True) + record = sub.add_parser("collect", help="Record OpenSky positions to SQLite; no map downloads") + record.add_argument("--db", type=Path, required=True) + record.add_argument( + "--bounds", + type=float, + nargs=4, + metavar=("LON1", "LAT1", "LON2", "LAT2"), + help="Opposite corners, longitude first; default: Uithoorn", + ) + record.add_argument( + "--interval", type=positive, help="Seconds between polls (default: 30, min: 10)" + ) + record.add_argument( + "--collection-margin", type=positive, help="Metres outside map (default: 2000)" + ) + record.add_argument( + "--duration", type=positive, help="Seconds to record; omit to run until Ctrl-C" + ) + record.add_argument( + "--resume", metavar="SESSION", help="Resume a stopped/interrupted recording" + ) + listing = sub.add_parser("sessions", help="List saved recording sessions") + listing.add_argument("--db", type=Path, required=True) + purge = sub.add_parser("purge", help="Permanently delete one session and its observations") + purge.add_argument("--db", type=Path, required=True) + purge.add_argument("--session", required=True, help="Exact recording session ID to delete") + purge.add_argument( + "--yes", + action="store_true", + required=True, + help="Confirm permanent deletion; stop collection for this session first", + ) + plot = sub.add_parser("plot", help="Draw from a saved session; never contacts OpenSky") + plot.add_argument("--db", type=Path, required=True) + plot.add_argument("--session", required=True) + plot.add_argument("--output", type=Path, required=True) + plot.add_argument("--from", dest="start", type=timestamp) + plot.add_argument("--to", dest="end", type=timestamp) + plot.add_argument("--width", type=positive, default=297, help="Page width in mm") + plot.add_argument("--height", type=positive, default=210, help="Page height in mm") + plot.add_argument("--margin", type=positive, default=12, help="Page margin in mm") + plot.add_argument("--extension", type=positive, default=4, help="Track extension in mm") + plot.add_argument("--pen-width", type=positive, default=0.3, help="Pen width in mm") + plot.add_argument("--font-height", type=positive, default=1.8, help="Letter height in mm") + plot.add_argument( + "--max-gap", type=positive, default=120, help="Split tracks at gaps above seconds" + ) + plot.add_argument("--cache", type=Path, default=Path(".flight-map-cache")) + plot.add_argument( + "--refresh-map", + action="store_true", + help="Download fresh vector map data and atomically replace its cache entry", + ) + plot.add_argument( + "--map-data", type=Path, help="Local Overpass JSON (out geom), for offline plotting" + ) + plot.add_argument( + "--terrain", type=Path, help="Local terrain GeoTIFF; elevations must be in metres" + ) + plot.add_argument( + "--terrain-reference", + choices=("unknown", "wgs84-ellipsoid", "orthometric"), + default="unknown", + ) + plot.add_argument( + "--geoid", + type=Path, + help="Metre-valued h-H raster matching the terrain's orthometric datum", + ) + plot.add_argument( + "--confirm-geometric-ellipsoid", + action="store_true", + help="Confirm OpenSky geometric heights can be treated as WGS84 ellipsoidal", + ) + plot.add_argument("--force", action="store_true", help="Replace an existing SVG") + return root + + +def main(argv: list[str] | None = None) -> int: + args = parser().parse_args(argv) + logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s") + try: + if args.command == "collect": + from flight_map.record import OpenSky, collect + + if args.resume and (args.bounds or args.interval or args.collection_margin): + raise ValueError("Resuming uses stored bounds and interval; do not override them.") + client = OpenSky() + if args.resume and not args.db.exists(): + raise ValueError("Cannot resume: recording database does not exist.") + with database(args.db) as connection: + initialize(connection) + if args.resume: + session = get_session(connection, args.resume) + else: + bounds = Bounds.corners(*args.bounds) if args.bounds else DEFAULT_BOUNDS + session = create_session( + connection, + bounds, + args.collection_margin or 2000, + args.interval or 30, + time.time(), + ) + print(f"Session: {session.id}", flush=True) + print(f"Map bounds (W,S,E,N): {session.bounds.values}", flush=True) + collect(connection, session, client, args.duration) + elif args.command == "sessions": + with database(args.db, readonly=True) as connection: + for row in connection.execute("SELECT id FROM sessions ORDER BY started"): + session = get_session(connection, row["id"]) + count = connection.execute( + "SELECT COUNT(*) FROM observations WHERE session_id = ?", (session.id,) + ).fetchone()[0] + started = datetime.fromtimestamp(session.started, UTC).isoformat() + ended = datetime.fromtimestamp(session.ended, UTC).isoformat() + print( + f"{session.id} {session.status} {started} to {ended} " + f"{count} observations bounds={session.bounds.values}" + ) + elif args.command == "purge": + with database(args.db, must_exist=True) as connection: + count = purge_session(connection, args.session) + print(f"Purged session {args.session} and {count} observation(s).") + else: + from flight_map.plot import PlotOptions, render + + options = PlotOptions( + args.width, + args.height, + args.margin, + args.extension, + args.pen_width, + args.font_height, + args.max_gap, + ) + count = render( + args.db, + args.session, + args.output, + start=args.start, + end=args.end, + options=options, + cache=args.cache, + local_map=args.map_data, + terrain=args.terrain, + terrain_reference=args.terrain_reference, + geoid=args.geoid, + confirm_geometric_ellipsoid=args.confirm_geometric_ellipsoid, + force=args.force, + refresh_map=args.refresh_map, + ) + print(f"Wrote {args.output} with {count} flight track(s).") + except (ValueError, RuntimeError, OSError, sqlite3.Error) as exc: + print(f"Error: {exc}", file=sys.stderr) + return 1 + return 0 diff --git a/flight_map/curves.py b/flight_map/curves.py new file mode 100644 index 0000000..41f3405 --- /dev/null +++ b/flight_map/curves.py @@ -0,0 +1,158 @@ +"""Interpolating cubic flight paths and geometric (not mask-based) clipping.""" + +import math +from dataclasses import dataclass + +Point2D = tuple[float, float] +Rectangle = tuple[float, float, float, float] + + +def mix(a: Point2D, b: Point2D, fraction: float) -> Point2D: + return a[0] + (b[0] - a[0]) * fraction, a[1] + (b[1] - a[1]) * fraction + + +@dataclass(frozen=True) +class Cubic: + start: Point2D + control1: Point2D + control2: Point2D + end: Point2D + + @property + def points(self) -> tuple[Point2D, Point2D, Point2D, Point2D]: + return self.start, self.control1, self.control2, self.end + + def point(self, fraction: float) -> Point2D: + a, b, c = ( + mix(p, q, fraction) + for p, q in ( + (self.start, self.control1), + (self.control1, self.control2), + (self.control2, self.end), + ) + ) + return mix(mix(a, b, fraction), mix(b, c, fraction), fraction) + + def split(self, fraction: float) -> tuple["Cubic", "Cubic"]: + a, b, c = ( + mix(p, q, fraction) + for p, q in ( + (self.start, self.control1), + (self.control1, self.control2), + (self.control2, self.end), + ) + ) + d, e = mix(a, b, fraction), mix(b, c, fraction) + middle = mix(d, e, fraction) + return Cubic(self.start, a, d, middle), Cubic(middle, e, c, self.end) + + def portion(self, first: float, last: float) -> "Cubic": + if not 0 <= first < last <= 1: + raise ValueError("Curve portion must be an increasing interval in [0, 1].") + left = self if last == 1 else self.split(last)[0] + return left if first == 0 else left.split(first / last)[1] + + def transform(self, matrix: list[float]) -> "Cubic": + a, b, d, e, xoff, yoff = matrix + points = [(a * x + b * y + xoff, d * x + e * y + yoff) for x, y in self.points] + return Cubic(*points) + + def boundary_roots(self, axis: int, boundary: float) -> list[float]: + p0, p1, p2, p3 = (point[axis] for point in self.points) + a = -p0 + 3 * p1 - 3 * p2 + p3 + b = 3 * p0 - 6 * p1 + 3 * p2 + c = -3 * p0 + 3 * p1 + critical = [0.0, 1.0] + if abs(a) < 1e-12: + if abs(b) >= 1e-12: + critical.append(-c / (2 * b)) + else: + discriminant = 4 * b * b - 12 * a * c + if discriminant >= 0: + root = math.sqrt(discriminant) + critical.extend(((-2 * b - root) / (6 * a), (-2 * b + root) / (6 * a))) + critical = sorted(t for t in critical if 0 <= t <= 1) + roots = [] + for first, last in zip(critical, critical[1:], strict=False): + low, high = first, last + f_low = self.point(low)[axis] - boundary + f_high = self.point(high)[axis] - boundary + if abs(f_low) < 1e-9: + roots.append(low) + if abs(f_high) < 1e-9: + roots.append(high) + if f_low * f_high >= 0: + continue + for _ in range(55): + middle = (low + high) / 2 + f_middle = self.point(middle)[axis] - boundary + if f_low * f_middle <= 0: + high = middle + else: + low, f_low = middle, f_middle + roots.append((low + high) / 2) + return roots + + def inside_intervals(self, rectangle: Rectangle) -> list[tuple[float, float]]: + west, south, east, north = rectangle + xs, ys = [p[0] for p in self.points], [p[1] for p in self.points] + if min(xs) >= west and max(xs) <= east and min(ys) >= south and max(ys) <= north: + return [(0.0, 1.0)] + if max(xs) < west or min(xs) > east or max(ys) < south or min(ys) > north: + return [] + parameters = [0.0, 1.0] + for axis, boundary in ((0, west), (0, east), (1, south), (1, north)): + parameters.extend(self.boundary_roots(axis, boundary)) + parameters.sort() + unique = [parameters[0]] + for value in parameters[1:]: + if value - unique[-1] > 1e-12: + unique.append(value) + result: list[tuple[float, float]] = [] + for first, last in zip(unique, unique[1:], strict=False): + x, y = self.point((first + last) / 2) + if west - 1e-9 <= x <= east + 1e-9 and south - 1e-9 <= y <= north + 1e-9: + if result and abs(result[-1][1] - first) < 1e-12: + result[-1] = result[-1][0], last + else: + result.append((first, last)) + return result + + def clipped(self, rectangle: Rectangle) -> list["Cubic"]: + return [self.portion(first, last) for first, last in self.inside_intervals(rectangle)] + + +def interpolate(points: list[Point2D]) -> list[Cubic]: + """Centripetal Catmull-Rom converted to cubic Beziers through every sample.""" + result = [] + for index in range(len(points) - 1): + p1, p2 = points[index], points[index + 1] + if p1 == p2: + result.append(Cubic(p1, p1, p2, p2)) + continue + p0 = points[index - 1] if index > 0 else p1 + p3 = points[index + 2] if index + 2 < len(points) else p2 + if p0 == p1: + p0 = (2 * p1[0] - p2[0], 2 * p1[1] - p2[1]) + if p3 == p2: + p3 = (2 * p2[0] - p1[0], 2 * p2[1] - p1[1]) + d01 = math.sqrt(math.dist(p0, p1)) + d12 = math.sqrt(math.dist(p1, p2)) + d23 = math.sqrt(math.dist(p2, p3)) + control1 = [] + control2 = [] + for axis in (0, 1): + tangent1 = d12 * ( + (p1[axis] - p0[axis]) / d01 + - (p2[axis] - p0[axis]) / (d01 + d12) + + (p2[axis] - p1[axis]) / d12 + ) + tangent2 = d12 * ( + (p2[axis] - p1[axis]) / d12 + - (p3[axis] - p1[axis]) / (d12 + d23) + + (p3[axis] - p2[axis]) / d23 + ) + control1.append(p1[axis] + tangent1 / 3) + control2.append(p2[axis] - tangent2 / 3) + result.append(Cubic(p1, (control1[0], control1[1]), (control2[0], control2[1]), p2)) + return result diff --git a/flight_map/geo.py b/flight_map/geo.py new file mode 100644 index 0000000..a6c28d8 --- /dev/null +++ b/flight_map/geo.py @@ -0,0 +1,63 @@ +import math +from dataclasses import dataclass + +from shapely.geometry import Polygon, box + +EARTH_RADIUS = 6_371_008.8 + + +@dataclass(frozen=True) +class Bounds: + west: float + south: float + east: float + north: float + + def __post_init__(self) -> None: + if not all(math.isfinite(v) for v in self.values): + raise ValueError("Map coordinates must be finite.") + if not (-180 <= self.west < self.east <= 180): + raise ValueError("Longitude bounds must increase within [-180, 180].") + if self.east - self.west > 180: + raise ValueError("Antimeridian-crossing maps are not supported.") + if not (-85 <= self.south < self.north <= 85): + raise ValueError("Latitude bounds must increase within [-85, 85].") + + @property + def values(self) -> tuple[float, float, float, float]: + return self.west, self.south, self.east, self.north + + @classmethod + def corners(cls, lon1: float, lat1: float, lon2: float, lat2: float) -> "Bounds": + return cls(min(lon1, lon2), min(lat1, lat2), max(lon1, lon2), max(lat1, lat2)) + + def expand(self, metres: float) -> "Bounds": + if not math.isfinite(metres) or metres <= 0: + raise ValueError("Collection margin must be a positive number of metres.") + dy = math.degrees(metres / EARTH_RADIUS) + dx = dy / math.cos(math.radians((self.south + self.north) / 2)) + return Bounds(self.west - dx, self.south - dy, self.east + dx, self.north + dy) + + def contains(self, lon: float, lat: float) -> bool: + return self.west <= lon <= self.east and self.south <= lat <= self.north + + +DEFAULT_BOUNDS = Bounds(4.784, 52.232, 4.864, 52.258) + + +class Projection: + """Local equirectangular projection: geographic bounds remain rectangular.""" + + def __init__(self, bounds: Bounds): + self.lon0 = (bounds.west + bounds.east) / 2 + self.lat0 = (bounds.south + bounds.north) / 2 + self.x_factor = EARTH_RADIUS * math.cos(math.radians(self.lat0)) + west, south = self.point(bounds.west, bounds.south) + east, north = self.point(bounds.east, bounds.north) + self.rectangle: Polygon = box(west, south, east, north) + + def point(self, lon: float, lat: float) -> tuple[float, float]: + return ( + self.x_factor * math.radians(lon - self.lon0), + EARTH_RADIUS * math.radians(lat - self.lat0), + ) diff --git a/flight_map/lettering.py b/flight_map/lettering.py new file mode 100644 index 0000000..358fe7e --- /dev/null +++ b/flight_map/lettering.py @@ -0,0 +1,71 @@ +"""Small original single-stroke alphabet; coordinates use a 4 by 6 grid.""" + +import unicodedata + +from shapely.geometry import LineString + +GLYPHS = { + "A": ["06 20 46", "14 34"], + "B": ["06 00 30 41 32 02", "32 44 35 06"], + "C": ["40 10 01 05 16 46"], + "D": ["06 00 20 41 45 26 06"], + "E": ["40 00 06 46", "03 33"], + "F": ["06 00 40", "03 33"], + "G": ["40 10 01 05 16 46 43 23"], + "H": ["00 06", "40 46", "03 43"], + "I": ["00 40", "20 26", "06 46"], + "J": ["00 40 45 36 16 05"], + "K": ["00 06", "40 03 46"], + "L": ["00 06 46"], + "M": ["06 00 23 40 46"], + "N": ["06 00 46 40"], + "O": ["10 30 41 45 36 16 05 01 10"], + "P": ["06 00 30 41 32 02"], + "Q": ["10 30 41 45 36 16 05 01 10", "24 46"], + "R": ["06 00 30 41 32 02", "22 46"], + "S": ["40 10 01 12 33 44 35 06"], + "T": ["00 40", "20 26"], + "U": ["00 05 16 36 45 40"], + "V": ["00 26 40"], + "W": ["00 16 23 36 40"], + "X": ["00 46", "40 06"], + "Y": ["00 23 40", "23 26"], + "Z": ["00 40 06 46"], + "0": ["10 30 41 45 36 16 05 01 10", "05 41"], + "1": ["02 20 26", "06 46"], + "2": ["01 10 30 41 32 06 46"], + "3": ["00 30 41 22 43 45 36 06"], + "4": ["00 03 43", "40 46"], + "5": ["40 00 03 33 44 35 06"], + "6": ["40 10 01 05 16 36 45 34 03"], + "7": ["00 40 16"], + "8": ["12 01 10 30 41 32 12 03 05 16 36 45 43 32"], + "9": ["43 12 01 10 30 41 45 36 06"], + "-": ["03 43"], + ".": ["25 26"], + ",": ["25 16"], + ":": ["21 22", "24 25"], + "/": ["06 40"], + "(": ["30 11 15 36"], + ")": ["10 31 35 16"], + "?": ["01 10 30 41 32 23 24", "25 26"], + "+": ["03 43", "20 26"], + "=": ["02 42", "04 44"], + "'": ["20 21"], + " ": [], +} + + +def lettering(text: str, x: float, y: float, height: float) -> list[LineString]: + text = unicodedata.normalize("NFKD", text).upper() + result = [] + for index, character in enumerate(text): + if character not in GLYPHS: + raise ValueError(f"Unsupported stroke lettering character: {character!r}") + for stroke in GLYPHS[character]: + points = [ + (x + index * height + int(p[0]) * height / 6, y + int(p[1]) * height / 6) + for p in stroke.split() + ] + result.append(LineString(points)) + return result diff --git a/flight_map/maps.py b/flight_map/maps.py new file mode 100644 index 0000000..d679853 --- /dev/null +++ b/flight_map/maps.py @@ -0,0 +1,127 @@ +import hashlib +import json +import logging +import math +import os +import tempfile +from pathlib import Path +from typing import Any +from urllib.parse import urlencode + +from shapely.geometry import LineString + +from flight_map.geo import Bounds, Projection +from flight_map.network import request_json + +OVERPASS_URL = "https://overpass-api.de/api/interpreter" +LOGGER = logging.getLogger(__name__) + + +def validate_map_data(data: Any) -> dict[str, Any]: + if not isinstance(data, dict) or not isinstance(data.get("elements"), list): + raise ValueError("Map data must be an Overpass JSON object containing elements.") + if data.get("remark"): + raise ValueError(f"Overpass could not complete the map query: {data['remark']}") + return data + + +def map_data( + bounds: Bounds, cache: Path, local_file: Path | None = None, *, refresh: bool = False +) -> dict[str, Any]: + if local_file and refresh: + raise ValueError("--refresh-map cannot be combined with --map-data.") + bbox = f"{bounds.south},{bounds.west},{bounds.north},{bounds.east}" + query = ( + "[out:json][timeout:60];(" + f'way["highway"]({bbox});way["waterway"]({bbox});' + f'way["building"]({bbox});way["natural"="water"]({bbox});' + ");out geom;" + ) + key = hashlib.sha256((OVERPASS_URL + "\n" + query).encode()).hexdigest() + path = local_file or cache / f"osm-{key}.json" + if path.exists() and (local_file or not refresh): + LOGGER.info("Reading %s map data: %s", "local" if local_file else "cached", path) + with path.open() as stream: + try: + data = json.load(stream) + except (ValueError, UnicodeError) as exc: + raise ValueError( + f"Invalid map JSON in {path}; correct the file or use --refresh-map " + "to replace a damaged cache entry." + ) from exc + elif local_file: + raise ValueError(f"Local map data file does not exist: {local_file}") + else: + LOGGER.info("Downloading vector map data from Overpass for %s", bounds.values) + data = request_json( + OVERPASS_URL, + data=urlencode({"data": query}).encode(), + headers={"Content-Type": "application/x-www-form-urlencoded"}, + ) + validate_map_data(data) + cache.mkdir(parents=True, exist_ok=True) + temporary: Path | None = None + try: + with tempfile.NamedTemporaryFile( + mode="w", encoding="utf-8", dir=cache, suffix=".tmp", delete=False + ) as stream: + temporary = Path(stream.name) + json.dump(data, stream) + stream.flush() + os.fsync(stream.fileno()) + os.replace(temporary, path) + finally: + if temporary and temporary.exists(): + temporary.unlink() + LOGGER.info("Cached vector map data: %s", path) + return validate_map_data(data) + + +def map_lines( + data: dict[str, Any], projection: Projection, *, min_building_size: float = 0 +) -> list[LineString]: + result = [] + for element in data["elements"]: + if not isinstance(element, dict): + raise ValueError("Map data contains a malformed feature.") + if element.get("type") != "way": + continue + tags = element.get("tags", {}) + if not isinstance(tags, dict): + raise ValueError("Map feature tags must be an object.") + if tags.get("highway") in { + "footway", + "path", + "steps", + "cycleway", + "bridleway", + "corridor", + "construction", + "proposed", + }: + continue + geometry = element.get("geometry") + if not isinstance(geometry, list) or len(geometry) < 2: + raise ValueError("Map way has no usable geometry; export using 'out geom'.") + try: + positions = [(float(p["lon"]), float(p["lat"])) for p in geometry] + except (KeyError, ValueError, TypeError) as exc: + raise ValueError("Map way contains invalid geographic coordinates.") from exc + if any( + not math.isfinite(lon) + or not math.isfinite(lat) + or not -180 <= lon <= 180 + or not -90 <= lat <= 90 + for lon, lat in positions + ): + raise ValueError("Map way coordinates must be finite and within geographic limits.") + points = [projection.point(lon, lat) for lon, lat in positions] + line = LineString(points) + if "building" in tags: + west, south, east, north = line.bounds + if min(east - west, north - south) < min_building_size: + continue + result.append(line) + if not result: + raise ValueError("No drawable map features were found for this area.") + return result diff --git a/flight_map/network.py b/flight_map/network.py new file mode 100644 index 0000000..038c689 --- /dev/null +++ b/flight_map/network.py @@ -0,0 +1,25 @@ +import json +from typing import Any +from urllib.error import HTTPError, URLError +from urllib.request import Request, urlopen + +USER_AGENT = "FlightMap/0.1 (local non-commercial flight mapping)" + + +def request_json( + url: str, *, data: bytes | None = None, headers: dict[str, str] | None = None +) -> Any: + request = Request(url, data=data, headers={"User-Agent": USER_AGENT, **(headers or {})}) + try: + with urlopen(request, timeout=60) as response: + return json.load(response) + except HTTPError as exc: + retry = exc.headers.get("Retry-After") + suffix = f"; retry after {retry}" if retry else "" + raise RuntimeError(f"Remote service returned HTTP {exc.code}{suffix}.") from exc + except (URLError, TimeoutError) as exc: + raise RuntimeError( + f"Remote service connection failed: {exc.reason if isinstance(exc, URLError) else exc}" + ) from exc + except (ValueError, UnicodeError) as exc: + raise RuntimeError("Remote service returned invalid JSON.") from exc diff --git a/flight_map/plot.py b/flight_map/plot.py new file mode 100644 index 0000000..f2d9ae8 --- /dev/null +++ b/flight_map/plot.py @@ -0,0 +1,315 @@ +import logging +import math +import os +import tempfile +from collections.abc import Iterator +from dataclasses import dataclass +from datetime import UTC, datetime +from pathlib import Path +from xml.etree import ElementTree as ET + +from shapely.affinity import affine_transform +from shapely.geometry import LineString, box +from shapely.geometry.base import BaseGeometry + +from flight_map.curves import Cubic +from flight_map.geo import Projection +from flight_map.lettering import lettering +from flight_map.maps import map_data, map_lines +from flight_map.storage import database, get_session, load_observations +from flight_map.terrain import Elevation +from flight_map.tracks import make_tracks + +LOGGER = logging.getLogger(__name__) +SVG_NS = "http://www.w3.org/2000/svg" +ET.register_namespace("", SVG_NS) + + +def tag(name: str) -> str: + return f"{{{SVG_NS}}}{name}" + + +def lines(geometry: BaseGeometry) -> Iterator[LineString]: + if isinstance(geometry, LineString): + if geometry.length > 0: + yield geometry + elif hasattr(geometry, "geoms"): + for part in geometry.geoms: + yield from lines(part) + + +def add_line(group: ET.Element, line: LineString) -> None: + coordinates = list(line.coords) + path = "M " + " L ".join(f"{x:.5f},{y:.5f}" for x, y in coordinates) + ET.SubElement(group, tag("path"), {"d": path}) + + +def add_curves(group: ET.Element, curves: list[Cubic]) -> None: + commands = [] + previous = None + for curve in curves: + if all(point == curve.start for point in curve.points): + continue + if previous is None or math.dist(previous, curve.start) > 1e-7: + commands.append(f"M {curve.start[0]:.5f},{curve.start[1]:.5f}") + commands.append( + "C " + + " ".join( + f"{x:.5f},{y:.5f}" + for x, y in ( + curve.control1, + curve.control2, + curve.end, + ) + ) + ) + previous = curve.end + if commands: + ET.SubElement(group, tag("path"), {"d": " ".join(commands)}) + + +def label(group: ET.Element, text: str, x: float, y: float, height: float) -> None: + for line in lettering(text, x, y, height): + add_line(group, line) + + +def utc(timestamp: float) -> str: + return datetime.fromtimestamp(timestamp, UTC).strftime("%Y-%m-%d %H:%M:%S") + + +@dataclass(frozen=True) +class PlotOptions: + width: float = 297 + height: float = 210 + margin: float = 12 + extension: float = 4 + pen_width: float = 0.3 + font_height: float = 1.8 + max_gap: float = 120 + + def validate(self) -> None: + if not all( + math.isfinite(value) and value > 0 + for value in ( + self.width, + self.height, + self.margin, + self.extension, + self.pen_width, + self.font_height, + self.max_gap, + ) + ): + raise ValueError( + "Page dimensions, margins, pen width and plot settings must be positive." + ) + if 2 * (self.margin + self.extension) >= min(self.width, self.height): + raise ValueError("Margins and extensions leave no room for the map.") + + +def render( + db: Path, + session_id: str, + output: Path, + *, + start: float | None = None, + end: float | None = None, + options: PlotOptions | None = None, + cache: Path = Path(".flight-map-cache"), + local_map: Path | None = None, + terrain: Path | None = None, + terrain_reference: str = "unknown", + geoid: Path | None = None, + confirm_geometric_ellipsoid: bool = False, + force: bool = False, + refresh_map: bool = False, +) -> int: + options = options or PlotOptions() + options.validate() + if output.exists() and not force: + raise ValueError(f"Output already exists: {output}. Use --force to replace it.") + if any( + output.resolve() == source.resolve() + for source in (db, local_map, terrain, geoid) + if source is not None + ): + raise ValueError("SVG output must not overwrite the recording database or input data.") + with database(db, readonly=True) as connection: + session = get_session(connection, session_id) + start = session.started if start is None else start + end = session.ended if end is None else end + if ( + not math.isfinite(start) + or not math.isfinite(end) + or start < session.started + or end > session.ended + or start >= end + ): + raise ValueError("Time selection must be a nonempty interval within the recording.") + observations = load_observations(connection, session_id, start, end) + if not observations: + raise ValueError("No recorded observations in the selected interval.") + projection = Projection(session.bounds) + with Elevation(terrain, terrain_reference, geoid, confirm_geometric_ellipsoid) as elevation: + ground = elevation.mean_ground(session.bounds) + tracks = make_tracks(observations, projection, elevation, ground, options.max_gap) + if not tracks: + raise ValueError("No usable airborne tracks cross the map in the selected interval.") + if ground is None: + LOGGER.warning( + "Average ground elevation unavailable: supply a complete terrain raster " + "and its vertical reference." + ) + missing = sum(track.average_height is None for track in tracks) + if missing: + LOGGER.warning( + "%d track height(s) unavailable: compatible geometric altitude, terrain " + "and vertical-reference conversion are required.", + missing, + ) + data = map_data(session.bounds, cache, local_map, refresh=refresh_map) + row_height = options.font_height * 1.8 + legend_height = (len(tracks) + 4) * row_height + available_width = options.width - 2 * (options.margin + options.extension) + available_height = options.height - 2 * (options.margin + options.extension) - legend_height + if available_height < 25: + raise ValueError("Too many tracks for this page. Select a shorter interval or larger page.") + west, south, east, north = projection.rectangle.bounds + scale = min(available_width / (east - west), available_height / (north - south)) + features = map_lines(data, projection, min_building_size=3 * options.pen_width / scale) + map_width, map_height = (east - west) * scale, (north - south) * scale + left = (options.width - map_width) / 2 + top = options.margin + options.extension + matrix = [scale, 0, 0, -scale, left - west * scale, top + north * scale] + map_rectangle = box(left, top, left + map_width, top + map_height) + track_rectangle = box( + left - options.extension, + top - options.extension, + left + map_width + options.extension, + top + map_height + options.extension, + ) + root = ET.Element( + tag("svg"), + { + "width": f"{options.width:g}mm", + "height": f"{options.height:g}mm", + "viewBox": f"0 0 {options.width:g} {options.height:g}", + "fill": "none", + "stroke": "black", + "stroke-width": f"{options.pen_width:g}", + "stroke-linecap": "round", + "stroke-linejoin": "round", + }, + ) + ET.SubElement(root, tag("title")).text = "Recorded flight map" + ET.SubElement(root, tag("desc")).text = ( + f"Session {session.id}; bounds W,S,E,N={session.bounds.values}; " + f"{utc(start)} to {utc(end)} UTC. Local equirectangular projection. " + "Height is time-weighted aircraft altitude within the map minus the " + "area-sampled mean ground elevation, not terrain clearance. " + f"Terrain reference: {terrain_reference}; " + f"mean ground elevation in metres: {ground}; " + f"geometric ellipsoid interpretation confirmed: {confirm_geometric_ellipsoid}. " + "Map data: OpenStreetMap contributors, ODbL, https://www.openstreetmap.org/copyright. " + "Flight data: The OpenSky Network, https://opensky-network.org. " + "Schaefer et al., Bringing Up OpenSky: A Large-scale ADS-B Sensor Network " + "for Research, IPSN 2014, pp. 83-94." + ) + map_group = ET.SubElement(root, tag("g"), {"id": "map"}) + flight_group = ET.SubElement(root, tag("g"), {"id": "flights"}) + text_group = ET.SubElement(root, tag("g"), {"id": "lettering"}) + add_line(map_group, LineString(map_rectangle.exterior.coords)) + seen = set() + for feature in features: + page_feature = affine_transform(feature, matrix).intersection(map_rectangle) + for part in lines(page_feature): + simplified = part.simplify(0.05) + key = simplified.normalize().wkb + if key not in seen: + for simplified_part in lines(simplified): + add_line(map_group, simplified_part) + seen.add(key) + occupied: list[BaseGeometry] = [] + for index, track in enumerate(tracks, 1): + page_curves = [ + part + for curve in track.curves + for part in curve.transform(matrix).clipped(track_rectangle.bounds) + ] + group = ET.SubElement( + flight_group, + tag("g"), + { + "id": f"flight-{index}", + "data-icao24": track.observations[0].icao24, + "data-callsign": track.observations[0].callsign or "unavailable", + "data-start-utc": utc(track.start), + "data-end-utc": utc(track.end), + "data-average-height-metres": ( + f"{track.average_height:.3f}" + if track.average_height is not None + else "unavailable" + ), + }, + ) + add_curves(group, page_curves) + short_label = str(index) + placed = False + for fraction in (0.9, 0.1, 0.7, 0.3, 0.5): + curve_index = min(int(fraction * len(page_curves)), len(page_curves) - 1) + anchor_x, anchor_y = page_curves[curve_index].point( + fraction * len(page_curves) - curve_index + ) + for dx, dy in ((1, 1), (1, -3), (-4, 1), (-4, -3)): + x, y = anchor_x + dx, anchor_y + dy + extent = box( + x - 0.3, + y - 0.3, + x + len(short_label) * options.font_height, + y + options.font_height + 0.3, + ) + if track_rectangle.covers(extent) and not any( + extent.intersects(previous) for previous in occupied + ): + label(text_group, short_label, x, y, options.font_height) + occupied.append(extent) + placed = True + break + if placed: + break + if not placed: + raise ValueError("Cannot place track identifiers legibly; select fewer tracks.") + legend_y = top + map_height + options.extension + row_height + ground_label = "N/A" if ground is None else f"{ground:.1f} M" + legend = [ + "FLIGHT TIMES UTC - HEIGHT ABOVE MAP AVERAGE GROUND, NOT TERRAIN CLEARANCE", + f"MAP AVERAGE GROUND ELEVATION: {ground_label} - N/A MEANS UNAVAILABLE", + ] + for index, track in enumerate(tracks, 1): + first = track.observations[0] + height_label = "N/A" if track.average_height is None else f"{track.average_height:.0f} M" + legend.append( + f"{index}: {utc(track.start)} TO {utc(track.end)} " + f"{first.callsign or 'CALLSIGN N/A'} ({first.icao24.upper()}) HEIGHT {height_label}" + ) + legend.extend( + [ + "MAP: OPENSTREETMAP CONTRIBUTORS (ODBL) - OPENSTREETMAP.ORG/COPYRIGHT", + "FLIGHTS: THE OPENSKY NETWORK - OPENSKY-NETWORK.ORG", + ] + ) + for index, text in enumerate(legend): + if len(text) * options.font_height > options.width - 2 * options.margin: + raise ValueError("Legend is too wide; reduce --font-height or increase page width.") + label(text_group, text, options.margin, legend_y + index * row_height, options.font_height) + output.parent.mkdir(parents=True, exist_ok=True) + temp_path: Path | None = None + try: + with tempfile.NamedTemporaryFile(dir=output.parent, suffix=".svg", delete=False) as stream: + temp_path = Path(stream.name) + ET.ElementTree(root).write(stream, encoding="utf-8", xml_declaration=True) + os.replace(temp_path, output) + finally: + if temp_path and temp_path.exists(): + temp_path.unlink() + return len(tracks) diff --git a/flight_map/record.py b/flight_map/record.py new file mode 100644 index 0000000..651e585 --- /dev/null +++ b/flight_map/record.py @@ -0,0 +1,116 @@ +import logging +import os +import sqlite3 +import time +from urllib.parse import urlencode + +from flight_map.geo import Bounds +from flight_map.network import request_json +from flight_map.storage import Observation, Session, checkpoint, save_observations + +LOGGER = logging.getLogger(__name__) +TOKEN_URL = ( + "https://auth.opensky-network.org/auth/realms/opensky-network/protocol/openid-connect/token" +) + + +class OpenSky: + def __init__(self) -> None: + self.client_id = os.environ.get("OPENSKY_CLIENT_ID") + self.secret = os.environ.get("OPENSKY_CLIENT_SECRET") + self.token = os.environ.get("OPENSKY_ACCESS_TOKEN") + self.expires = float("inf") if self.token else 0.0 + if bool(self.client_id) != bool(self.secret): + raise ValueError("Set both OPENSKY_CLIENT_ID and OPENSKY_CLIENT_SECRET.") + if not self.client_id and not self.token: + LOGGER.warning("Using anonymous OpenSky access; daily credit limits apply.") + + def observations(self, bounds: Bounds) -> list[Observation]: + if self.client_id and time.time() >= self.expires: + response = request_json( + TOKEN_URL, + data=urlencode( + { + "grant_type": "client_credentials", + "client_id": self.client_id, + "client_secret": self.secret, + } + ).encode(), + headers={"Content-Type": "application/x-www-form-urlencoded"}, + ) + if not isinstance(response, dict) or not isinstance(response.get("access_token"), str): + raise RuntimeError("OpenSky authentication returned no access token.") + self.token = response["access_token"] + self.expires = time.time() + float(response.get("expires_in", 300)) - 30 + headers = {"Authorization": f"Bearer {self.token}"} if self.token else {} + query = urlencode( + { + "lomin": bounds.west, + "lamin": bounds.south, + "lomax": bounds.east, + "lamax": bounds.north, + } + ) + response = request_json( + "https://opensky-network.org/api/states/all?" + query, headers=headers + ) + if not isinstance(response, dict) or "states" not in response: + raise RuntimeError("OpenSky response is missing aircraft states.") + states = response["states"] + if states is None: + return [] + if not isinstance(states, list): + raise RuntimeError("OpenSky aircraft states are malformed.") + observations = [] + for state in states: + if not isinstance(state, list) or len(state) < 14: + raise RuntimeError("OpenSky returned a malformed aircraft state.") + if any(state[index] is None for index in (3, 5, 6)): + continue + try: + observations.append( + Observation( + str(state[0]).lower(), + str(state[1]).strip() or None if state[1] else None, + float(state[3]), + float(state[5]), + float(state[6]), + float(state[7]) if state[7] is not None else None, + float(state[13]) if state[13] is not None else None, + bool(state[8]), + ) + ) + except (ValueError, TypeError) as exc: + raise RuntimeError("OpenSky returned invalid aircraft observation values.") from exc + return observations + + +def collect( + connection: sqlite3.Connection, + session: Session, + client: OpenSky, + duration: float | None = None, +) -> None: + started = time.monotonic() + checkpoint(connection, session.id, time.time(), "recording") + try: + while True: + tick = time.monotonic() + observations = client.observations(session.collection_bounds) + save_observations(connection, session, observations, time.time()) + LOGGER.info("Session %s: received %d observations", session.id, len(observations)) + remaining = None if duration is None else duration - (time.monotonic() - started) + if remaining is not None and remaining <= 0: + break + delay = max(0, session.interval - (time.monotonic() - tick)) + time.sleep(delay if remaining is None else min(delay, remaining)) + if duration is not None and time.monotonic() - started >= duration: + break + except KeyboardInterrupt: + checkpoint(connection, session.id, time.time(), "interrupted") + LOGGER.info("Recording interrupted; committed observations are preserved.") + return + except (RuntimeError, ValueError, OSError, sqlite3.Error): + checkpoint(connection, session.id, time.time(), "failed") + raise + checkpoint(connection, session.id, time.time(), "stopped") diff --git a/flight_map/storage.py b/flight_map/storage.py new file mode 100644 index 0000000..07e9322 --- /dev/null +++ b/flight_map/storage.py @@ -0,0 +1,204 @@ +import json +import math +import sqlite3 +import uuid +from collections.abc import Iterator +from contextlib import contextmanager +from dataclasses import dataclass +from pathlib import Path + +from flight_map.geo import Bounds + + +@dataclass(frozen=True) +class Observation: + icao24: str + callsign: str | None + timestamp: float + longitude: float + latitude: float + baro_altitude: float | None + geo_altitude: float | None + on_ground: bool = False + + def __post_init__(self) -> None: + if len(self.icao24) != 6 or any(c not in "0123456789abcdef" for c in self.icao24): + raise ValueError("Invalid ICAO24 aircraft identifier.") + if not all(math.isfinite(v) for v in (self.timestamp, self.longitude, self.latitude)): + raise ValueError("Observation contains non-finite coordinates or timestamp.") + if not (-180 <= self.longitude <= 180 and -90 <= self.latitude <= 90): + raise ValueError("Observation coordinates are outside geographic limits.") + for altitude in (self.baro_altitude, self.geo_altitude): + if altitude is not None and not math.isfinite(altitude): + raise ValueError("Altitude must be finite or missing.") + + +@dataclass(frozen=True) +class Session: + id: str + bounds: Bounds + collection_bounds: Bounds + interval: float + started: float + ended: float + status: str + source: str + + +@contextmanager +def database( + path: Path, *, readonly: bool = False, must_exist: bool = False +) -> Iterator[sqlite3.Connection]: + if readonly or must_exist: + mode = "ro" if readonly else "rw" + connection = sqlite3.connect(path.resolve().as_uri() + f"?mode={mode}", uri=True) + else: + connection = sqlite3.connect(path) + connection.row_factory = sqlite3.Row + connection.execute("PRAGMA foreign_keys = ON") + try: + yield connection + finally: + connection.close() + + +def initialize(connection: sqlite3.Connection) -> None: + version = connection.execute("PRAGMA user_version").fetchone()[0] + if version not in (0, 1): + raise ValueError(f"Unsupported flight database version: {version}") + connection.executescript(""" + CREATE TABLE IF NOT EXISTS sessions ( + id TEXT PRIMARY KEY, bounds TEXT NOT NULL, collection_bounds TEXT NOT NULL, + interval REAL NOT NULL, started REAL NOT NULL, ended REAL NOT NULL, + status TEXT NOT NULL, source TEXT NOT NULL + ); + CREATE TABLE IF NOT EXISTS observations ( + session_id TEXT NOT NULL REFERENCES sessions(id), + icao24 TEXT NOT NULL, callsign TEXT, timestamp REAL NOT NULL, + longitude REAL NOT NULL, latitude REAL NOT NULL, + baro_altitude REAL, geo_altitude REAL, on_ground INTEGER NOT NULL, + PRIMARY KEY (session_id, icao24, timestamp) + ); + CREATE INDEX IF NOT EXISTS observation_time + ON observations(session_id, timestamp); + PRAGMA user_version = 1; + """) + + +def create_session( + connection: sqlite3.Connection, bounds: Bounds, margin: float, interval: float, now: float +) -> Session: + if not math.isfinite(interval) or interval < 10: + raise ValueError("Recording interval must be at least 10 seconds.") + session = Session( + uuid.uuid4().hex, + bounds, + bounds.expand(margin), + interval, + now, + now, + "recording", + "OpenSky Network", + ) + with connection: + connection.execute( + "INSERT INTO sessions VALUES (?, ?, ?, ?, ?, ?, ?, ?)", + ( + session.id, + json.dumps(bounds.values), + json.dumps(session.collection_bounds.values), + interval, + now, + now, + session.status, + session.source, + ), + ) + return session + + +def get_session(connection: sqlite3.Connection, session_id: str) -> Session: + row = connection.execute("SELECT * FROM sessions WHERE id = ?", (session_id,)).fetchone() + if row is None: + raise ValueError(f"Recording session not found: {session_id}") + return Session( + row["id"], + Bounds(*json.loads(row["bounds"])), + Bounds(*json.loads(row["collection_bounds"])), + row["interval"], + row["started"], + row["ended"], + row["status"], + row["source"], + ) + + +def checkpoint(connection: sqlite3.Connection, session_id: str, now: float, status: str) -> None: + with connection: + connection.execute( + "UPDATE sessions SET ended = ?, status = ? WHERE id = ?", (now, status, session_id) + ) + + +def purge_session(connection: sqlite3.Connection, session_id: str) -> int: + """Delete one session and its observations in a single write transaction.""" + with connection: + connection.execute("BEGIN IMMEDIATE") + exists = connection.execute("SELECT 1 FROM sessions WHERE id = ?", (session_id,)).fetchone() + if exists is None: + raise ValueError(f"Recording session not found: {session_id}") + deleted = connection.execute( + "DELETE FROM observations WHERE session_id = ?", (session_id,) + ).rowcount + connection.execute("DELETE FROM sessions WHERE id = ?", (session_id,)) + return deleted + + +def save_observations( + connection: sqlite3.Connection, session: Session, observations: list[Observation], now: float +) -> None: + with connection: + for item in observations: + if ( + item.timestamp < session.started + or item.timestamp > now + or not session.collection_bounds.contains(item.longitude, item.latitude) + ): + continue + connection.execute( + "INSERT OR IGNORE INTO observations VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?)", + ( + session.id, + item.icao24, + item.callsign, + item.timestamp, + item.longitude, + item.latitude, + item.baro_altitude, + item.geo_altitude, + item.on_ground, + ), + ) + connection.execute("UPDATE sessions SET ended = ? WHERE id = ?", (now, session.id)) + + +def load_observations( + connection: sqlite3.Connection, session_id: str, start: float, end: float +) -> list[Observation]: + return [ + Observation( + row["icao24"], + row["callsign"], + row["timestamp"], + row["longitude"], + row["latitude"], + row["baro_altitude"], + row["geo_altitude"], + bool(row["on_ground"]), + ) + for row in connection.execute( + """SELECT * FROM observations WHERE session_id = ? AND timestamp BETWEEN ? AND ? + ORDER BY icao24, timestamp""", + (session_id, start, end), + ) + ] diff --git a/flight_map/terrain.py b/flight_map/terrain.py new file mode 100644 index 0000000..c5932d8 --- /dev/null +++ b/flight_map/terrain.py @@ -0,0 +1,102 @@ +import math +from contextlib import ExitStack +from pathlib import Path +from typing import Any + +import rasterio +from pyproj import Transformer + +from flight_map.geo import Bounds +from flight_map.storage import Observation + + +class Elevation: + """Metre-valued terrain and optional h-H geoid separation rasters.""" + + def __init__( + self, + terrain: Path | None, + reference: str = "unknown", + geoid: Path | None = None, + confirm_geometric_ellipsoid: bool = False, + ) -> None: + self.terrain_path = terrain + self.reference = reference + self.geoid_path = geoid + self.confirm_geometric_ellipsoid = confirm_geometric_ellipsoid + self.stack = ExitStack() + self.ground: Any = None + self.geoid: Any = None + self.ground_transform: Transformer | None = None + self.geoid_transform: Transformer | None = None + + def __enter__(self) -> "Elevation": + with ExitStack() as stack: + if self.terrain_path: + self.ground = stack.enter_context(rasterio.open(self.terrain_path)) + self.ground_transform = self._transform(self.ground) + if self.geoid_path: + self.geoid = stack.enter_context(rasterio.open(self.geoid_path)) + self.geoid_transform = self._transform(self.geoid) + self.stack = stack.pop_all() + return self + + def __exit__(self, *args: object) -> None: + self.stack.close() + + @staticmethod + def _transform(dataset: Any) -> Transformer: + if dataset.crs is None: + raise ValueError("Elevation raster has no coordinate reference system.") + if dataset.count < 1: + raise ValueError("Elevation raster contains no data band.") + return Transformer.from_crs("EPSG:4326", dataset.crs, always_xy=True) + + @staticmethod + def sample(dataset: Any, transform: Transformer, lon: float, lat: float) -> float | None: + x, y = transform.transform(lon, lat) + row, col = dataset.index(x, y) + if row < 0 or col < 0 or row >= dataset.height or col >= dataset.width: + return None + value = next(dataset.sample([(x, y)], indexes=1, masked=True))[0] + if bool(getattr(value, "mask", False)): + return None + result = float(value) * dataset.scales[0] + dataset.offsets[0] + return result if math.isfinite(result) else None + + def mean_ground(self, bounds: Bounds, samples: int = 64) -> float | None: + if not self.ground or self.ground_transform is None or self.reference == "unknown": + return None + if samples < 2: + raise ValueError("Terrain sampling requires at least two samples per axis.") + sin_south, sin_north = ( + math.sin(math.radians(bounds.south)), + math.sin(math.radians(bounds.north)), + ) + values = [] + # Uniform longitude and sin(latitude) samples represent equal spherical areas. + for row in range(samples): + lat = math.degrees( + math.asin(sin_south + (sin_north - sin_south) * (row + 0.5) / samples) + ) + for col in range(samples): + lon = bounds.west + (bounds.east - bounds.west) * (col + 0.5) / samples + value = self.sample(self.ground, self.ground_transform, lon, lat) + if value is None: + return None + values.append(value) + return sum(values) / len(values) + + def aircraft_altitude(self, observation: Observation) -> float | None: + if not self.confirm_geometric_ellipsoid or observation.geo_altitude is None: + return None + if self.reference == "wgs84-ellipsoid": + return observation.geo_altitude + if self.reference != "orthometric" or not self.geoid or self.geoid_transform is None: + return None + separation = self.sample( + self.geoid, self.geoid_transform, observation.longitude, observation.latitude + ) + if separation is None: + return None + return observation.geo_altitude - separation diff --git a/flight_map/tracks.py b/flight_map/tracks.py new file mode 100644 index 0000000..cfb8396 --- /dev/null +++ b/flight_map/tracks.py @@ -0,0 +1,85 @@ +import math +from dataclasses import dataclass +from itertools import groupby, pairwise + +from flight_map.curves import Cubic, interpolate +from flight_map.geo import Projection +from flight_map.storage import Observation +from flight_map.terrain import Elevation + + +@dataclass +class Track: + observations: list[Observation] + curves: list[Cubic] + average_height: float | None + start: float + end: float + + +def split_tracks(observations: list[Observation], max_gap: float) -> list[list[Observation]]: + if not math.isfinite(max_gap) or max_gap <= 0: + raise ValueError("Maximum track gap must be positive.") + result = [] + for _, aircraft in groupby( + sorted(observations, key=lambda o: (o.icao24, o.timestamp)), key=lambda o: o.icao24 + ): + current: list[Observation] = [] + for item in aircraft: + if item.on_ground: + if len(current) >= 2: + result.append(current) + current = [] + continue + if current and ( + item.timestamp - current[-1].timestamp > max_gap + or item.timestamp <= current[-1].timestamp + or item.callsign != current[-1].callsign + ): + if len(current) >= 2: + result.append(current) + current = [] + current.append(item) + if len(current) >= 2: + result.append(current) + return result + + +def make_tracks( + observations: list[Observation], + projection: Projection, + elevation: Elevation, + ground: float | None, + max_gap: float, +) -> list[Track]: + result = [] + for items in split_tracks(observations, max_gap): + points = [projection.point(o.longitude, o.latitude) for o in items] + curves = interpolate(points) + if all(curve.start == curve.end for curve in curves): + continue + total_time = 0.0 + altitude_integral = 0.0 + missing_height = ground is None + start, end = math.inf, -math.inf + for index, (first, last) in enumerate(pairwise(items)): + duration = last.timestamp - first.timestamp + h0, h1 = elevation.aircraft_altitude(first), elevation.aircraft_altitude(last) + for f0, f1 in curves[index].inside_intervals(projection.rectangle.bounds): + inside_duration = duration * (f1 - f0) + total_time += inside_duration + start = min(start, first.timestamp + f0 * duration) + end = max(end, first.timestamp + f1 * duration) + if h0 is None or h1 is None: + missing_height = True + else: + entry = h0 + (h1 - h0) * f0 + exit_height = h0 + (h1 - h0) * f1 + altitude_integral += (entry + exit_height) / 2 * inside_duration + if total_time == 0: + continue + average = None + if not missing_height and ground is not None: + average = altitude_integral / total_time - ground + result.append(Track(items, curves, average, start, end)) + return sorted(result, key=lambda track: (track.start, track.observations[0].icao24)) diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..bd0900f --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,35 @@ +[build-system] +requires = ["hatchling"] +build-backend = "hatchling.build" + +[project] +name = "flight-map" +version = "0.1.0" +description = "Record OpenSky flights and draw pen-plotter SVG maps" +readme = "README.md" +requires-python = ">=3.11" +dependencies = ["shapely>=2.0,<3", "pyproj>=3.6,<4", "rasterio>=1.4,<2"] + +[project.optional-dependencies] +dev = ["pytest>=8", "ruff>=0.11", "mypy>=1.15", "types-shapely>=2"] + +[project.scripts] +flight-map = "flight_map.cli:main" + +[tool.pytest.ini_options] +testpaths = ["tests"] + +[tool.ruff] +line-length = 100 + +[tool.ruff.lint] +select = ["E4", "E7", "E9", "F", "I", "B", "UP"] + +[tool.mypy] +python_version = "3.11" +check_untyped_defs = true +disallow_untyped_defs = true + +[[tool.mypy.overrides]] +module = ["rasterio", "rasterio.*"] +ignore_missing_imports = true diff --git a/tests/test_curves.py b/tests/test_curves.py new file mode 100644 index 0000000..e651fc0 --- /dev/null +++ b/tests/test_curves.py @@ -0,0 +1,89 @@ +import math +import random + +import pytest + +from flight_map.curves import Cubic, interpolate + + +def test_spline_interpolates_samples_with_continuous_tangent() -> None: + points = [(0.0, 0.0), (1.0, 2.0), (5.0, 3.0), (6.0, -1.0)] + curves = interpolate(points) + assert len(curves) == len(points) - 1 + for index, curve in enumerate(curves): + assert curve.point(0) == pytest.approx(points[index]) + assert curve.point(1) == pytest.approx(points[index + 1]) + for first, second in zip(curves, curves[1:], strict=False): + incoming = (first.end[0] - first.control2[0], first.end[1] - first.control2[1]) + outgoing = (second.control1[0] - second.start[0], second.control1[1] - second.start[1]) + cross = incoming[0] * outgoing[1] - incoming[1] * outgoing[0] + dot = incoming[0] * outgoing[0] + incoming[1] * outgoing[1] + assert abs(cross) < 1e-10 and dot > 0 + # A non-collinear path must actually curve, not just encode straight lines as cubics. + assert curves[0].point(0.5) != pytest.approx((0.5, 1.0)) + + +def test_two_points_remain_a_straight_path_and_duplicates_are_finite() -> None: + curve = interpolate([(0, 0), (9, 6)])[0] + assert curve.point(0.5) == pytest.approx((4.5, 3)) + for segment in interpolate([(0, 0), (0, 0), (1, 2), (1, 2), (3, 0)]): + assert all(math.isfinite(coordinate) for p in segment.points for coordinate in p) + assert interpolate([]) == [] + assert interpolate([(0, 0)]) == [] + + +def test_curve_portion_preserves_original_shape() -> None: + curve = Cubic((-2, 0), (1, 5), (3, -5), (6, 0)) + portion = curve.portion(0.2, 0.8) + for fraction in (0, 0.1, 0.5, 0.8, 1): + assert portion.point(fraction) == pytest.approx(curve.point(0.2 + fraction * 0.6)) + with pytest.raises(ValueError): + curve.portion(0.8, 0.2) + + +def test_curve_can_enter_bounds_when_both_endpoints_are_outside() -> None: + curve = Cubic((-2, 0), (-1, 3), (1, 3), (2, 0)) + rectangle = (-1.0, 1.0, 1.0, 3.0) + portions = curve.clipped(rectangle) + assert len(portions) == 1 + assert portions[0].start[0] == pytest.approx(-1) + assert portions[0].end[0] == pytest.approx(1) + assert portions[0].point(0.5)[1] > 2 + + +def test_clipping_handles_multiple_entries_tangency_and_boundary_lines() -> None: + curve = Cubic((-2, 0), (4, 0), (-4, 0), (2, 0)) + rectangle = (-0.1, -1.0, 0.1, 1.0) + assert len(curve.inside_intervals(rectangle)) == 3 + boundary = Cubic((-2, 0), (-1, 0), (1, 0), (2, 0)) + assert len(boundary.clipped((-1, 0, 1, 1))) == 1 + # Only touches the lower edge; it has no visible interval inside the rectangle. + tangent = Cubic((0, -3), (1, 1), (2, 1), (3, -3)) + assert tangent.clipped((-1, 0, 4, 4)) == [] + + +def test_random_clipping_matches_dense_curve_membership() -> None: + randomizer = random.Random(42) + rectangle = (-1.0, -1.0, 1.0, 1.0) + for _ in range(100): + points = [(randomizer.uniform(-3, 3), randomizer.uniform(-3, 3)) for _ in range(4)] + curve = Cubic(*points) + intervals = curve.inside_intervals(rectangle) + for index in range(201): + t = index / 200 + x, y = curve.point(t) + inside = -1 <= x <= 1 and -1 <= y <= 1 + selected = any(first - 1e-10 <= t <= last + 1e-10 for first, last in intervals) + assert selected == inside + for portion in curve.clipped(rectangle): + for index in range(51): + x, y = portion.point(index / 50) + assert -1 - 1e-8 <= x <= 1 + 1e-8 and -1 - 1e-8 <= y <= 1 + 1e-8 + + +def test_transformation_preserves_bezier_shape() -> None: + curve = Cubic((0, 0), (1, 2), (2, -1), (3, 0)) + transformed = curve.transform([2, 0, 0, -2, 10, 20]) + for t in (0, 0.1, 0.6, 1): + x, y = curve.point(t) + assert transformed.point(t) == pytest.approx((2 * x + 10, -2 * y + 20)) diff --git a/tests/test_flight_map.py b/tests/test_flight_map.py new file mode 100644 index 0000000..76cbd4b --- /dev/null +++ b/tests/test_flight_map.py @@ -0,0 +1,563 @@ +import hashlib +import json +import re +import sqlite3 +from dataclasses import replace +from pathlib import Path +from typing import Any +from unittest.mock import Mock +from xml.etree import ElementTree as ET + +import numpy as np +import pytest +import rasterio +from rasterio.transform import from_bounds + +from flight_map.cli import main, parser, timestamp +from flight_map.geo import DEFAULT_BOUNDS, Bounds, Projection +from flight_map.lettering import GLYPHS, lettering +from flight_map.maps import map_data, map_lines +from flight_map.network import request_json +from flight_map.plot import SVG_NS, PlotOptions, render +from flight_map.record import OpenSky, collect +from flight_map.storage import ( + Observation, + Session, + checkpoint, + create_session, + database, + get_session, + initialize, + load_observations, + save_observations, +) +from flight_map.terrain import Elevation +from flight_map.tracks import make_tracks, split_tracks + +T0 = 1_792_497_600.0 + + +def observation(**changes: Any) -> Observation: + return replace(Observation("abc123", "KLM123", T0 + 10, 4.80, 52.24, 900, 1043), **changes) + + +def raster(path: Path, value: float, *, bounds: Bounds = DEFAULT_BOUNDS) -> Path: + with rasterio.open( + path, + "w", + driver="GTiff", + height=10, + width=10, + count=1, + dtype="float32", + crs="EPSG:4326", + transform=from_bounds(*bounds.values, 10, 10), + nodata=-9999, + ) as dataset: + dataset.write(np.full((1, 10, 10), value, dtype="float32")) + return path + + +@pytest.fixture +def recorded(tmp_path: Path) -> tuple[Path, Session]: + path = tmp_path / "flights.sqlite" + with database(path) as connection: + initialize(connection) + session = create_session(connection, DEFAULT_BOUNDS, 2000, 30, T0) + items = [ + observation(longitude=4.77, timestamp=T0 + 10, geo_altitude=843), + observation(longitude=4.82, timestamp=T0 + 40, geo_altitude=1043), + observation(longitude=4.88, timestamp=T0 + 70, geo_altitude=1243), + ] + save_observations(connection, session, items, T0 + 80) + checkpoint(connection, session.id, T0 + 90, "stopped") + return path, session + + +@pytest.fixture +def local_map(tmp_path: Path) -> Path: + path = tmp_path / "map.json" + path.write_text( + json.dumps( + { + "elements": [ + { + "type": "way", + "id": 1, + "geometry": [ + {"lon": 4.77, "lat": 52.235}, + {"lon": 4.90, "lat": 52.245}, + ], + }, + { + "type": "way", + "id": 2, + "geometry": [ + {"lon": 4.80, "lat": 52.22}, + {"lon": 4.83, "lat": 52.26}, + ], + }, + ] + } + ) + ) + return path + + +def test_default_and_reversed_bounds() -> None: + assert DEFAULT_BOUNDS.values == (4.784498, 52.232025, 4.863525, 52.253985) + assert Bounds.corners(4.863525, 52.253985, 4.784498, 52.232025) == DEFAULT_BOUNDS + assert DEFAULT_BOUNDS.contains(4.83, 52.24) + expanded = DEFAULT_BOUNDS.expand(2000) + assert expanded.west < DEFAULT_BOUNDS.west and expanded.north > DEFAULT_BOUNDS.north + + +@pytest.mark.parametrize( + "values", + [ + (1, 1, 1, 2), + (180, 1, 181, 2), + (-179, 1, 179, 2), + (1, -90, 2, 90), + (1, 2, float("nan"), 3), + ], +) +def test_invalid_bounds(values: tuple[float, float, float, float]) -> None: + with pytest.raises(ValueError): + Bounds(*values) + + +def test_projection_corners() -> None: + projection = Projection(DEFAULT_BOUNDS) + southwest = projection.point(DEFAULT_BOUNDS.west, DEFAULT_BOUNDS.south) + northeast = projection.point(DEFAULT_BOUNDS.east, DEFAULT_BOUNDS.north) + assert projection.rectangle.bounds == (*southwest, *northeast) + assert northeast[0] > southwest[0] and northeast[1] > southwest[1] + + +def test_persistence_and_deduplication(recorded: tuple[Path, Session]) -> None: + path, original = recorded + with database(path) as connection: + session = get_session(connection, original.id) + assert session.bounds == DEFAULT_BOUNDS + assert session.collection_bounds == original.collection_bounds + assert session.status == "stopped" + items = load_observations(connection, session.id, T0, T0 + 100) + assert len(items) == 3 + save_observations(connection, session, items, T0 + 100) + assert len(load_observations(connection, session.id, T0, T0 + 100)) == 3 + assert items[0].geo_altitude == 843 + with pytest.raises(ValueError, match="not found"): + get_session(connection, "missing") + + +def test_observations_outside_session_are_not_saved(recorded: tuple[Path, Session]) -> None: + path, session = recorded + with database(path) as connection: + save_observations( + connection, + session, + [ + observation(timestamp=T0 - 1), + observation(timestamp=T0 + 500), + observation(longitude=5.5), + ], + T0 + 100, + ) + assert len(load_observations(connection, session.id, T0, T0 + 100)) == 3 + + +def test_database_readonly(recorded: tuple[Path, Session]) -> None: + path, _ = recorded + with database(path, readonly=True) as connection: + with pytest.raises(sqlite3.OperationalError, match="readonly"): + connection.execute("DELETE FROM observations") + + +def test_terrain_and_datum_conversion(tmp_path: Path) -> None: + terrain = raster(tmp_path / "terrain.tif", -2) + geoid = raster(tmp_path / "geoid.tif", 43) + with Elevation(terrain, "orthometric", geoid, True) as elevation: + assert elevation.mean_ground(DEFAULT_BOUNDS, 4) == -2 + assert elevation.aircraft_altitude(observation()) == 1000 + with Elevation(terrain, "orthometric", geoid, False) as elevation: + assert elevation.aircraft_altitude(observation()) is None + with Elevation(terrain, "orthometric", None, True) as elevation: + assert elevation.aircraft_altitude(observation()) is None + with Elevation(terrain, "wgs84-ellipsoid", None, True) as elevation: + assert elevation.aircraft_altitude(observation()) == 1043 + assert elevation.aircraft_altitude(observation(geo_altitude=None)) is None + + +def test_terrain_missing_data(tmp_path: Path) -> None: + terrain = raster(tmp_path / "terrain.tif", -9999) + with Elevation(terrain, "wgs84-ellipsoid") as elevation: + assert elevation.mean_ground(DEFAULT_BOUNDS, 4) is None + small = Bounds(4.80, 52.24, 4.81, 52.25) + raster(terrain, 0, bounds=small) + with Elevation(terrain, "wgs84-ellipsoid") as elevation: + assert elevation.mean_ground(DEFAULT_BOUNDS, 4) is None + + +def test_equal_area_ground_average(tmp_path: Path) -> None: + terrain = raster(tmp_path / "terrain.tif", 0) + with rasterio.open(terrain, "r+") as dataset: + values = np.zeros((10, 10), dtype="float32") + values[:, 5:] = 10 + dataset.write(values, 1) + with Elevation(terrain, "wgs84-ellipsoid") as elevation: + assert elevation.mean_ground(DEFAULT_BOUNDS, 10) == pytest.approx(5) + + +def test_height_is_time_weighted_and_excludes_extensions() -> None: + bounds = Bounds(0, 0, 2, 1) + observations = [ + observation(longitude=-1, latitude=0.5, timestamp=T0, geo_altitude=0), + observation(longitude=1, latitude=0.5, timestamp=T0 + 20, geo_altitude=200), + observation(longitude=3, latitude=0.5, timestamp=T0 + 60, geo_altitude=600), + ] + with Elevation(None, "wgs84-ellipsoid", None, True) as elevation: + tracks = make_tracks(observations, Projection(bounds), elevation, 10, 120) + # Inside: 10s at mean altitude 150 + 20s at mean altitude 300. + assert len(tracks) == 1 + assert tracks[0].average_height == pytest.approx(240) + assert tracks[0].start == pytest.approx(T0 + 10) + assert tracks[0].end == pytest.approx(T0 + 40) + + +def test_stationary_positions_and_missing_altitude() -> None: + items = [ + observation(), + observation(timestamp=T0 + 25), + observation(timestamp=T0 + 40, longitude=4.82, geo_altitude=None), + ] + with Elevation(None, "wgs84-ellipsoid", None, True) as elevation: + tracks = make_tracks(items, Projection(DEFAULT_BOUNDS), elevation, -2, 120) + assert len(tracks) == 1 and tracks[0].average_height is None + with Elevation(None) as elevation: + assert make_tracks(items[:2], Projection(DEFAULT_BOUNDS), elevation, None, 120) == [] + + +def test_track_splitting_and_ground_filter() -> None: + items = [observation(timestamp=T0 + i * 30) for i in range(6)] + items[2] = replace(items[2], on_ground=True) + items[4] = replace(items[4], callsign="OTHER") + items[5] = replace(items[5], callsign="OTHER") + assert [len(track) for track in split_tracks(items, 120)] == [2, 2] + assert split_tracks([observation(), observation(timestamp=T0 + 500)], 120) == [] + + +def test_stroke_font_has_no_fills_or_missing_common_characters() -> None: + for character in GLYPHS: + assert all(line.length > 0 for line in lettering(character, 0, 0, 2)) + with pytest.raises(ValueError, match="Unsupported"): + lettering("@", 0, 0, 2) + + +def test_end_to_end_offline_repeatable( + recorded: tuple[Path, Session], + local_map: Path, + tmp_path: Path, + monkeypatch: pytest.MonkeyPatch, +) -> None: + path, session = recorded + monkeypatch.setattr("flight_map.record.OpenSky.observations", Mock(side_effect=AssertionError)) + monkeypatch.setattr("flight_map.maps.request_json", Mock(side_effect=AssertionError)) + terrain = raster(tmp_path / "terrain.tif", 41, bounds=session.collection_bounds) + before = hashlib.sha256(path.read_bytes()).hexdigest() + output = tmp_path / "map.svg" + for name in ("map.svg", "second.svg"): + assert ( + render( + path, + session.id, + tmp_path / name, + local_map=local_map, + terrain=terrain, + terrain_reference="wgs84-ellipsoid", + confirm_geometric_ellipsoid=True, + ) + == 1 + ) + assert output.read_bytes() == (tmp_path / "second.svg").read_bytes() + assert hashlib.sha256(path.read_bytes()).hexdigest() == before + root = ET.parse(output).getroot() + assert root.attrib["width"] == "297mm" and root.attrib["height"] == "210mm" + assert root.attrib["fill"] == "none" + assert root.attrib["viewBox"] == "0 0 297 210" + forbidden = {"image", "text", "clipPath", "mask", "use"} + assert not any(element.tag.split("}")[-1] in forbidden for element in root.iter()) + assert {g.attrib.get("id") for g in root.findall(f"{{{SVG_NS}}}g")} == { + "map", + "flights", + "lettering", + } + flight_paths = root.findall(f".//{{{SVG_NS}}}g[@id='flights']//{{{SVG_NS}}}path") + assert flight_paths and all("C " in element.attrib["d"] for element in flight_paths) + flight = root.find(f".//{{{SVG_NS}}}g[@id='flight-1']") + assert flight is not None and flight.attrib["data-callsign"] == "KLM123" + assert flight.attrib["data-average-height-metres"] != "unavailable" + for element in root.iter(f"{{{SVG_NS}}}path"): + assert "clip-path" not in element.attrib + for x, y in re.findall(r"(-?\d+\.\d+),(-?\d+\.\d+)", element.attrib["d"]): + assert 0 <= float(x) <= 297 and 0 <= float(y) <= 210 + assert "HEIGHT" not in output.read_text() # Lettering is geometry, not SVG text. + with pytest.raises(ValueError, match="already exists"): + render(path, session.id, output, local_map=local_map) + + +def test_invalid_plot_time_and_empty_data(recorded: tuple[Path, Session], tmp_path: Path) -> None: + path, session = recorded + with pytest.raises(ValueError, match="within"): + render(path, session.id, tmp_path / "bad.svg", start=T0 - 10) + with pytest.raises(ValueError, match="No recorded"): + render(path, session.id, tmp_path / "bad.svg", start=T0 + 81, end=T0 + 85) + with pytest.raises(ValueError, match="positive"): + PlotOptions(width=float("nan")).validate() + with pytest.raises(ValueError, match="no room"): + PlotOptions(width=10).validate() + with pytest.raises(ValueError, match="must not overwrite"): + render(path, session.id, path, force=True) + with pytest.raises(ValueError, match="must not overwrite"): + render( + path, session.id, tmp_path / "terrain.tif", terrain=tmp_path / "terrain.tif", force=True + ) + + +def test_map_cache_and_errors(tmp_path: Path, monkeypatch: pytest.MonkeyPatch) -> None: + payload = { + "elements": [ + { + "type": "way", + "geometry": [ + {"lon": 4.8, "lat": 52.24}, + {"lon": 4.82, "lat": 52.24}, + ], + } + ] + } + fetch = Mock(return_value=payload) + monkeypatch.setattr("flight_map.maps.request_json", fetch) + assert map_data(DEFAULT_BOUNDS, tmp_path) == payload + assert map_data(DEFAULT_BOUNDS, tmp_path) == payload + fetch.assert_called_once() + with pytest.raises(ValueError, match="No drawable"): + map_lines({"elements": []}, Projection(DEFAULT_BOUNDS)) + with pytest.raises(ValueError, match="out geom"): + map_lines({"elements": [{"type": "way"}]}, Projection(DEFAULT_BOUNDS)) + fetch.return_value = {"elements": [], "remark": "runtime error"} + with pytest.raises(ValueError, match="could not complete"): + map_data(DEFAULT_BOUNDS, tmp_path / "other") + + +def test_cache_refresh_is_atomic_and_keys_include_area( + tmp_path: Path, + monkeypatch: pytest.MonkeyPatch, +) -> None: + old = {"elements": [], "version": 1} + new = {"elements": [], "version": 2} + fetch = Mock(return_value=old) + monkeypatch.setattr("flight_map.maps.request_json", fetch) + assert map_data(DEFAULT_BOUNDS, tmp_path) == old + entry = next(tmp_path.glob("osm-*.json")) + fetch.return_value = new + assert map_data(DEFAULT_BOUNDS, tmp_path) == old + assert fetch.call_count == 1 + assert map_data(DEFAULT_BOUNDS, tmp_path, refresh=True) == new + assert json.loads(entry.read_text()) == new + assert not list(tmp_path.glob("*.tmp")) + fetch.side_effect = RuntimeError("offline") + with pytest.raises(RuntimeError, match="offline"): + map_data(DEFAULT_BOUNDS, tmp_path, refresh=True) + assert json.loads(entry.read_text()) == new + assert map_data(DEFAULT_BOUNDS, tmp_path) == new + fetch.side_effect = None + map_data(Bounds(4.8, 52.2, 4.9, 52.3), tmp_path) + assert len(list(tmp_path.glob("osm-*.json"))) == 2 + entry.write_text("{") + with pytest.raises(ValueError, match="refresh-map"): + map_data(DEFAULT_BOUNDS, tmp_path) + assert map_data(DEFAULT_BOUNDS, tmp_path, refresh=True) == new + with pytest.raises(ValueError, match="cannot be combined"): + map_data(DEFAULT_BOUNDS, tmp_path, entry, refresh=True) + + +def test_interrupted_cache_write_preserves_existing_entry( + tmp_path: Path, + monkeypatch: pytest.MonkeyPatch, +) -> None: + payload = {"elements": []} + monkeypatch.setattr("flight_map.maps.request_json", Mock(return_value=payload)) + map_data(DEFAULT_BOUNDS, tmp_path) + entry = next(tmp_path.glob("osm-*.json")) + before = entry.read_bytes() + monkeypatch.setattr("flight_map.maps.os.replace", Mock(side_effect=OSError("write failed"))) + with pytest.raises(OSError, match="write failed"): + map_data(DEFAULT_BOUNDS, tmp_path, refresh=True) + assert entry.read_bytes() == before + assert not list(tmp_path.glob("*.tmp")) + + +@pytest.mark.parametrize( + "data", + [ + {"elements": [None]}, + { + "elements": [ + { + "type": "way", + "geometry": [ + {"lon": float("nan"), "lat": 52.2}, + {"lon": 4.8, "lat": 52.3}, + ], + } + ] + }, + {"elements": [{"type": "way", "geometry": [None, None]}]}, + ], +) +def test_invalid_map_geometry_is_reported(data: dict[str, Any]) -> None: + with pytest.raises(ValueError, match="Map"): + map_lines(data, Projection(DEFAULT_BOUNDS)) + + +def test_map_detail_is_selected_for_pen_width() -> None: + geometry = [ + {"lon": 4.80, "lat": 52.24}, + {"lon": 4.8001, "lat": 52.2401}, + {"lon": 4.80, "lat": 52.24}, + ] + data = { + "elements": [ + {"type": "way", "tags": {"building": "yes"}, "geometry": geometry}, + {"type": "way", "tags": {"highway": "footway"}, "geometry": geometry}, + {"type": "way", "tags": {"waterway": "stream"}, "geometry": geometry}, + ] + } + projection = Projection(DEFAULT_BOUNDS) + assert len(map_lines(data, projection, min_building_size=20)) == 1 + assert len(map_lines(data, projection, min_building_size=0)) == 2 + + +def test_opensky_parsing_and_auth(monkeypatch: pytest.MonkeyPatch) -> None: + for key in ("OPENSKY_CLIENT_ID", "OPENSKY_CLIENT_SECRET", "OPENSKY_ACCESS_TOKEN"): + monkeypatch.delenv(key, raising=False) + state = ["abc123", " KLM123 ", "NL", T0, T0, 4.82, 52.24, 900, False, 100, 30, 0, None, 1043] + fetch = Mock(return_value={"states": [state]}) + monkeypatch.setattr("flight_map.record.request_json", fetch) + client = OpenSky() + items = client.observations(DEFAULT_BOUNDS) + assert items[0].callsign == "KLM123" and items[0].geo_altitude == 1043 + assert "lamin=52.232025" in fetch.call_args.args[0] + state[5] = None + assert client.observations(DEFAULT_BOUNDS) == [] + fetch.return_value = {"states": None} + assert client.observations(DEFAULT_BOUNDS) == [] + fetch.return_value = {"wrong": []} + with pytest.raises(RuntimeError, match="missing"): + client.observations(DEFAULT_BOUNDS) + monkeypatch.setenv("OPENSKY_CLIENT_ID", "test-client") + with pytest.raises(ValueError, match="both"): + OpenSky() + monkeypatch.setenv("OPENSKY_CLIENT_SECRET", "test-secret") + fetch.side_effect = [{"access_token": "test-token", "expires_in": 300}, {"states": []}] + OpenSky().observations(DEFAULT_BOUNDS) + assert fetch.call_args.kwargs["headers"]["Authorization"] == "Bearer test-token" + + +def test_collection_checkpoints_without_map_services( + recorded: tuple[Path, Session], + monkeypatch: pytest.MonkeyPatch, +) -> None: + path, session = recorded + monkeypatch.setattr("flight_map.maps.request_json", Mock(side_effect=AssertionError)) + client = Mock(spec=OpenSky) + client.observations.return_value = [] + with database(path) as connection: + collect(connection, session, client, duration=0) + assert get_session(connection, session.id).status == "stopped" + client.observations.side_effect = RuntimeError("outage") + with pytest.raises(RuntimeError, match="outage"): + collect(connection, session, client, duration=0) + assert get_session(connection, session.id).status == "failed" + client.observations.side_effect = KeyboardInterrupt + collect(connection, session, client) + assert get_session(connection, session.id).status == "interrupted" + + +def test_cli_defaults_custom_bounds_and_resume( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch +) -> None: + path = tmp_path / "flights.sqlite" + collector = Mock() + monkeypatch.setattr("flight_map.record.collect", collector) + monkeypatch.setattr("flight_map.record.OpenSky", Mock()) + assert main(["collect", "--db", str(path), "--duration", "1"]) == 0 + session = collector.call_args.args[1] + assert session.bounds == DEFAULT_BOUNDS + custom = ["--bounds", "4.90", "52.30", "4.80", "52.20"] + assert main(["collect", "--db", str(path), *custom, "--duration", "1"]) == 0 + custom_session = collector.call_args.args[1] + assert custom_session.bounds == Bounds(4.8, 52.2, 4.9, 52.3) + assert main(["collect", "--db", str(path), "--resume", session.id, "--duration", "1"]) == 0 + assert collector.call_args.args[1].bounds == DEFAULT_BOUNDS + assert main(["collect", "--db", str(path), "--resume", session.id, *custom]) == 1 + assert main(["sessions", "--db", str(path)]) == 0 + + +def test_cli_validation_and_plot( + recorded: tuple[Path, Session], + local_map: Path, + tmp_path: Path, +) -> None: + path, session = recorded + assert ( + main( + [ + "plot", + "--db", + str(path), + "--session", + session.id, + "--map-data", + str(local_map), + "--output", + str(tmp_path / "cli.svg"), + ] + ) + == 0 + ) + assert ( + main( + [ + "plot", + "--db", + str(path), + "--session", + "missing", + "--output", + str(tmp_path / "bad.svg"), + ] + ) + == 1 + ) + with pytest.raises(SystemExit): + parser().parse_args(["collect", "--db", str(path), "--interval", "nan"]) + with pytest.raises(Exception, match="time zone"): + timestamp("2026-10-20T12:00:00") + assert timestamp("2026-10-20T12:00:00Z") == timestamp("2026-10-20T14:00:00+02:00") + + +def test_network_errors(monkeypatch: pytest.MonkeyPatch) -> None: + from email.message import Message + from urllib.error import HTTPError, URLError + + headers = Message() + headers["Retry-After"] = "60" + fetch = Mock(side_effect=HTTPError("https://example.invalid", 429, "limited", headers, None)) + monkeypatch.setattr("flight_map.network.urlopen", fetch) + with pytest.raises(RuntimeError, match="429.*60"): + request_json("https://example.invalid") + fetch.side_effect = URLError("offline") + with pytest.raises(RuntimeError, match="connection failed"): + request_json("https://example.invalid") diff --git a/tests/test_purge.py b/tests/test_purge.py new file mode 100644 index 0000000..92daa4e --- /dev/null +++ b/tests/test_purge.py @@ -0,0 +1,92 @@ +import sqlite3 +from pathlib import Path + +import pytest + +from flight_map.cli import main +from flight_map.geo import DEFAULT_BOUNDS +from flight_map.storage import ( + Observation, + create_session, + database, + get_session, + initialize, + load_observations, + purge_session, + save_observations, +) + + +def seed(path: Path) -> tuple[str, str]: + with database(path) as connection: + initialize(connection) + sessions = [create_session(connection, DEFAULT_BOUNDS, 2000, 30, 100) for _ in range(2)] + for session in sessions: + save_observations( + connection, + session, + [ + Observation("abc123", "TEST", 110, 4.82, 52.24, 900, 1000), + ], + 120, + ) + return sessions[0].id, sessions[1].id + + +def test_purge_only_selected_session(tmp_path: Path, capsys: pytest.CaptureFixture[str]) -> None: + path = tmp_path / "test.sqlite" + target, other = seed(path) + assert main(["purge", "--db", str(path), "--session", target, "--yes"]) == 0 + assert f"Purged session {target} and 1 observation(s)." in capsys.readouterr().out + with database(path, readonly=True) as connection: + with pytest.raises(ValueError, match="not found"): + get_session(connection, target) + assert load_observations(connection, target, 100, 120) == [] + assert get_session(connection, other).bounds == DEFAULT_BOUNDS + assert len(load_observations(connection, other, 100, 120)) == 1 + assert connection.execute("PRAGMA foreign_key_check").fetchall() == [] + + +def test_confirmation_and_unknown_session_preserve_data(tmp_path: Path) -> None: + path = tmp_path / "test.sqlite" + target, _ = seed(path) + before = path.read_bytes() + with pytest.raises(SystemExit): + main(["purge", "--db", str(path), "--session", target]) + assert path.read_bytes() == before + assert main(["purge", "--db", str(path), "--session", "missing", "--yes"]) == 1 + assert path.read_bytes() == before + + +def test_missing_database_is_not_created(tmp_path: Path) -> None: + path = tmp_path / "missing.sqlite" + assert main(["purge", "--db", str(path), "--session", "missing", "--yes"]) == 1 + assert not path.exists() + + +def test_purge_rolls_back_observations_if_session_delete_fails(tmp_path: Path) -> None: + path = tmp_path / "test.sqlite" + target, _ = seed(path) + with database(path) as connection: + connection.execute(""" + CREATE TRIGGER prevent_session_delete BEFORE DELETE ON sessions + BEGIN SELECT RAISE(ABORT, 'test failure'); END + """) + connection.commit() + with pytest.raises(sqlite3.IntegrityError, match="test failure"): + purge_session(connection, target) + assert get_session(connection, target).id == target + assert len(load_observations(connection, target, 100, 120)) == 1 + connection.execute("DROP TRIGGER prevent_session_delete") + connection.commit() + assert purge_session(connection, target) == 1 + + +def test_empty_session_can_be_purged(tmp_path: Path) -> None: + path = tmp_path / "test.sqlite" + with database(path) as connection: + initialize(connection) + session = create_session(connection, DEFAULT_BOUNDS, 2000, 30, 100) + assert purge_session(connection, session.id) == 0 + with pytest.raises(ValueError, match="not found"): + get_session(connection, session.id)