This repository has no description
57 kB
1411 lines
1"""Render a live world map wallpaper: oblique Mercator with a lens on home.
2
3The globe is rotated so home (default Berlin) is in the middle of the screen
4on the projection's equator, then shown in Mercator filling the screen edge
5to edge. Around home a smooth radial lens magnifies the map and blends back
6into the plain map at its rim.
7
8Imagery is fetched over WMS at the resolution each pixel needs: a few local
9levels around home, one global level, and downsampled global mipmaps. The
10night side is the day map darkened, plus city lights as a glow. Every fetch
11is cached with a TTL; a failing source falls back to its last good copy or
12is skipped, so one dead upstream no longer kills the whole wallpaper.
13"""
14
15from __future__ import annotations
16
17import argparse
18import dataclasses
19import concurrent.futures
20import datetime as dt
21import hashlib
22import io
23import json
24import math
25import os
26import sys
27import time
28import urllib.parse
29import urllib.request
30import warnings
31from dataclasses import dataclass, field
32from pathlib import Path
33from typing import Callable
34
35import numpy as np
36from PIL import Image, ImageDraw, ImageFilter, ImageFont
37
38R_EARTH = 6371.0088 # km
39KM_PER_DEG = math.pi * R_EARTH / 180
40GLOBAL_WIDTH = 8192 # px of the global equirectangular level
41# px per WMS request side. The servers render on the fly at a roughly fixed
42# px/s per request (EOX: ~0.35 Mpx/s) but serve requests in parallel, so
43# many medium tiles fetched concurrently are much faster than few big ones
44MAX_TILE = 2048
45TILE_WORKERS = 4 # concurrent requests per level
46MAX_LEVEL_SIDE = 8192 # px per local level side
47USER_AGENT = "realwallpaper (+https://github.com/Lassulus/superconfig)"
48
49GIBS = "https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi"
50EUMETSAT = "https://view.eumetsat.int/geoserver/wms"
51EOX = "https://tiles.maps.eox.at/wms"
52
53
54def log(msg: str) -> None:
55 print(msg, file=sys.stderr, flush=True)
56
57
58def smoothstep(e0: float, e1: float, x: np.ndarray) -> np.ndarray:
59 t = np.clip((x - e0) / (e1 - e0), 0.0, 1.0)
60 return t * t * (3.0 - 2.0 * t)
61
62
63# --------------------------------------------------------------------------
64# projection
65
66
67def solve_cylinder(half_height: float) -> float:
68 """Parameter a of the cylindrical y(lat) = ln tan(pi/4 + a lat/2) / a
69 whose poles sit exactly half_height (radians of longitude) from the
70 equator. a -> 0 is plate carrée (pi/2), a = 1 is Mercator (infinite);
71 Miller is a = 0.8. Wider screens than 2:1 fall back to plate carrée."""
72 lo, hi = 1e-6, 1 - 1e-9
73 if half_height <= math.pi / 2:
74 return lo
75 for _ in range(200):
76 a = (lo + hi) / 2
77 if math.log(math.tan(math.pi / 4 + a * math.pi / 4)) / a > half_height:
78 hi = a
79 else:
80 lo = a
81 return (lo + hi) / 2
82
83
84def rotate_to_home(lat0: float, lat, dlon):
85 """Geographic (lat, lon - lon0) -> frame where home is lat' = lon' = 0.
86
87 Rotating about the east-west axis keeps local north pointing up at home.
88 All angles in radians.
89 """
90 x = np.cos(lat) * np.cos(dlon)
91 y = np.cos(lat) * np.sin(dlon)
92 z = np.sin(lat)
93 s0, c0 = math.sin(lat0), math.cos(lat0)
94 xr = x * c0 + z * s0
95 zr = -x * s0 + z * c0
96 return np.arcsin(np.clip(zr, -1.0, 1.0)), np.arctan2(y, xr)
97
98
99def rotate_from_home(lat0: float, lat_r, lon_r):
100 """Inverse of rotate_to_home: returns (lat, lon - lon0) in radians."""
101 xr = np.cos(lat_r) * np.cos(lon_r)
102 y = np.cos(lat_r) * np.sin(lon_r)
103 zr = np.sin(lat_r)
104 s0, c0 = math.sin(lat0), math.cos(lat0)
105 x = xr * c0 - zr * s0
106 z = xr * s0 + zr * c0
107 return np.arcsin(np.clip(z, -1.0, 1.0)), np.arctan2(y, x)
108
109
110class LensMap:
111 """Screen-filling oblique cylindrical map around home, with a lens.
112
113 The globe is rotated so home sits on the projection's equator in the
114 middle of the screen (local north up); its poles land top and bottom
115 only when home is on the real equator. The cylinder is Mercator-like
116 (conformal along that equator) but with finite poles, y = ln tan(pi/4 +
117 a lat/2) / a, with a picked so the whole globe fits the screen exactly:
118 nothing is cut off, only squeezed towards the top and bottom edges.
119
120 Projections (--projection):
121 - lens: the above plus the magnifying lens around home (default)
122 - plain: the above without the lens, no distortion beyond the cylinder
123 - mercator: the classic map, north up, true Mercator centred on the
124 equator and the home longitude; the poles are cut off (about +-71
125 degrees fit 16:9) and there is no lens
126
127 The lens is focus + glue + context (Yamamoto et al. 2009; the same idea
128 as Snyder's 1987 "magnifying-glass" azimuthal projections): a lens
129 "radius" n is measured as a circle, or following the screen's aspect
130 as an ellipse (p = 2) or rounded rectangle (p = 4), the norm
131
132 n = ((|dx| / ax)^p + (|dy| / ay)^p)^(1/p)
133
134 (ax, ay = half the screen size times --lens). A screen point at n shows
135 the map point at h(n) = n / z(n) along the same ray, where z is the
136 local zoom: M in the focus (n < f, undistorted), 1 from n = 1 on (plain
137 map), and in between, with t = (n - f) / (1 - f),
138
139 ln z(t) = ln M * (1 - I(t) / I(1)), I(t) = int_0^t u^(a-1) (1 - u) du
140
141 How squashed shapes are is the ratio of radial to tangential scale,
142 1 - t z'/z; its excess over 1 is proportional to t^a (1 - t): none at home, largest at
143 t = a / (a + 1) =: --falloff, back to 1 at the lens edge. So --falloff
144 says where the distortion sits: low keeps it near home, high pushes it
145 out towards the border and leaves the middle closer to a plain zoom.
146 z never increases outwards, so h' > 0 and the lens never folds.
147 """
148
149 def __init__(
150 self,
151 lat: float,
152 lon: float,
153 width: int,
154 height: int,
155 magnification: float,
156 lens: float,
157 shape: str,
158 focus: float,
159 falloff: float,
160 projection: str = "lens",
161 ) -> None:
162 if magnification < 1:
163 raise SystemExit("--magnification must be >= 1")
164 if not 0 <= focus < 1:
165 raise SystemExit("--focus must be in [0, 1)")
166 if not 0 < falloff < 1:
167 raise SystemExit("--falloff must be in (0, 1)")
168 self.lat0 = math.radians(lat)
169 self.lon0 = math.radians(lon)
170 self.width = width
171 self.height = height
172 self.m = magnification if projection == "lens" else 1.0
173 self.k = 2 * math.pi / width # radians of longitude per px
174 if projection == "mercator":
175 self.rot_lat = 0.0 # no rotation: the real equator is the axis
176 self.a = 1.0 # true Mercator
177 else:
178 self.rot_lat = self.lat0 # home on the projection's equator
179 self.a = solve_cylinder(math.pi * height / width)
180 self.px = width / 2
181 self.py = height / 2
182 self.shape = shape
183 if shape == "circle":
184 self.ax = self.ay = lens * height / 2
185 else:
186 self.ax = lens * width / 2
187 self.ay = lens * height / 2
188 # tabulate the profile; forward() inverts it by interpolation
189 self._s = np.linspace(0.0, 1.0, 8001)
190 t = np.clip((self._s - focus) / (1 - focus), 0.0, 1.0)
191 a = falloff / (1 - falloff)
192 with np.errstate(divide="ignore", invalid="ignore"):
193 f_int = np.where(t > 0, t ** (a - 1) * (1 - t), 0.0)
194 integral = np.concatenate([[0.0], np.cumsum((f_int[1:] + f_int[:-1]) / 2 * np.diff(t))])
195 zoom = np.exp(math.log(self.m) * (1 - integral / max(integral[-1], 1e-12)))
196 self._h = self._s / zoom
197 self._dh = np.gradient(self._h, self._s)
198 # sample the screen to plan which resolutions are needed where
199 ys, xs = np.mgrid[0:height:6, 0:width:6] + 3.0
200 lat_g, lon_g, q = self.inverse(xs.ravel(), ys.ravel())
201 s = self.norm(xs.ravel() - self.px, ys.ravel() - self.py)
202 self._plan_q = q[s < 1]
203 self._plan_d = great_circle_km(lat, lon, lat_g[s < 1], lon_g[s < 1])
204 self.q0 = float(q.min())
205 self.q_max = float(q.max())
206
207 def norm(self, dx, dy):
208 if self.shape == "circle":
209 return np.hypot(dx, dy) / self.ax
210 if self.shape == "ellipse":
211 return np.hypot(dx / self.ax, dy / self.ay)
212 return ((np.abs(dx) / self.ax) ** 4 + (np.abs(dy) / self.ay) ** 4) ** 0.25
213
214 def h(self, s):
215 return np.interp(s, self._s, self._h)
216
217 def dh(self, s):
218 return np.interp(s, self._s, self._dh)
219
220 def distance_for_demand(self, q: float) -> float:
221 """Radius (km) around home inside which some pixel needs finer than q."""
222 need = self._plan_d[self._plan_q < q]
223 return float(min(need.max(), 3500.0)) if len(need) else 0.0
224
225 def inverse(self, xs: np.ndarray, ys: np.ndarray):
226 """Screen px -> lat, lon (degrees) and ground size of the px (km)."""
227 dx, dy = xs - self.px, ys - self.py
228 s = self.norm(dx, dy)
229 inside = s < 1
230 si = np.where(inside, s, 1.0)
231 with np.errstate(divide="ignore", invalid="ignore"):
232 tangential = np.where(si > 1e-9, self.h(si) / si, 1 / self.m)
233 radial = self.dh(si)
234 scale = np.where(inside, tangential, 1.0)
235 wx = self.px + dx * scale
236 wy = self.py + dy * scale
237 lon_r = (wx - self.px) * self.k
238 u = (self.py - wy) * self.k
239 lat_r = (2 / self.a) * (np.arctan(np.exp(self.a * u)) - math.pi / 4)
240 lat_r = np.clip(lat_r, -math.pi / 2, math.pi / 2)
241 lat, dlon = rotate_from_home(self.rot_lat, lat_r, lon_r)
242 lon = (self.lon0 + dlon + math.pi) % (2 * math.pi) - math.pi
243 # the lens squeezes radially more than tangentially in its ring; use
244 # the finer scale for sharpness, but not below a quarter of the other
245 fine = np.minimum(tangential, radial)
246 coarse = np.maximum(tangential, radial)
247 lens_f = np.where(inside, np.maximum(fine, coarse / 4), 1.0)
248 # cylinder: east-west px size R k cos(lat), north-south R k cos(a lat)
249 base_f = np.maximum(np.cos(lat_r), np.cos(self.a * lat_r) / 4)
250 q = R_EARTH * self.k * base_f * lens_f
251 return np.degrees(lat), np.degrees(lon), np.maximum(q, 1e-3)
252
253 def forward(self, lat_deg, lon_deg):
254 lat = np.radians(np.asarray(lat_deg, dtype=np.float64))
255 dlon = np.radians(np.asarray(lon_deg, dtype=np.float64)) - self.lon0
256 lat_r, lon_r = rotate_to_home(self.rot_lat, lat, dlon)
257 lat_r = np.clip(lat_r, -math.radians(89.5), math.radians(89.5))
258 wx = self.px + lon_r / self.k
259 wy = self.py - np.log(np.tan(math.pi / 4 + self.a * lat_r / 2)) / (self.a * self.k)
260 dx, dy = wx - self.px, wy - self.py
261 r = self.norm(dx, dy)
262 s = np.interp(np.minimum(r, 1.0), self._h, self._s)
263 with np.errstate(divide="ignore", invalid="ignore"):
264 scale = np.where(r < 1, np.where(r > 1e-12, s / r, self.m), 1.0)
265 return self.px + dx * scale, self.py + dy * scale
266
267
268def great_circle_km(lat0: float, lon0: float, lat, lon) -> np.ndarray:
269 la0, lo0 = math.radians(lat0), math.radians(lon0)
270 la, lo = np.radians(lat), np.radians(lon)
271 c = math.sin(la0) * np.sin(la) + math.cos(la0) * np.cos(la) * np.cos(lo - lo0)
272 return np.arccos(np.clip(c, -1.0, 1.0)) * R_EARTH
273
274
275# --------------------------------------------------------------------------
276# sources and levels
277
278
279@dataclass(frozen=True)
280class Wms:
281 key: str
282 base: str
283 layer: str
284 native_km: float
285 ttl: float
286 fmt: str = "image/png"
287 time: str | None = None
288 # TIME=<today - n days> for each n, newest first; gaps are filled from
289 # older days (daily polar-orbiter products are full of swath gaps)
290 days_back: tuple[int, ...] = ()
291 # fetch this many times finer and box-filter down: some servers (EOX)
292 # pick single pixels when asked for a coarse scale, which aliases
293 oversample: int = 1
294 # per-channel gain to match the colours of the other day source
295 gain: tuple[float, float, float] | None = None
296
297 def times(self, when: dt.datetime) -> list[str | None]:
298 if self.days_back:
299 return [(when - dt.timedelta(days=n)).strftime("%Y-%m-%d") for n in self.days_back]
300 return [self.time]
301
302 def url(self, bbox: tuple[float, float, float, float], w: int, h: int, t: str | None) -> str:
303 s, west, n, e = bbox
304 params = {
305 "SERVICE": "WMS",
306 "VERSION": "1.3.0",
307 "REQUEST": "GetMap",
308 "LAYERS": self.layer,
309 "STYLES": "",
310 "CRS": "EPSG:4326",
311 "BBOX": f"{s:.6f},{west:.6f},{n:.6f},{e:.6f}",
312 "WIDTH": str(w),
313 "HEIGHT": str(h),
314 "FORMAT": self.fmt,
315 "TRANSPARENT": "TRUE",
316 }
317 if t:
318 params["TIME"] = t
319 return self.base + "?" + urllib.parse.urlencode(params)
320
321
322@dataclass
323class Level:
324 res: float # km per px
325 south: float
326 west: float
327 north: float
328 east: float
329 width: int
330 height: int
331 wrap: bool
332 data: np.ndarray | None = None
333
334 @property
335 def bbox(self) -> tuple[float, float, float, float]:
336 return (self.south, self.west, self.north, self.east)
337
338
339def local_level(proj: LensMap, res: float, radius_km: float) -> Level:
340 lat = math.degrees(proj.lat0)
341 lon = math.degrees(proj.lon0)
342 dlat = math.degrees(radius_km / R_EARTH)
343 south, north = max(-90.0, lat - dlat), min(90.0, lat + dlat)
344 coslat = math.cos(proj.lat0)
345 s = math.sin(radius_km / R_EARTH) / coslat if coslat > 1e-6 else 2.0
346 if north >= 90.0 or south <= -90.0 or s >= 1.0:
347 west, east = -180.0, 180.0
348 else:
349 dlon = math.degrees(math.asin(s))
350 west, east = lon - dlon, lon + dlon
351 mid = math.cos(math.radians((south + north) / 2))
352 width_km = (east - west) * KM_PER_DEG * max(mid, 0.05)
353 height_km = (north - south) * KM_PER_DEG
354 res = max(res, width_km / MAX_LEVEL_SIDE, height_km / MAX_LEVEL_SIDE)
355 return Level(
356 res=res,
357 south=south,
358 west=west,
359 north=north,
360 east=east,
361 width=max(2, math.ceil(width_km / res)),
362 height=max(2, math.ceil(height_km / res)),
363 wrap=False,
364 )
365
366
367def global_level(width: int) -> Level:
368 return Level(
369 res=360 / width * KM_PER_DEG,
370 south=-90.0,
371 west=-180.0,
372 north=90.0,
373 east=180.0,
374 width=width,
375 height=width // 2,
376 wrap=True,
377 )
378
379
380def plan_local(proj: LensMap, native_km: float) -> list[Level]:
381 """Local levels, finest first, each 4x coarser than the previous one."""
382 g = global_level(GLOBAL_WIDTH).res
383 resolutions = [max(proj.q0, native_km)]
384 while resolutions[-1] * 4 < g / 1.5:
385 resolutions.append(resolutions[-1] * 4)
386 if resolutions[0] >= g / 1.5:
387 return []
388 levels = []
389 for i, res in enumerate(resolutions):
390 nxt = resolutions[i + 1] if i + 1 < len(resolutions) else g
391 radius = proj.distance_for_demand(nxt) * 1.1
392 if radius > 0:
393 levels.append(local_level(proj, res, radius))
394 return levels
395
396
397def mipmaps(level: Level, q_max: float) -> list[Level]:
398 out = []
399 cur = level
400 while cur.res < q_max and cur.height >= 64 and cur.data is not None:
401 d = cur.data.astype(np.uint16)
402 d = (d[0::2, 0::2] + d[1::2, 0::2] + d[0::2, 1::2] + d[1::2, 1::2] + 2) // 4
403 cur = Level(
404 res=cur.res * 2,
405 south=-90.0,
406 west=-180.0,
407 north=90.0,
408 east=180.0,
409 width=cur.width // 2,
410 height=cur.height // 2,
411 wrap=True,
412 data=d.astype(np.uint8),
413 )
414 out.append(cur)
415 return out
416
417
418# --------------------------------------------------------------------------
419# raster preparation: RGBA -> channels, last channel is data coverage, value
420# channels are premultiplied by coverage. Bump PREP_VERSION whenever a prep
421# function changes so cached prepared rasters are rebuilt from the raw ones.
422
423PREP_VERSION = 3
424
425
426class Colormap:
427 """GIBS colormap: maps the rendered colours back to a 0..1 data value."""
428
429 def __init__(self, path: Path) -> None:
430 import xml.etree.ElementTree as ET
431
432 entries = [
433 e
434 for e in ET.parse(path).iter()
435 if e.tag.endswith("ColorMapEntry")
436 and e.get("transparent") != "true"
437 and e.get("nodata") != "true"
438 ]
439 if not entries:
440 raise ValueError(f"{path}: no colormap entries")
441 self.rgb = np.array([[int(c) for c in e.get("rgb", "").split(",")] for e in entries])
442 self.value = np.linspace(0.0, 1.0, len(entries))
443
444 def decode(self, rgba: np.ndarray) -> np.ndarray:
445 out = np.zeros(rgba.shape[:2], np.float32)
446 mask = rgba[..., 3] > 0
447 px = rgba[mask][:, :3].astype(np.int64)
448 if len(px) == 0:
449 return out
450 uniq, inv = np.unique((px[:, 0] << 16) | (px[:, 1] << 8) | px[:, 2], return_inverse=True)
451 colors = np.stack([uniq >> 16, (uniq >> 8) & 255, uniq & 255], axis=-1)
452 values = np.empty(len(uniq))
453 for i in range(0, len(uniq), 4096):
454 d = ((colors[i : i + 4096, None, :] - self.rgb[None]) ** 2).sum(-1)
455 values[i : i + 4096] = self.value[d.argmin(1)]
456 out[mask] = values[inv]
457 return out
458
459
460def to_u8(x: np.ndarray) -> np.ndarray:
461 return np.round(np.clip(x, 0, 1) * 255).astype(np.uint8)
462
463
464def prep_rgb(layer: "Layer", rgba: np.ndarray, cmap: Colormap | None) -> np.ndarray:
465 a = rgba[..., 3:4].astype(np.uint16)
466 rgb = (rgba[..., :3].astype(np.uint16) * a + 127) // 255
467 return np.concatenate([rgb.astype(np.uint8), rgba[..., 3:4]], axis=-1)
468
469
470def prep_lights(layer: "Layer", rgba: np.ndarray, cmap: Colormap | None) -> np.ndarray:
471 """Black Marble -> city light intensity; drops its dim blue land texture."""
472 red = rgba[..., 0].astype(np.float32) / 255
473 full = np.full(rgba.shape[:2], 255, np.uint8)
474 return np.stack([to_u8(np.clip((red - 0.12) / 0.75, 0, 1) ** 1.2), full], axis=-1)
475
476
477def prep_cloud(layer: "Layer", rgba: np.ndarray, cmap: Colormap | None) -> np.ndarray:
478 """Infrared 10.8 µm brightness (cold = bright) -> cloud opacity."""
479 v = rgba[..., 0].astype(np.float32) / 255
480 cov = rgba[..., 3].astype(np.float32) / 255
481 assert layer.transform is not None
482 return np.stack([to_u8(layer.transform(v) * cov), rgba[..., 3]], axis=-1)
483
484
485def prep_scalar(layer: "Layer", rgba: np.ndarray, cmap: Colormap | None) -> np.ndarray:
486 """Colormapped data value (or plain presence) -> opacity, full coverage."""
487 v = cmap.decode(rgba) if cmap else rgba[..., 3].astype(np.float32) / 255
488 if layer.transform is not None:
489 v = layer.transform(v) * (rgba[..., 3] > 0)
490 full = np.full(rgba.shape[:2], 255, np.uint8)
491 return np.stack([to_u8(v), full], axis=-1)
492
493
494def prep_tint(layer: "Layer", rgba: np.ndarray, cmap: Colormap | None) -> np.ndarray:
495 """Colormapped value plus where it is defined, full coverage."""
496 assert cmap is not None
497 present = (rgba[..., 3] > 0).astype(np.float32)
498 full = np.full(rgba.shape[:2], 255, np.uint8)
499 return np.stack([to_u8(cmap.decode(rgba) * present), to_u8(present), full], axis=-1)
500
501
502def combine_max(layers: list[np.ndarray]) -> np.ndarray:
503 out = layers[0].copy()
504 for older in layers[1:]:
505 np.maximum(out, older, out=out)
506 return out
507
508
509def combine_fill(layers: list[np.ndarray]) -> np.ndarray:
510 """Fill where the newest has no data (channel -2 = presence) from older."""
511 out = layers[0].copy()
512 for older in layers[1:]:
513 gap = out[..., -2] == 0
514 out[gap] = older[gap]
515 return out
516
517
518# --------------------------------------------------------------------------
519# fetching
520
521
522class FetchError(Exception):
523 pass
524
525
526def http_get(url: str, timeout: float = 120) -> tuple[bytes, str]:
527 last: Exception | None = None
528 for attempt in range(3):
529 try:
530 req = urllib.request.Request(url, headers={"User-Agent": USER_AGENT})
531 with urllib.request.urlopen(req, timeout=timeout) as resp:
532 return resp.read(), resp.headers.get("Content-Type", "")
533 except Exception as e: # noqa: BLE001 - network errors come in many types
534 last = e
535 time.sleep(2 * (attempt + 1))
536 raise FetchError(f"{url}: {last}")
537
538
539def fetch_image(url: str) -> np.ndarray:
540 body, ctype = http_get(url)
541 if not ctype.startswith("image/"):
542 raise FetchError(f"{url}: got {ctype}: {body[:200]!r}")
543 return np.asarray(Image.open(io.BytesIO(body)).convert("RGBA"))
544
545
546def split_columns(level: Level) -> list[tuple[int, int]]:
547 px = (level.east - level.west) / level.width
548 cuts = {0, level.width}
549 for meridian in (-180.0, 180.0):
550 if level.west < meridian < level.east and not level.wrap:
551 cuts.add(round((meridian - level.west) / px))
552 c = 0
553 while c < level.width:
554 cuts.add(c)
555 c += MAX_TILE
556 edges = sorted(x for x in cuts if 0 <= x <= level.width)
557 return [(a, b) for a, b in zip(edges, edges[1:]) if b > a]
558
559
560def fetch_level(src: Wms, level: Level, t: str | None) -> np.ndarray:
561 factor = max(1, min(src.oversample, 2 * MAX_LEVEL_SIDE // max(level.width, level.height)))
562 big = dataclasses.replace(level, width=level.width * factor, height=level.height * factor,
563 data=None)
564 img = fetch_tiles(src, big, t)
565 if factor > 1:
566 img = np.asarray(Image.fromarray(img, "RGBA").reduce(factor))
567 return img
568
569
570def apply_gain(rgba: np.ndarray, gain: tuple[float, float, float] | None) -> np.ndarray:
571 if gain is None:
572 return rgba
573 out = rgba.copy()
574 out[..., :3] = to_u8(rgba[..., :3] / 255.0 * np.array(gain))
575 return out
576
577
578def fetch_tiles(src: Wms, level: Level, t: str | None) -> np.ndarray:
579 px_lon = (level.east - level.west) / level.width
580 px_lat = (level.north - level.south) / level.height
581 rows = [(r, min(r + MAX_TILE, level.height)) for r in range(0, level.height, MAX_TILE)]
582 img = np.zeros((level.height, level.width, 4), np.uint8)
583
584 def one(r0: int, r1: int, c0: int, c1: int) -> None:
585 w = level.west + c0 * px_lon
586 e = level.west + c1 * px_lon
587 mid = (w + e) / 2
588 shift = -360.0 if mid > 180 else 360.0 if mid < -180 else 0.0
589 w, e = max(-180.0, w + shift), min(180.0, e + shift)
590 n = level.north - r0 * px_lat
591 s = level.north - r1 * px_lat
592 tile = fetch_image(src.url((s, w, n, e), c1 - c0, r1 - r0, t))
593 if tile.shape[:2] != (r1 - r0, c1 - c0):
594 tile = np.asarray(Image.fromarray(tile).resize((c1 - c0, r1 - r0), Image.BILINEAR))
595 img[r0:r1, c0:c1] = tile
596
597 with concurrent.futures.ThreadPoolExecutor(max_workers=TILE_WORKERS) as pool:
598 futures = [pool.submit(one, r0, r1, c0, c1)
599 for r0, r1 in rows for c0, c1 in split_columns(level)]
600 for fut in futures:
601 fut.result() # re-raise the first failure
602 return img
603
604
605@dataclass
606class Layer:
607 name: str
608 local: Wms | None
609 global_: Wms
610 prep: Callable[["Layer", np.ndarray, Colormap | None], np.ndarray]
611 combine: Callable[[list[np.ndarray]], np.ndarray] = combine_max
612 colormap: str | None = None # GIBS colormap name
613 transform: Callable[[np.ndarray], np.ndarray] | None = None
614 required: bool = False
615 cmap: Colormap | None = None
616 levels: list[Level] = field(default_factory=list)
617 log_res: np.ndarray | None = None
618
619
620class Cache:
621 """Raw WMS rasters as PNG (refetched after the source TTL) and prepared
622 rasters as .npy (rebuilt when a raw one changed or PREP_VERSION moved)."""
623
624 def __init__(self, root: Path) -> None:
625 self.root = root
626 self.root.mkdir(parents=True, exist_ok=True)
627 self.used: set[Path] = set()
628
629 def _path(self, prefix: str, ident: object, suffix: str) -> Path:
630 digest = hashlib.sha1(json.dumps(ident, sort_keys=True).encode()).hexdigest()[:16]
631 path = self.root / f"{prefix}-{digest}{suffix}"
632 self.used.add(path)
633 return path
634
635 def raw(self, src: Wms, level: Level, t: str | None) -> Path | None:
636 ident = [src.layer, t, level.bbox, level.width, level.height, src.oversample]
637 path = self._path(src.key, ident, ".png")
638 if path.exists() and time.time() - path.stat().st_mtime < src.ttl:
639 return path
640 try:
641 t0 = time.monotonic()
642 img = fetch_level(src, level, t)
643 tmp = path.with_name(f".{path.name}.tmp")
644 Image.fromarray(img, "RGBA").save(tmp, format="PNG", compress_level=1)
645 os.replace(tmp, path)
646 log(f"fetched {src.key} {t or ''} {level.width}x{level.height} "
647 f"@{level.res:.2f}km in {time.monotonic() - t0:.1f}s")
648 except Exception as e: # noqa: BLE001 - any failure falls back to the cache
649 log(f"warning: {src.key} {t or ''}: {e}"
650 + ("; using cached copy" if path.exists() else ""))
651 return path if path.exists() else None
652
653 def get(self, layer: Layer, src: Wms, level: Level, when: dt.datetime) -> np.ndarray | None:
654 raws = [p for t in src.times(when) if (p := self.raw(src, level, t)) is not None]
655 if not raws:
656 return None
657 path = self._path(
658 f"{layer.name}-prep",
659 [PREP_VERSION, src.key, src.gain, [p.name for p in raws]],
660 ".npy",
661 )
662 newest = max(p.stat().st_mtime for p in raws)
663 if not path.exists() or path.stat().st_mtime < newest:
664 try:
665 preps = [
666 layer.prep(layer, apply_gain(np.asarray(Image.open(p).convert("RGBA")), src.gain),
667 layer.cmap)
668 for p in raws
669 ]
670 tmp = path.with_name(f".{path.name}.tmp.npy")
671 np.save(tmp, layer.combine(preps))
672 os.replace(tmp, path)
673 except Exception as e: # noqa: BLE001 - a broken raster drops the level
674 log(f"warning: {layer.name}: preparing {src.key}: {e}")
675 return None
676 return np.load(path, mmap_mode="r")
677
678 def get_url(self, name: str, url: str, ttl: float) -> Path | None:
679 path = self.root / name
680 self.used.add(path)
681 if path.exists() and time.time() - path.stat().st_mtime < ttl:
682 return path
683 try:
684 body, _ = http_get(url, timeout=30)
685 if not body:
686 raise FetchError(f"{url}: empty response")
687 tmp = path.with_name(f".{path.name}.tmp")
688 tmp.write_bytes(body)
689 os.replace(tmp, path)
690 except Exception as e: # noqa: BLE001
691 log(f"warning: {name}: {e}")
692 return path if path.exists() else None
693
694 def cleanup(self, max_age: float = 3 * 86400) -> None:
695 """Delete files no run has used for max_age.
696
697 File mtimes are the fetch times the TTLs run on, so last use is kept
698 in a separate index; that also lets several homes share one cache
699 without deleting each other's local imagery.
700 """
701 index_path = self.root / "last-used.json"
702 try:
703 last_used = json.loads(index_path.read_text())
704 except (OSError, ValueError):
705 last_used = {}
706 now = time.time()
707 for p in self.used:
708 last_used[p.name] = now
709 for p in self.root.iterdir():
710 if p == index_path or p.name.startswith("."):
711 continue
712 seen = last_used.get(p.name, p.stat().st_mtime)
713 if now - seen > max_age:
714 p.unlink()
715 last_used.pop(p.name, None)
716 last_used = {name: t for name, t in last_used.items() if (self.root / name).exists()}
717 tmp = index_path.with_name(f".{index_path.name}.tmp")
718 tmp.write_text(json.dumps(last_used))
719 os.replace(tmp, index_path)
720
721
722def snow_fraction(ndsi: np.ndarray) -> np.ndarray:
723 # MODIS/VIIRS fractional snow cover from NDSI, ignoring the faint
724 # detections that are mostly cloud edges and bright soil
725 return smoothstep(0.35, 0.6, ndsi)
726
727
728def cloud_opacity(ir: np.ndarray) -> np.ndarray:
729 return smoothstep(0.42, 0.78, ir) ** 1.2
730
731
732def layers() -> list[Layer]:
733 day = 86400.0
734 return [
735 Layer(
736 "day",
737 # s2cloudless is bluer than Blue Marble over the same land
738 Wms("s2cloudless", EOX, "s2cloudless-2025", 0.01, 180 * day, "image/jpeg",
739 oversample=4, gain=(1.0, 0.96, 0.85)),
740 Wms("bluemarble", GIBS, "BlueMarble_NextGeneration", 0.5, 180 * day, "image/jpeg"),
741 prep_rgb,
742 required=True,
743 ),
744 Layer(
745 "lights",
746 None,
747 Wms("blackmarble", GIBS, "VIIRS_Black_Marble", 0.5, 180 * day, "image/jpeg",
748 time="2016-01-01"),
749 prep_lights,
750 ),
751 Layer(
752 "clouds",
753 Wms("mtg-ir105", EUMETSAT, "mtg_fd:ir105_hrfi", 1.0, 600),
754 Wms("georing-ir108", EUMETSAT, "mumi:worldcloudmap_ir108", 3.0, 3600),
755 prep_cloud,
756 transform=cloud_opacity,
757 ),
758 Layer(
759 "snow",
760 None,
761 Wms("snow", GIBS, "VIIRS_NOAA20_NDSI_Snow_Cover", 0.375, 6 * 3600,
762 days_back=(1, 2, 3)),
763 prep_scalar,
764 colormap="MODIS_NDSI_Snow_Cover",
765 transform=snow_fraction,
766 ),
767 Layer(
768 "ice",
769 None,
770 Wms("seaice", GIBS, "GHRSST_L4_MUR_Sea_Ice_Concentration", 1.0, 6 * 3600,
771 days_back=(1, 2, 3)),
772 prep_scalar,
773 colormap="GHRSST_Sea_Ice_Concentration",
774 transform=lambda c: smoothstep(0.15, 0.9, c),
775 ),
776 Layer(
777 "fire",
778 None,
779 Wms("fire", GIBS, "VIIRS_NOAA20_Thermal_Anomalies_375m_All", 0.375, 3600,
780 days_back=(0, 1)),
781 prep_scalar,
782 ),
783 Layer(
784 "chlorophyll",
785 None,
786 Wms("chlorophyll", GIBS, "VIIRS_NOAA20_Chlorophyll_a", 1.0, 6 * 3600,
787 days_back=(1, 2, 3, 4)),
788 prep_tint,
789 combine=combine_fill,
790 colormap="VIIRS_Chlorophyll",
791 ),
792 ]
793
794
795def load_layers(proj: LensMap, cache: Cache, when: dt.datetime) -> dict[str, Layer]:
796 all_layers = []
797 for layer in layers():
798 if layer.colormap:
799 path = cache.get_url(f"colormap-{layer.colormap}.xml",
800 f"https://gibs.earthdata.nasa.gov/colormaps/v1.3/{layer.colormap}.xml",
801 30 * 86400)
802 try:
803 layer.cmap = Colormap(path) if path else None
804 except Exception as e: # noqa: BLE001
805 log(f"warning: colormap {layer.colormap}: {e}")
806 if layer.cmap is None:
807 log(f"warning: no colormap for {layer.name}, leaving it out")
808 continue
809 all_layers.append(layer)
810 jobs = []
811 for layer in all_layers:
812 local_src = layer.local or layer.global_
813 for lvl in plan_local(proj, local_src.native_km):
814 jobs.append((layer, local_src, lvl))
815 jobs.append((layer, layer.global_, global_level(GLOBAL_WIDTH)))
816 with concurrent.futures.ThreadPoolExecutor(max_workers=6) as pool:
817 futures = [(layer, lvl, pool.submit(cache.get, layer, src, lvl, when))
818 for layer, src, lvl in jobs]
819 for layer, lvl, fut in futures:
820 lvl.data = fut.result()
821 if lvl.data is not None:
822 layer.levels.append(lvl)
823 out = {}
824 for layer in all_layers:
825 if not layer.levels or not layer.levels[-1].wrap:
826 if layer.required:
827 raise SystemExit(f"no global data for required layer {layer.name}")
828 log(f"warning: layer {layer.name} has no global data, leaving it out")
829 continue
830 layer.levels.extend(mipmaps(layer.levels[-1], proj.q_max))
831 # np.interp needs increasing resolutions; a size-capped local level
832 # can end up coarser than the next one
833 layer.log_res = np.maximum.accumulate(np.log([lvl.res for lvl in layer.levels]))
834 out[layer.name] = layer
835 return out
836
837
838# --------------------------------------------------------------------------
839# sampling
840
841
842def sample_level(level: Level, lat: np.ndarray, lon: np.ndarray):
843 data = level.data
844 assert data is not None
845 h, w = data.shape[:2]
846 col = ((lon - level.west) % 360.0) / (level.east - level.west) * w - 0.5
847 row = (level.north - lat) / (level.north - level.south) * h - 0.5
848 if level.wrap:
849 valid = np.ones(lat.shape, bool)
850 c = np.floor(col)
851 fc = (col - c)[:, None]
852 c0 = c.astype(np.int64) % w
853 c1 = (c0 + 1) % w
854 else:
855 valid = (col >= -0.5) & (col <= w - 0.5) & (row >= -0.5) & (row <= h - 0.5)
856 col = np.clip(col, 0, w - 1)
857 c0 = np.floor(col).astype(np.int64)
858 c1 = np.minimum(c0 + 1, w - 1)
859 fc = (col - c0)[:, None]
860 row = np.clip(row, 0, h - 1)
861 r0 = np.floor(row).astype(np.int64)
862 r1 = np.minimum(r0 + 1, h - 1)
863 fr = (row - r0)[:, None]
864 top = data[r0, c0] * (1 - fc) + data[r0, c1] * fc
865 bottom = data[r1, c0] * (1 - fc) + data[r1, c1] * fc
866 return (top * (1 - fr) + bottom * fr) / 255.0, valid
867
868
869def sample(layer: Layer, lat: np.ndarray, lon: np.ndarray, logq: np.ndarray):
870 """Trilinear sample: blend the two levels whose resolution brackets q.
871
872 Returns premultiplied values and the total data coverage.
873 """
874 n = len(layer.levels)
875 t = np.interp(logq, layer.log_res, np.arange(n, dtype=np.float64))
876 channels = layer.levels[0].data.shape[-1] - 1
877 out = np.zeros((len(lat), channels))
878 rem = np.ones(len(lat))
879 for i, lvl in enumerate(layer.levels):
880 w = np.clip(i + 1 - t, 0.0, 1.0)
881 idx = np.nonzero((w > 0) & (rem > 1e-4))[0]
882 if len(idx) == 0:
883 continue
884 v, valid = sample_level(lvl, lat[idx], lon[idx])
885 wi = rem[idx] * w[idx] * valid
886 out[idx] += v[:, :-1] * wi[:, None]
887 rem[idx] -= wi * v[:, -1]
888 return out, 1.0 - rem
889
890
891def lerp(a: np.ndarray, b, t: np.ndarray) -> np.ndarray:
892 return a + (np.asarray(b) - a) * t[:, None]
893
894
895def subsolar_cos(lat, lon, sun_lat: float, sun_lon: float) -> np.ndarray:
896 la, lo = np.radians(lat), np.radians(lon)
897 sl, so = math.radians(sun_lat), math.radians(sun_lon)
898 return np.sin(la) * math.sin(sl) + np.cos(la) * math.cos(sl) * np.cos(lo - so)
899
900
901# --------------------------------------------------------------------------
902# cloud post-processing: the best live cloud data is ~1 km, far coarser than
903# the zoomed-in pixels around home. Synthesise the missing detail (noise
904# anchored to the ground, so it moves with the clouds, not the screen), and
905# light the clouds from the sun: relief shading plus shadows on the ground.
906
907CLOUD_DATA_KM = 1.0 # resolution of the finest cloud source
908CLOUD_HEIGHT_KM = 3.0 # assumed cloud height for the shadow offset
909
910
911def value_noise(p: np.ndarray) -> np.ndarray:
912 """Smooth 3D value noise in [-1, 1] at integer lattice spacing."""
913 i = np.floor(p).astype(np.int64)
914 f = p - i
915 u = f * f * (3 - 2 * f)
916 out = np.zeros(len(p))
917 for dx in (0, 1):
918 for dy in (0, 1):
919 for dz in (0, 1):
920 h = ((i[:, 0] + dx) * 73856093) ^ ((i[:, 1] + dy) * 19349663) ^ ((i[:, 2] + dz) * 83492791)
921 h = (h * 0x2545F491) & 0xFFFFFF
922 w = (u[:, 0] if dx else 1 - u[:, 0]) * (u[:, 1] if dy else 1 - u[:, 1]) \
923 * (u[:, 2] if dz else 1 - u[:, 2])
924 out += w * (h / 0x7FFFFF - 1.0)
925 return out
926
927
928def cloud_detail(a: np.ndarray, lat: np.ndarray, lon: np.ndarray, q: np.ndarray) -> np.ndarray:
929 """Turn the soft, upsampled cloud edges into crisp fractal ones.
930
931 Octaves (in km on the ground) are only added where they span at least
932 two output px and lie below the data resolution, so the far, small parts
933 of the map stay as they are.
934 """
935 out = a.copy()
936 flat_a, flat_q = a.ravel(), q.ravel()
937 todo = np.nonzero((flat_a > 0.02) & (flat_a < 0.98) & (flat_q < 2 * CLOUD_DATA_KM))[0]
938 octaves = [(0.4, 0.35), (0.8, 0.5), (1.6, 0.75), (3.2, 1.0)] # wavelength km, amplitude
939 total = sum(amp for _, amp in octaves)
940 flat_out = out.ravel()
941 for c0 in range(0, len(todo), 1_000_000):
942 idx = todo[c0 : c0 + 1_000_000]
943 la, lo = np.radians(lat.ravel()[idx]), np.radians(lon.ravel()[idx])
944 p = R_EARTH * np.stack([np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)], -1)
945 qi, ai = flat_q[idx], flat_a[idx]
946 n = np.zeros(len(idx))
947 weight = np.zeros(len(idx))
948 for wavelength, amp in octaves:
949 w = smoothstep(2.0, 4.0, wavelength / qi) * amp
950 n += w * value_noise(p / wavelength)
951 weight += w
952 detail = np.clip(weight / total * 1.5, 0, 1) * smoothstep(2 * CLOUD_DATA_KM, 0.5 * CLOUD_DATA_KM, qi)
953 nn = n / np.maximum(weight, 1e-6) # fBm in [-1, 1]
954 # modulate rather than threshold: thin parts get the most texture,
955 # dense cores only some, and translucency survives
956 textured = np.clip(ai * (1 + 0.9 * nn * (1 - 0.6 * ai)), 0, 1)
957 textured = textured + (smoothstep(0.15, 0.85, textured) - textured) * 0.5
958 flat_out[idx] = ai + (textured - ai) * detail
959 return out
960
961
962def box_blur(a: np.ndarray, r: int) -> np.ndarray:
963 """Approximate Gaussian blur (three box passes per axis)."""
964 if r < 1:
965 return a
966 out = a.astype(np.float32)
967 for axis in (0, 1):
968 for _ in range(3):
969 pad = [(0, 0), (0, 0)]
970 pad[axis] = (r + 1, r)
971 c = np.cumsum(np.pad(out, pad, mode="edge"), axis=axis, dtype=np.float64)
972 hi = np.take(c, np.arange(2 * r + 1, c.shape[axis]), axis=axis)
973 lo = np.take(c, np.arange(0, c.shape[axis] - 2 * r - 1), axis=axis)
974 out = ((hi - lo) / (2 * r + 1)).astype(np.float32)
975 return out
976
977
978def gather(a: np.ndarray, x: np.ndarray, y: np.ndarray) -> np.ndarray:
979 """Bilinear lookup of a 2D array at fractional pixel positions."""
980 h, w = a.shape
981 x = np.clip(x, 0, w - 1.001)
982 y = np.clip(y, 0, h - 1.001)
983 x0, y0 = np.floor(x).astype(np.int64), np.floor(y).astype(np.int64)
984 fx, fy = x - x0, y - y0
985 top = a[y0, x0] * (1 - fx) + a[y0, x0 + 1] * fx
986 bottom = a[y0 + 1, x0] * (1 - fx) + a[y0 + 1, x0 + 1] * fx
987 return top * (1 - fy) + bottom * fy
988
989
990def render_base(proj: LensMap, lay: dict[str, Layer], sun: tuple[float, float]) -> np.ndarray:
991 W, H = proj.width, proj.height
992 ground = np.zeros((H, W, 3), np.float32)
993 cloud = np.zeros((H, W), np.float32)
994 cosz_f = np.zeros((H, W), np.float32)
995 lights = np.zeros((H, W), np.float32)
996 lat_f = np.zeros((H, W), np.float32)
997 lon_f = np.zeros((H, W), np.float32)
998 q_f = np.zeros((H, W), np.float32)
999 chunk = max(1, 1_000_000 // W)
1000 xs = np.arange(W) + 0.5
1001 snow_rgb = np.array([0.87, 0.87, 0.87])
1002 ice_rgb = np.array([0.90, 0.94, 1.00])
1003 fire_rgb = np.array([0.94, 0.52, 0.05])
1004 # chlorophyll: turquoise for moderately productive water, green for blooms
1005 plankton_rgb = np.array([0.06, 0.32, 0.40])
1006 bloom_rgb = np.array([0.14, 0.40, 0.16])
1007 # the night side is the (high resolution) day map in faint moonlight
1008 moonlight = np.array([0.20, 0.24, 0.36])
1009 for y0 in range(0, H, chunk):
1010 ys = np.arange(y0, min(H, y0 + chunk)) + 0.5
1011 rows = slice(y0, y0 + len(ys))
1012 X, Y = np.meshgrid(xs, ys)
1013 lat, lon, q = proj.inverse(X.ravel(), Y.ravel())
1014 logq = np.log(q)
1015
1016 def get(name: str):
1017 return sample(lay[name], lat, lon, logq) if name in lay else None
1018
1019 day, cov = get("day")
1020 day /= np.maximum(cov, 1e-6)[:, None]
1021
1022 if (s := get("chlorophyll")) is not None:
1023 # the value is log-scaled 0.01..20 mg/m³: 0.30 is ~0.1 mg/m³
1024 # (productive shelf water), 0.75 is ~3 mg/m³ (a bloom)
1025 v, _ = s
1026 present = v[:, 1]
1027 value = v[:, 0] / np.maximum(present, 1e-6)
1028 tint = lerp(np.tile(plankton_rgb, (len(value), 1)), bloom_rgb,
1029 smoothstep(0.55, 0.85, value))
1030 day = lerp(day, tint, present * smoothstep(0.30, 0.75, value) * 0.7)
1031 if (s := get("ice")) is not None:
1032 day = lerp(day, ice_rgb, s[0][:, 0] * 0.9)
1033 if (s := get("snow")) is not None:
1034 day = lerp(day, snow_rgb, s[0][:, 0] * 0.85)
1035
1036 cosz = subsolar_cos(lat, lon, *sun)
1037 light = smoothstep(-0.10, 0.08, cosz)
1038 day_lit = day * (0.55 + 0.45 * np.sqrt(np.clip(cosz, 0, 1)))[:, None]
1039 color = lerp(day * moonlight, day_lit, light)
1040 if (s := get("fire")) is not None:
1041 color = np.maximum(color, fire_rgb * 0.8 * s[0][:, 0][:, None])
1042 if (s := get("clouds")) is not None:
1043 cloud[rows] = np.clip(s[0][:, 0], 0, 1).reshape(len(ys), W)
1044 if (s := get("lights")) is not None:
1045 lights[rows] = (s[0][:, 0] * (1 - light)).reshape(len(ys), W)
1046
1047 ground[rows] = color.reshape(len(ys), W, 3)
1048 cosz_f[rows] = cosz.reshape(len(ys), W)
1049 lat_f[rows] = lat.reshape(len(ys), W)
1050 lon_f[rows] = lon.reshape(len(ys), W)
1051 q_f[rows] = q.reshape(len(ys), W)
1052
1053 scale = H / 1800
1054 light = smoothstep(-0.10, 0.08, cosz_f)
1055 if "clouds" in lay:
1056 cloud = cloud_detail(cloud, lat_f, lon_f, q_f)
1057 del lat_f, lon_f
1058
1059 # sun direction on screen: cos(zenith) grows towards the sun
1060 gy, gx = np.gradient(cosz_f)
1061 norm = np.hypot(gx, gy)
1062 lx = np.where(norm > 1e-9, gx / np.maximum(norm, 1e-9), 0.0).astype(np.float32)
1063 ly = np.where(norm > 1e-9, gy / np.maximum(norm, 1e-9), 0.0).astype(np.float32)
1064
1065 # shadows: the cloud seen from the ground towards the sun, displaced by
1066 # height / tan(elevation), only where the sun is up
1067 elev = np.arcsin(np.clip(cosz_f, 0.0, 1.0))
1068 offset = np.clip(CLOUD_HEIGHT_KM / np.maximum(np.tan(elev), 0.1) / q_f, 0, 40 * scale)
1069 yy, xx = np.mgrid[0:H, 0:W].astype(np.float32)
1070 soft = box_blur(cloud, max(1, round(1.5 * scale)))
1071 shadow = gather(soft, xx + lx * offset, yy + ly * offset)
1072 del xx, yy
1073 ground *= (1 - 0.45 * shadow * light * smoothstep(0.0, 0.15, cosz_f))[..., None]
1074
1075 # relief: treat the (blurred) opacity as cloud height and light its
1076 # slopes; a slope facing the sun has its height falling towards the sun
1077 height = box_blur(cloud, max(1, round(2 * scale)))
1078 hy, hx = np.gradient(height)
1079 relief = np.clip(1 - 6.0 * (hx * lx + hy * ly) / scale, 0.6, 1.12)
1080
1081 sunlit = light * (0.6 + 0.4 * np.sqrt(np.clip(cosz_f, 0, 1))) * relief
1082 cloud_rgb = (
1083 np.array([0.11, 0.12, 0.16], np.float32) * (1 - light)[..., None]
1084 + np.array([0.88, 0.89, 0.90], np.float32) * sunlit[..., None]
1085 )
1086 out = ground + (cloud_rgb - ground) * cloud[..., None]
1087
1088 # city lights as a glow drawn in screen space: smooth at any zoom, so the
1089 # ~1 km squares of the source never show; clouds dim them
1090 lights *= 1 - 0.75 * cloud
1091 img = Image.fromarray(to_u8(lights), "L")
1092 core = np.asarray(img.filter(ImageFilter.GaussianBlur(0.9 * scale)), np.float32) / 255
1093 halo = np.asarray(img.filter(ImageFilter.GaussianBlur(6 * scale)), np.float32) / 255
1094 glow = np.minimum(core * 1.1 + halo * 0.8, 1.5)[..., None]
1095 out += glow * np.array([1.0, 0.80, 0.50], np.float32)
1096 return to_u8(out)
1097
1098
1099# --------------------------------------------------------------------------
1100# sky: sun, moon, planets and constellations at their sub-points
1101
1102
1103PLANETS = {
1104 "mercury": "☿",
1105 "venus": "♀",
1106 "mars": "♂",
1107 "jupiter": "♃",
1108 "saturn": "♄",
1109 "uranus": "♅",
1110 "neptune": "♆",
1111}
1112
1113
1114@dataclass
1115class Sky:
1116 bodies: dict[str, tuple[float, float]]
1117 moon_elongation: float
1118 lines: list[tuple[tuple[float, float], tuple[float, float]]]
1119
1120
1121def wrap180(x):
1122 return (np.asarray(x) + 180.0) % 360.0 - 180.0
1123
1124
1125def compute_sky(when: dt.datetime, constellations: Path | None) -> Sky:
1126 warnings.filterwarnings("ignore", module="astropy")
1127 warnings.filterwarnings("ignore", module="erfa")
1128 from astropy.utils import iers
1129
1130 # offline: bundled IERS tables, UT1-UTC drift is far below a pixel
1131 iers.conf.auto_download = False
1132 iers.conf.auto_max_age = None
1133 iers.conf.iers_degraded_accuracy = "warn"
1134 import astropy.units as u
1135 from astropy.coordinates import TETE, SkyCoord, get_body
1136 from astropy.time import Time
1137
1138 t = Time(when)
1139 gast = t.sidereal_time("apparent", "greenwich").deg
1140 frame = TETE(obstime=t)
1141
1142 def subpoints(coord):
1143 c = coord.transform_to(frame)
1144 return np.asarray(c.dec.deg), wrap180(np.asarray(c.ra.deg) - gast)
1145
1146 bodies = {}
1147 raw = {}
1148 for name in ["sun", "moon", *PLANETS]:
1149 raw[name] = get_body(name, t)
1150 la, lo = subpoints(raw[name])
1151 bodies[name] = (float(la), float(lo))
1152 elong = float(raw["sun"].separation(raw["moon"]).deg)
1153
1154 lines = []
1155 if constellations is not None:
1156 feats = json.loads(constellations.read_text())["features"]
1157 strips = [
1158 line for f in feats for line in f["geometry"]["coordinates"] if len(line) > 1
1159 ]
1160 pts = np.array([p for line in strips for p in line], dtype=np.float64)
1161 la, lo = subpoints(SkyCoord(ra=pts[:, 0] * u.deg, dec=pts[:, 1] * u.deg, frame="icrs"))
1162 i = 0
1163 for line in strips:
1164 for k in range(len(line) - 1):
1165 lines.append(((la[i + k], lo[i + k]), (la[i + k + 1], lo[i + k + 1])))
1166 i += len(line)
1167 return Sky(bodies, elong, lines)
1168
1169
1170def to_vec(lat: float, lon: float) -> np.ndarray:
1171 la, lo = math.radians(lat), math.radians(lon)
1172 return np.array([math.cos(la) * math.cos(lo), math.cos(la) * math.sin(lo), math.sin(la)])
1173
1174
1175def great_circle(a, b, proj: LensMap, max_px: float = 3.0) -> list[list[tuple[float, float]]]:
1176 """Project the great circle arc a->b, split where it wraps around the map."""
1177 va, vb = to_vec(*a), to_vec(*b)
1178 omega = math.acos(max(-1.0, min(1.0, float(va @ vb))))
1179 # sample densely, then drop points that add nothing on screen
1180 n = 2
1181 while n < 4096:
1182 ts = np.linspace(0, 1, n)
1183 if omega < 1e-9:
1184 vs = np.outer(1 - ts, va) + np.outer(ts, vb)
1185 else:
1186 vs = (np.outer(np.sin((1 - ts) * omega), va) + np.outer(np.sin(ts * omega), vb)) / math.sin(omega)
1187 lat = np.degrees(np.arcsin(np.clip(vs[:, 2], -1, 1)))
1188 lon = np.degrees(np.arctan2(vs[:, 1], vs[:, 0]))
1189 x, y = proj.forward(lat, lon)
1190 step = np.hypot(np.diff(x), np.diff(y))
1191 if step.max() <= max_px:
1192 break
1193 n *= 2
1194 pieces, cur = [], [(x[0], y[0])]
1195 for k, s in enumerate(step):
1196 if s > proj.width * 0.25:
1197 if len(cur) > 1:
1198 pieces.append(cur)
1199 cur = []
1200 cur.append((x[k + 1], y[k + 1]))
1201 if len(cur) > 1:
1202 pieces.append(cur)
1203 return pieces
1204
1205
1206def sun_icon(path: Path | None, size: int) -> Image.Image:
1207 if path is not None:
1208 try:
1209 img = Image.open(path).convert("RGB").resize((size, size), Image.LANCZOS)
1210 yy, xx = np.mgrid[:size, :size]
1211 r = np.hypot(xx - size / 2 + 0.5, yy - size / 2 + 0.5) / (size / 2)
1212 alpha = np.clip((0.86 - r) / 0.06, 0, 1)
1213 rgba = np.dstack([np.asarray(img), (alpha * 255).astype(np.uint8)])
1214 return Image.fromarray(rgba, "RGBA")
1215 except Exception as e: # noqa: BLE001
1216 log(f"warning: sun image: {e}")
1217 img = Image.new("RGBA", (size, size))
1218 ImageDraw.Draw(img).ellipse((size * 0.2, size * 0.2, size * 0.8, size * 0.8), fill=(255, 210, 80, 255))
1219 return img
1220
1221
1222def moon_icon(size: int, elongation: float, sun_dir: tuple[float, float]) -> Image.Image:
1223 """Moon disc lit from screen direction sun_dir with the real phase."""
1224 phase = math.radians(180.0 - elongation) # sun-moon-earth angle
1225 ex, ey = sun_dir
1226 yy, xx = np.mgrid[:size, :size]
1227 x = (xx + 0.5 - size / 2) / (size / 2)
1228 y = -(yy + 0.5 - size / 2) / (size / 2)
1229 r2 = x * x + y * y
1230 z = np.sqrt(np.clip(1 - r2, 0, 1))
1231 lit = (x * ex + y * ey) * math.sin(phase) + z * math.cos(phase)
1232 shade = np.clip(lit * 4, 0, 1) * 0.85 + 0.12
1233 alpha = np.clip((1 - np.sqrt(r2)) * size / 2, 0, 1)
1234 g = (shade * 235).astype(np.uint8)
1235 rgba = np.dstack([g, g, (shade * 225).astype(np.uint8), (alpha * 255).astype(np.uint8)])
1236 return Image.fromarray(rgba, "RGBA")
1237
1238
1239def paste_centered(canvas: Image.Image, icon: Image.Image, x: float, y: float) -> None:
1240 canvas.alpha_composite(icon, (int(round(x - icon.width / 2)), int(round(y - icon.height / 2))))
1241
1242
1243def draw_lines(canvas: Image.Image, proj: LensMap, sky: Sky, scale: float) -> None:
1244 layer = Image.new("RGBA", canvas.size)
1245 d = ImageDraw.Draw(layer)
1246 width = max(1, round(1.2 * scale))
1247 for a, b in sky.lines:
1248 for piece in great_circle(a, b, proj):
1249 d.line(piece, fill=(255, 255, 255, 110), width=width, joint="curve")
1250 canvas.alpha_composite(layer)
1251
1252
1253def draw_bodies(canvas, proj, sky: Sky, scale: float, font_path, sun_image) -> None:
1254 sx, sy = proj.forward(*sky.bodies["sun"])
1255 mx, my = proj.forward(*sky.bodies["moon"])
1256 font = ImageFont.truetype(str(font_path), max(8, round(34 * scale))) if font_path else None
1257 d = ImageDraw.Draw(canvas)
1258 for name, glyph in PLANETS.items():
1259 x, y = proj.forward(*sky.bodies[name])
1260 if font is not None:
1261 d.text((float(x), float(y)), glyph, font=font, anchor="mm", fill=(254, 128, 25, 255),
1262 stroke_width=max(1, round(2 * scale)), stroke_fill=(0, 0, 0, 200))
1263 else:
1264 r = 6 * scale
1265 d.ellipse((x - r, y - r, x + r, y + r), fill=(254, 128, 25, 255))
1266 dx, dy = float(sx - mx), float(my - sy)
1267 norm = math.hypot(dx, dy) or 1.0
1268 paste_centered(canvas, moon_icon(max(8, round(50 * scale)), sky.moon_elongation, (dx / norm, dy / norm)),
1269 float(mx), float(my))
1270 paste_centered(canvas, sun_icon(sun_image, max(8, round(100 * scale))), float(sx), float(sy))
1271
1272
1273def draw_markers(canvas, proj, markers: list[tuple[float, float]], icon_path, scale: float) -> None:
1274 if not markers or icon_path is None:
1275 return
1276 size = max(8, round(20 * scale))
1277 icon = Image.open(icon_path).convert("RGBA").resize((size, size), Image.LANCZOS)
1278 for lat, lon in markers:
1279 x, y = proj.forward(lat, lon)
1280 paste_centered(canvas, icon, float(x), float(y))
1281
1282
1283def load_markers(path: Path | None) -> list[tuple[float, float]]:
1284 if path is None:
1285 return []
1286 try:
1287 data = json.loads(path.read_text())
1288 except Exception as e: # noqa: BLE001
1289 log(f"warning: markers: {e}")
1290 return []
1291 entries = data.values() if isinstance(data, dict) else data
1292 out = []
1293 for v in entries:
1294 if isinstance(v, dict) and v.get("latitude") is not None and v.get("longitude") is not None:
1295 out.append((float(v["latitude"]), float(v["longitude"])))
1296 return out
1297
1298
1299def save(img: Image.Image, path: Path, **kw) -> None:
1300 tmp = path.with_name(f".{path.name}.tmp")
1301 img.save(tmp, format=path.suffix.lstrip(".").replace("jpg", "jpeg").upper(), **kw)
1302 os.replace(tmp, path)
1303
1304
1305# --------------------------------------------------------------------------
1306
1307
1308def main() -> None:
1309 env = os.environ.get
1310 p = argparse.ArgumentParser(prog="realwallpaper", description=__doc__.splitlines()[0])
1311 p.add_argument("--dir", type=Path, default=Path("."),
1312 help="output and cache directory (default: .)")
1313 p.add_argument("--cache-dir", type=Path, default=None,
1314 help="downloaded data, may be shared between runs for different homes "
1315 "(default: DIR/cache)")
1316 p.add_argument("--lat", type=float, default=52.520008,
1317 help="home latitude in degrees, south negative (default Berlin)")
1318 p.add_argument("--lon", type=float, default=13.404954,
1319 help="home longitude in degrees, west negative (default Berlin)")
1320 p.add_argument("--width", type=int, default=3840)
1321 p.add_argument("--height", type=int, default=2160)
1322 p.add_argument("--magnification", type=float, default=6.0,
1323 help="zoom at home relative to the plain map (default 6)")
1324 p.add_argument("--lens", type=float, default=1.0,
1325 help="lens size relative to the screen; 1 reaches the edges (default 1)")
1326 p.add_argument("--lens-shape", choices=["circle", "ellipse", "rect"], default="ellipse",
1327 help="circle; ellipse or rect follow the screen's aspect (default ellipse)")
1328 p.add_argument("--focus", type=float, default=0.0,
1329 help="share of the lens zoomed undistorted before the transition (default 0)")
1330 p.add_argument("--falloff", type=float, default=0.65,
1331 help="where across the lens (0..1) the distortion peaks; higher keeps "
1332 "the middle closer to a plain zoom (default 0.65)")
1333 p.add_argument("--projection", choices=["lens", "plain", "mercator"], default="lens",
1334 help="lens: centred on home with the zoom lens; plain: centred on home "
1335 "without distortion; mercator: classic north-up map (default lens)")
1336 p.add_argument("--marker-url", default=None,
1337 help="JSON with {name: {latitude, longitude}} for krebs markers")
1338 p.add_argument("--time", type=dt.datetime.fromisoformat, default=None,
1339 help="render for this UTC time instead of now (ISO 8601)")
1340 p.add_argument("--no-archive", action="store_true", help="do not keep a jpg in archive/")
1341 args = p.parse_args()
1342
1343 when = args.time or dt.datetime.now(dt.timezone.utc)
1344 if when.tzinfo is None:
1345 when = when.replace(tzinfo=dt.timezone.utc)
1346 when = when.astimezone(dt.timezone.utc)
1347
1348 out_dir = args.dir.resolve()
1349 out_dir.mkdir(parents=True, exist_ok=True)
1350 cache = Cache(args.cache_dir.resolve() if args.cache_dir else out_dir / "cache")
1351 proj = LensMap(args.lat, args.lon, args.width, args.height,
1352 args.magnification, args.lens, args.lens_shape, args.focus, args.falloff,
1353 args.projection)
1354 home_x, home_y = proj.forward(args.lat, args.lon)
1355 home_q = proj.inverse(np.array([home_x]), np.array([home_y]))[2][0]
1356 log(f"projection: {args.projection}, home at {home_q * 1000:.0f} m/px, "
1357 f"finest {proj.q0 * 1000:.0f} m/px")
1358
1359 font = env("REALWALLPAPER_FONT")
1360 constellations = env("REALWALLPAPER_CONSTELLATIONS")
1361 krebs_icon = env("REALWALLPAPER_KREBS_ICON")
1362
1363 t0 = time.monotonic()
1364 sun_image = cache.get_url("sun.jpg", "https://sdo.gsfc.nasa.gov/assets/img/latest/latest_512_0171.jpg", 3600)
1365 marker_file = cache.get_url("markers.json", args.marker_url, 0) if args.marker_url else None
1366 lay = load_layers(proj, cache, when)
1367 log(f"data ready in {time.monotonic() - t0:.1f}s")
1368
1369 try:
1370 sky = compute_sky(when, Path(constellations) if constellations else None)
1371 except Exception as e: # noqa: BLE001 - the map is still useful without the sky
1372 log(f"warning: sky: {e}")
1373 sky = None
1374 sun = sky.bodies["sun"] if sky else (0.0, 0.0)
1375
1376 t0 = time.monotonic()
1377 base = Image.fromarray(render_base(proj, lay, sun)).convert("RGBA")
1378 log(f"rendered in {time.monotonic() - t0:.1f}s")
1379
1380 scale = args.height / 1800
1381 markers = load_markers(marker_file)
1382
1383 plain = base
1384 marker = plain.copy()
1385 if sky:
1386 draw_bodies(marker, proj, sky, scale, font, sun_image)
1387 krebs = marker.copy()
1388 draw_markers(krebs, proj, markers, krebs_icon, scale)
1389 stars = plain.copy()
1390 if sky:
1391 draw_lines(stars, proj, sky, scale)
1392 draw_bodies(stars, proj, sky, scale, font, sun_image)
1393 draw_markers(stars, proj, markers, krebs_icon, scale)
1394
1395 for name, img in [
1396 ("realwallpaper.png", plain),
1397 ("realwallpaper-marker.png", marker),
1398 ("realwallpaper-krebs.png", krebs),
1399 ("realwallpaper-krebs-stars.png", stars),
1400 ]:
1401 save(img.convert("RGB"), out_dir / name, compress_level=3)
1402 if not args.no_archive:
1403 (out_dir / "archive").mkdir(exist_ok=True)
1404 save(krebs.convert("RGB"), out_dir / "archive" / f"{when.strftime('%Y-%m-%dT%H:%M:%SZ')}.jpg",
1405 quality=88)
1406 cache.cleanup()
1407 log("done")
1408
1409
1410if __name__ == "__main__":
1411 main()