bandsaunter/bandsaunter/basemap.py
The Dust Council 87b0954f0c A sharper map, and aircraft that fade rather than blink out
Three things asked for, and a fourth found while doing them.

The map looked like a photograph of a map, and did so twice over.  The
window fetched at its own pixel size but for a box a third larger in each
direction -- the margin added to stop the ground blinking -- and then cut
the middle out, so every pixel was enlarged by two thirds.  It now asks
for enough pixels to cover the bigger box at the window's own detail.
Underneath that, both the window and the animation took the nearest source
pixel: the tile mosaic is commonly half again the size of the picture, so
most of every tile was thrown away and what survived was the aliasing.
Both now average the source pixels that fall in each output cell, done as
the difference of a running total rather than a loop.

An aircraft that goes quiet now fades instead of vanishing.  Taking it off
between one frame and the next says it stopped existing; fading says it
stopped talking, which is what happened.  It fades where it was last
actually seen and never along a reckoned track, because the reason for
giving up on it is that where it would be by now is a guess.  --fade sets
how long, and it is in the menu.  The window has alpha and fades smoothly;
an indexed picture cannot blend, so the animation gained a fourth ramp at
a seventh of full and fades in four steps, which at a second apart reads
as a fade.  A trail fades with the aircraft it belongs to, and the box
goes before the symbol does.

The aircraft's country of registration carries a flag now as well as the
two ends of its route, from the register where one answered and from the
address block otherwise.

And the fourth: the window's header counted an aircraft that had gone
quiet as overhead, which was saying more than had been heard.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016PsWPTweCT6pwxKngvVxcg
2026-09-04 19:04:09 -07:00

506 lines
21 KiB
Python

"""The ground under the aircraft: map tiles, fetched, cached and dimmed.
A flight path over a black rectangle says how the aircraft moved and nothing
about where it was. Over a coastline it says which airport it left. So the
map behind the animation is a real one: standard raster tiles, fetched once,
kept on disk, reprojected onto the picture and dimmed until the aircraft are
the brightest thing on it.
Three things follow from using somebody else's tile server, and all three are
obligations rather than options. Tiles are **cached** and never fetched
twice. Every request identifies the program in its User-Agent. And the
attribution the licence requires is drawn onto the picture, not left to a
readme nobody ships with a GIF. A drawing is capped at a few dozen tiles: an
aircraft map is a hobby drawing, not a reason to hammer a volunteer-funded
service.
The PNG decoding is here for the same reason the PNG writing is in
:mod:`bandsaunter.images`: a scanner that cannot draw a map because an
imaging library is missing is worse than one that draws it from zlib and
numpy, which is all this needs.
"""
from __future__ import annotations
import json
import math
import os
import struct
import time
import urllib.parse
import urllib.request
import zlib
from pathlib import Path
import numpy as np
from . import __version__
__all__ = ["decode_png", "tile_of", "choose_zoom", "fetch_tile", "mosaic",
"ground_under", "TILE_URL", "ATTRIBUTION", "MAX_TILES", "MAX_ZOOM",
"cache_dir", "PNGError", "airports_in", "AIRPORTS_URL"]
# The standard OpenStreetMap tiles. Any {z}/{x}/{y} server can be put here
# instead; nothing below knows anything about this one in particular.
TILE_URL = "https://tile.openstreetmap.org/{z}/{x}/{y}.png"
# Drawn onto every picture that used the tiles. The licence requires it and
# a GIF travels without its readme.
ATTRIBUTION = "MAP DATA (C) OPENSTREETMAP CONTRIBUTORS"
# A drawing is worth a few dozen tiles and no more. Past that the zoom is
# reduced instead: a coarser map still says where the coastline is.
MAX_TILES = 90
MAX_ZOOM = 13
MIN_ZOOM = 2
TILE_PIXELS = 256
USER_AGENT = (f"bandsaunter/{__version__} "
"(+https://github.com/topics/rtl-sdr; aircraft map drawing)")
# Politeness between requests to a volunteer-funded service. Only paid on a
# tile that was not already on the disk.
FETCH_PAUSE = 0.12
# Where to ask what aerodromes are in a piece of the world. The same OSM
# data the tiles are drawn from, asked as a question rather than a picture.
AIRPORTS_URL = "https://overpass-api.de/api/interpreter"
# One question covers a whole view and is kept for a month, because runways
# do not move. The same politeness as the tiles: cache it, say who is
# asking, and do not ask twice for the same thing.
AIRPORT_CACHE_DAYS = 30
# Bumped when the question changes, so that answers to the old one are asked
# again rather than believed.
AIRPORT_CACHE_VERSION = 2
# Enough to name what is under an aircraft and not so many that the map is a
# list of airstrips.
MOST_AIRPORTS = 40
class PNGError(ValueError):
"""A PNG this decoder cannot read."""
# ---------------------------------------------------------------------------
# Reading a PNG
# ---------------------------------------------------------------------------
_CHANNELS = {0: 1, 2: 3, 3: 1, 4: 2, 6: 4}
def decode_png(data: bytes) -> np.ndarray:
"""A PNG's pixels as an ``(h, w, 3)`` array of bytes.
Eight bits a channel and no interlacing, which is what every tile server
sends; anything else raises rather than being guessed at.
"""
if data[:8] != b"\x89PNG\r\n\x1a\x08"[:8] and data[:8] != \
b"\x89PNG\r\n\x1a\n":
raise PNGError("not a PNG")
width = height = depth = colour = interlace = 0
palette = None
body = bytearray()
at = 8
while at + 8 <= len(data):
length, tag = struct.unpack(">I4s", data[at:at + 8])
chunk = data[at + 8:at + 8 + length]
at += 12 + length
if tag == b"IHDR":
(width, height, depth, colour, _compression, _filter,
interlace) = struct.unpack(">IIBBBBB", chunk)
elif tag == b"PLTE":
palette = np.frombuffer(chunk, dtype=np.uint8).reshape(-1, 3)
elif tag == b"IDAT":
body += chunk
elif tag == b"IEND":
break
if depth != 8:
raise PNGError(f"{depth}-bit PNG: only 8 bits a channel is read here")
if interlace:
raise PNGError("interlaced PNG")
if colour not in _CHANNELS:
raise PNGError(f"colour type {colour}")
if not width or not height:
raise PNGError("no image")
channels = _CHANNELS[colour]
raw = _unfilter(zlib.decompress(bytes(body)), width, height, channels)
if colour == 3:
if palette is None:
raise PNGError("palette image with no palette")
return palette[np.clip(raw[:, :, 0], 0, len(palette) - 1)]
if colour == 0:
return np.repeat(raw, 3, axis=2)
if colour == 4:
return np.repeat(raw[:, :, :1], 3, axis=2)
return raw[:, :, :3]
def _unfilter(data: bytes, width: int, height: int,
channels: int) -> np.ndarray:
"""Undo the per-row filters PNG applies before compressing.
The five of them, from the specification. None and Up are whole-row
arithmetic; Sub is a running total along the row, which is a cumulative
sum once the bytes are grouped by which channel they belong to; Average
and Paeth each need the byte before them to have been worked out
already, so those two are the only ones that walk the row.
"""
stride = width * channels
if len(data) < height * (stride + 1):
raise PNGError("truncated image data")
out = np.zeros((height, stride), dtype=np.uint8)
previous = np.zeros(stride, dtype=np.uint8)
at = 0
for row in range(height):
kind = data[at]
line = np.frombuffer(data, dtype=np.uint8, count=stride,
offset=at + 1).astype(np.uint16)
at += stride + 1
if kind == 0:
current = line.astype(np.uint8)
elif kind == 1:
current = np.empty(stride, dtype=np.uint8)
for offset in range(channels):
current[offset::channels] = np.cumsum(
line[offset::channels], dtype=np.uint32) % 256
elif kind == 2:
current = ((line + previous) % 256).astype(np.uint8)
elif kind in (3, 4):
current = _walk_row(bytes(line.astype(np.uint8)), previous,
channels, kind)
else:
raise PNGError(f"filter type {kind}")
out[row] = current
previous = current
return out.reshape(height, width, channels)
def _walk_row(line: bytes, previous: np.ndarray, channels: int,
kind: int) -> np.ndarray:
"""Average and Paeth: each byte needs the one ``channels`` back."""
up = previous.tolist()
out = [0] * len(line)
for i, value in enumerate(line):
left = out[i - channels] if i >= channels else 0
above = up[i]
if kind == 3:
out[i] = (value + ((left + above) >> 1)) & 0xFF
else:
upleft = up[i - channels] if i >= channels else 0
base = left + above - upleft
da, db, dc = abs(base - left), abs(base - above), abs(base - upleft)
near = left if (da <= db and da <= dc) else (
above if db <= dc else upleft)
out[i] = (value + near) & 0xFF
return np.array(out, dtype=np.uint8)
# ---------------------------------------------------------------------------
# Which tiles
# ---------------------------------------------------------------------------
def tile_of(lat: float, lon: float, zoom: int) -> tuple[float, float]:
"""Where a coordinate falls in the tile grid, in fractional tiles.
Web Mercator, which is what every {z}/{x}/{y} tile server serves and is
not what this program's maps are drawn in -- hence the resampling
further down rather than a straight paste.
"""
lat = max(-85.05112878, min(85.05112878, lat))
n = float(2 ** zoom)
x = (lon + 180.0) / 360.0 * n
radians = math.radians(lat)
y = (1.0 - math.asinh(math.tan(radians)) / math.pi) / 2.0 * n
return x, y
def choose_zoom(south: float, west: float, north: float, east: float,
max_tiles: int = MAX_TILES, most: int = MAX_ZOOM) -> int:
"""The most detail that fits inside the tile budget."""
for zoom in range(min(most, MAX_ZOOM), MIN_ZOOM - 1, -1):
x0, y0 = tile_of(north, west, zoom)
x1, y1 = tile_of(south, east, zoom)
wide = int(math.floor(x1)) - int(math.floor(x0)) + 1
tall = int(math.floor(y1)) - int(math.floor(y0)) + 1
if wide * tall <= max_tiles:
return zoom
return MIN_ZOOM
# ---------------------------------------------------------------------------
# Fetching them, once
# ---------------------------------------------------------------------------
def cache_dir() -> Path:
root = os.environ.get("XDG_CACHE_HOME") or "~/.cache"
return Path(root).expanduser() / "bandsaunter" / "tiles"
def fetch_tile(zoom: int, x: int, y: int, url: str = TILE_URL,
timeout: float = 10.0, cache: Path | None = None) -> bytes | None:
"""One tile, from the disk if it has ever been fetched before.
Returns the PNG bytes, or None if it could not be had. A missing tile is
not an error: the map is drawn with a hole in it, which is better than no
map and much better than an exception in the middle of an animation.
"""
where = (cache if cache is not None else cache_dir()) / str(zoom) / str(x)
path = where / f"{y}.png"
try:
return path.read_bytes()
except OSError:
pass
request = urllib.request.Request(url.format(z=zoom, x=x, y=y),
headers={"User-Agent": USER_AGENT})
try:
with urllib.request.urlopen(request, timeout=timeout) as answer:
body = answer.read(2_000_000)
except Exception:
return None
try:
where.mkdir(parents=True, exist_ok=True)
tmp = path.with_suffix(".tmp")
tmp.write_bytes(body)
tmp.replace(path)
except OSError:
pass # an unwritable cache is not a reason to stop
return body
def mosaic(south: float, west: float, north: float, east: float, zoom: int,
fetch=fetch_tile, pause: float = FETCH_PAUSE, **kw):
"""Every tile the box touches, stitched into one image.
Returns the pixels and where their top-left corner sits in the world, in
tile-grid pixels at this zoom, so the resampling below can place them.
"""
x0, y0 = tile_of(north, west, zoom)
x1, y1 = tile_of(south, east, zoom)
left, top = int(math.floor(x0)), int(math.floor(y0))
right, bottom = int(math.floor(x1)), int(math.floor(y1))
span = 2 ** zoom
wide, tall = right - left + 1, bottom - top + 1
if wide <= 0 or tall <= 0 or wide * tall > MAX_TILES * 4:
return None, 0, 0
canvas = np.zeros((tall * TILE_PIXELS, wide * TILE_PIXELS, 3),
dtype=np.uint8)
got = 0
for row in range(tall):
for column in range(wide):
tx, ty = (left + column) % span, top + row
if not 0 <= ty < span:
continue
body = fetch(zoom, tx, ty, **kw)
if body is None:
continue
try:
tile = decode_png(body)
except (PNGError, zlib.error, ValueError):
continue
if tile.shape[0] != TILE_PIXELS or tile.shape[1] != TILE_PIXELS:
continue
canvas[row * TILE_PIXELS:(row + 1) * TILE_PIXELS,
column * TILE_PIXELS:(column + 1) * TILE_PIXELS] = tile
got += 1
if pause:
time.sleep(pause)
if not got:
return None, 0, 0
return canvas, left * TILE_PIXELS, top * TILE_PIXELS
# ---------------------------------------------------------------------------
# What is on the ground
# ---------------------------------------------------------------------------
def _resample(values: np.ndarray, edges: np.ndarray, axis: int) -> np.ndarray:
"""Average each output cell over the source pixels that fall in it.
A box filter, done as the difference of a running total so the whole
axis is two passes rather than a loop. Where the source is coarser than
the output -- a map zoomed in further than the tiles go -- a cell covers
less than one source pixel, and it takes that one.
"""
length = values.shape[axis]
starts = np.clip(np.floor(edges[:-1]).astype(np.int64), 0, length - 1)
ends = np.clip(np.ceil(edges[1:]).astype(np.int64), 1, length)
ends = np.maximum(ends, starts + 1)
running = np.cumsum(values, axis=axis, dtype=np.float64)
pad = np.zeros_like(np.take(running, [0], axis=axis))
running = np.concatenate([pad, running], axis=axis)
total = (np.take(running, ends, axis=axis)
- np.take(running, starts, axis=axis))
counts = (ends - starts).astype(np.float64)
shape = [1] * values.ndim
shape[axis] = counts.size
return (total / counts.reshape(shape)).astype(np.float32)
def _airport_cache(box) -> Path:
root = os.environ.get("XDG_CACHE_HOME") or "~/.cache"
name = "_".join(f"{round(v, 1):+06.1f}" for v in box)
return Path(root).expanduser() / "bandsaunter" / "airports" / f"{name}.json"
def airports_in(south: float, west: float, north: float, east: float,
url: str = AIRPORTS_URL, timeout: float = 45.0,
cache: Path | None = None, ask=None) -> list[dict]:
"""Every aerodrome in a piece of the world, with a code and a position.
Asked of the same OpenStreetMap data the tiles are drawn from, as a
question rather than a picture, and kept on disk afterwards: a runway
does not move, so one question covers a view for a month.
An aerodrome with no code is left out. A map wants to say *which*
airport an aircraft is over, and there are a great many landing strips
with a name and nothing else; the ones worth marking have a code.
"""
box = (round(south, 1), round(west, 1), round(north, 1), round(east, 1))
path = cache if cache is not None else _airport_cache(box)
try:
body = json.loads(path.read_text())
if int(body.get("version") or 0) >= AIRPORT_CACHE_VERSION and \
time.time() - float(body.get("fetched_at") or 0) < \
AIRPORT_CACHE_DAYS * 86_400:
return body.get("airports") or []
except (OSError, ValueError):
pass
found: list[dict] = []
try:
raw = (ask or _ask_overpass)(box, url, timeout)
for element in (raw or {}).get("elements", []):
tags = element.get("tags") or {}
code = _airport_code(tags)
lat = element.get("lat") or (element.get("center") or {}).get("lat")
lon = element.get("lon") or (element.get("center") or {}).get("lon")
if not code or lat is None or lon is None:
continue
found.append({"code": code,
"name": (tags.get("name") or "").strip(),
"icao": bool(tags.get("icao")),
"latitude": float(lat), "longitude": float(lon)})
except Exception:
return [] # a map with no airports on it, not a crash
# The ones with a real ICAO code first, since those are the ones an
# aircraft is likely to be flying to or from.
found.sort(key=lambda a: (0 if a.get("icao") else 1, a["code"]))
# One airport is often tagged twice -- a point for the terminal and an
# outline for the field -- and marking it twice writes its name over
# itself.
once: dict[str, dict] = {}
for one in found:
once.setdefault(one["code"], one)
found = list(once.values())[:MOST_AIRPORTS]
try:
path.parent.mkdir(parents=True, exist_ok=True)
path.write_text(json.dumps({"fetched_at": time.time(),
"version": AIRPORT_CACHE_VERSION,
"airports": found}))
except OSError:
pass
return found
def _airport_code(tags: dict) -> str:
"""What to call an aerodrome, or nothing if it has no code at all.
The ICAO code where it has one. Otherwise a reference, and only when it
is four letters: a great many landing strips carry a local identifier
like "14AZ" or "MX-0492", which names nothing anybody would recognise
and turns a map into a list of airstrips.
"""
icao = (tags.get("icao") or "").strip().upper()
if len(icao) == 4 and icao.isalpha():
return icao
ref = (tags.get("ref") or "").strip().upper()
if len(ref) == 4 and ref.isalpha():
return ref
return ""
def _ask_overpass(box, url: str, timeout: float) -> dict:
"""The one question this program asks Overpass."""
south, west, north, east = box
where = f"({south},{west},{north},{east})"
# Relations as well as nodes and ways: a big airport is a relation more
# often than not -- Tucson International and Davis-Monthan both are --
# so asking only for the other two finds every airstrip in the county
# and misses the two the county is known for.
query = ("[out:json][timeout:25];"
f'(node["aeroway"="aerodrome"]{where};'
f' way["aeroway"="aerodrome"]{where};'
f' relation["aeroway"="aerodrome"]{where};);'
"out center tags;")
request = urllib.request.Request(
url, data=urllib.parse.urlencode({"data": query}).encode(),
headers={"User-Agent": USER_AGENT})
with urllib.request.urlopen(request, timeout=timeout) as answer:
return json.loads(answer.read(4_000_000).decode("utf8", "replace"))
# ---------------------------------------------------------------------------
# Putting it under the picture
# ---------------------------------------------------------------------------
def ground_under(south: float, west: float, north: float, east: float,
width: int, height: int, shades: int = 32,
fetch=fetch_tile, zoom: int | None = None,
**kw) -> np.ndarray | None:
"""The map for one picture, as ``shades`` levels of brightness.
The tiles are Web Mercator and the picture is not, so every output pixel
asks the mosaic where its own latitude and longitude landed rather than
the mosaic being pasted in. Over the couple of hundred miles a receiver
hears, the difference is a few pixels of drift at the top of the frame --
which is a few pixels an aircraft would be drawn wrong by, and the whole
point of putting a coastline under it is that the coastline is where the
aircraft was.
Returns None when nothing could be fetched, which the caller draws as the
plain grid it drew before.
"""
if width < 1 or height < 1 or north <= south or east <= west:
return None
zoom = choose_zoom(south, west, north, east) if zoom is None else zoom
tiles, origin_x, origin_y = mosaic(south, west, north, east, zoom,
fetch=fetch, **kw)
if tiles is None:
return None
# The edges of each output pixel rather than its middle, so that what
# lands in it can be averaged. Taking the nearest source pixel instead
# throws away most of a tile -- the mosaic is commonly half again the
# size of the picture -- and what survives is the aliasing: hard, broken
# lettering and roads that come and go along their length.
lons = west + (east - west) * np.arange(width + 1) / width
lats = north - (north - south) * np.arange(height + 1) / height
span = float(2 ** zoom) * TILE_PIXELS
xs = (lons + 180.0) / 360.0 * span - origin_x
clipped = np.clip(lats, -85.05112878, 85.05112878)
ys = (1.0 - np.arcsinh(np.tan(np.radians(clipped))) / math.pi) / 2.0 \
* span - origin_y
# Brightness only, inverted, and dimmed. Inverted because a printed map
# is ink on white paper and this picture is the other way round: the
# things drawn on the map -- coastlines, roads, the names of towns --
# are the dark parts of a tile, and they are what should show against a
# night background. Dimmed because the map is the ground under the
# aircraft rather than the subject: anything drawn on top has to stay the
# brightest thing on the picture.
whole = (0.299 * tiles[:, :, 0] + 0.587 * tiles[:, :, 1]
+ 0.114 * tiles[:, :, 2]).astype(np.float32)
luma = _resample(_resample(whole, ys, axis=0), xs, axis=1)
low, high = float(luma.min()), float(luma.max())
if high - low < 1.0:
levels = np.zeros_like(luma)
else:
levels = 1.0 - (luma - low) / (high - low)
return np.clip((levels * (shades - 1)).round(), 0,
shades - 1).astype(np.uint8)