This repository has no description
56 kB
1386 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 for p in self.root.iterdir():
696 if p not in self.used and time.time() - p.stat().st_mtime > max_age:
697 p.unlink()
698
699
700def snow_fraction(ndsi: np.ndarray) -> np.ndarray:
701 # MODIS/VIIRS fractional snow cover from NDSI, ignoring the faint
702 # detections that are mostly cloud edges and bright soil
703 return smoothstep(0.35, 0.6, ndsi)
704
705
706def cloud_opacity(ir: np.ndarray) -> np.ndarray:
707 return smoothstep(0.42, 0.78, ir) ** 1.2
708
709
710def layers() -> list[Layer]:
711 day = 86400.0
712 return [
713 Layer(
714 "day",
715 # s2cloudless is bluer than Blue Marble over the same land
716 Wms("s2cloudless", EOX, "s2cloudless-2025", 0.01, 180 * day, "image/jpeg",
717 oversample=4, gain=(1.0, 0.96, 0.85)),
718 Wms("bluemarble", GIBS, "BlueMarble_NextGeneration", 0.5, 180 * day, "image/jpeg"),
719 prep_rgb,
720 required=True,
721 ),
722 Layer(
723 "lights",
724 None,
725 Wms("blackmarble", GIBS, "VIIRS_Black_Marble", 0.5, 180 * day, "image/jpeg",
726 time="2016-01-01"),
727 prep_lights,
728 ),
729 Layer(
730 "clouds",
731 Wms("mtg-ir105", EUMETSAT, "mtg_fd:ir105_hrfi", 1.0, 600),
732 Wms("georing-ir108", EUMETSAT, "mumi:worldcloudmap_ir108", 3.0, 3600),
733 prep_cloud,
734 transform=cloud_opacity,
735 ),
736 Layer(
737 "snow",
738 None,
739 Wms("snow", GIBS, "VIIRS_NOAA20_NDSI_Snow_Cover", 0.375, 6 * 3600,
740 days_back=(1, 2, 3)),
741 prep_scalar,
742 colormap="MODIS_NDSI_Snow_Cover",
743 transform=snow_fraction,
744 ),
745 Layer(
746 "ice",
747 None,
748 Wms("seaice", GIBS, "GHRSST_L4_MUR_Sea_Ice_Concentration", 1.0, 6 * 3600,
749 days_back=(1, 2, 3)),
750 prep_scalar,
751 colormap="GHRSST_Sea_Ice_Concentration",
752 transform=lambda c: smoothstep(0.15, 0.9, c),
753 ),
754 Layer(
755 "fire",
756 None,
757 Wms("fire", GIBS, "VIIRS_NOAA20_Thermal_Anomalies_375m_All", 0.375, 3600,
758 days_back=(0, 1)),
759 prep_scalar,
760 ),
761 Layer(
762 "chlorophyll",
763 None,
764 Wms("chlorophyll", GIBS, "VIIRS_NOAA20_Chlorophyll_a", 1.0, 6 * 3600,
765 days_back=(1, 2, 3, 4)),
766 prep_tint,
767 combine=combine_fill,
768 colormap="VIIRS_Chlorophyll",
769 ),
770 ]
771
772
773def load_layers(proj: LensMap, cache: Cache, when: dt.datetime) -> dict[str, Layer]:
774 all_layers = []
775 for layer in layers():
776 if layer.colormap:
777 path = cache.get_url(f"colormap-{layer.colormap}.xml",
778 f"https://gibs.earthdata.nasa.gov/colormaps/v1.3/{layer.colormap}.xml",
779 30 * 86400)
780 try:
781 layer.cmap = Colormap(path) if path else None
782 except Exception as e: # noqa: BLE001
783 log(f"warning: colormap {layer.colormap}: {e}")
784 if layer.cmap is None:
785 log(f"warning: no colormap for {layer.name}, leaving it out")
786 continue
787 all_layers.append(layer)
788 jobs = []
789 for layer in all_layers:
790 local_src = layer.local or layer.global_
791 for lvl in plan_local(proj, local_src.native_km):
792 jobs.append((layer, local_src, lvl))
793 jobs.append((layer, layer.global_, global_level(GLOBAL_WIDTH)))
794 with concurrent.futures.ThreadPoolExecutor(max_workers=6) as pool:
795 futures = [(layer, lvl, pool.submit(cache.get, layer, src, lvl, when))
796 for layer, src, lvl in jobs]
797 for layer, lvl, fut in futures:
798 lvl.data = fut.result()
799 if lvl.data is not None:
800 layer.levels.append(lvl)
801 out = {}
802 for layer in all_layers:
803 if not layer.levels or not layer.levels[-1].wrap:
804 if layer.required:
805 raise SystemExit(f"no global data for required layer {layer.name}")
806 log(f"warning: layer {layer.name} has no global data, leaving it out")
807 continue
808 layer.levels.extend(mipmaps(layer.levels[-1], proj.q_max))
809 # np.interp needs increasing resolutions; a size-capped local level
810 # can end up coarser than the next one
811 layer.log_res = np.maximum.accumulate(np.log([lvl.res for lvl in layer.levels]))
812 out[layer.name] = layer
813 return out
814
815
816# --------------------------------------------------------------------------
817# sampling
818
819
820def sample_level(level: Level, lat: np.ndarray, lon: np.ndarray):
821 data = level.data
822 assert data is not None
823 h, w = data.shape[:2]
824 col = ((lon - level.west) % 360.0) / (level.east - level.west) * w - 0.5
825 row = (level.north - lat) / (level.north - level.south) * h - 0.5
826 if level.wrap:
827 valid = np.ones(lat.shape, bool)
828 c = np.floor(col)
829 fc = (col - c)[:, None]
830 c0 = c.astype(np.int64) % w
831 c1 = (c0 + 1) % w
832 else:
833 valid = (col >= -0.5) & (col <= w - 0.5) & (row >= -0.5) & (row <= h - 0.5)
834 col = np.clip(col, 0, w - 1)
835 c0 = np.floor(col).astype(np.int64)
836 c1 = np.minimum(c0 + 1, w - 1)
837 fc = (col - c0)[:, None]
838 row = np.clip(row, 0, h - 1)
839 r0 = np.floor(row).astype(np.int64)
840 r1 = np.minimum(r0 + 1, h - 1)
841 fr = (row - r0)[:, None]
842 top = data[r0, c0] * (1 - fc) + data[r0, c1] * fc
843 bottom = data[r1, c0] * (1 - fc) + data[r1, c1] * fc
844 return (top * (1 - fr) + bottom * fr) / 255.0, valid
845
846
847def sample(layer: Layer, lat: np.ndarray, lon: np.ndarray, logq: np.ndarray):
848 """Trilinear sample: blend the two levels whose resolution brackets q.
849
850 Returns premultiplied values and the total data coverage.
851 """
852 n = len(layer.levels)
853 t = np.interp(logq, layer.log_res, np.arange(n, dtype=np.float64))
854 channels = layer.levels[0].data.shape[-1] - 1
855 out = np.zeros((len(lat), channels))
856 rem = np.ones(len(lat))
857 for i, lvl in enumerate(layer.levels):
858 w = np.clip(i + 1 - t, 0.0, 1.0)
859 idx = np.nonzero((w > 0) & (rem > 1e-4))[0]
860 if len(idx) == 0:
861 continue
862 v, valid = sample_level(lvl, lat[idx], lon[idx])
863 wi = rem[idx] * w[idx] * valid
864 out[idx] += v[:, :-1] * wi[:, None]
865 rem[idx] -= wi * v[:, -1]
866 return out, 1.0 - rem
867
868
869def lerp(a: np.ndarray, b, t: np.ndarray) -> np.ndarray:
870 return a + (np.asarray(b) - a) * t[:, None]
871
872
873def subsolar_cos(lat, lon, sun_lat: float, sun_lon: float) -> np.ndarray:
874 la, lo = np.radians(lat), np.radians(lon)
875 sl, so = math.radians(sun_lat), math.radians(sun_lon)
876 return np.sin(la) * math.sin(sl) + np.cos(la) * math.cos(sl) * np.cos(lo - so)
877
878
879# --------------------------------------------------------------------------
880# cloud post-processing: the best live cloud data is ~1 km, far coarser than
881# the zoomed-in pixels around home. Synthesise the missing detail (noise
882# anchored to the ground, so it moves with the clouds, not the screen), and
883# light the clouds from the sun: relief shading plus shadows on the ground.
884
885CLOUD_DATA_KM = 1.0 # resolution of the finest cloud source
886CLOUD_HEIGHT_KM = 3.0 # assumed cloud height for the shadow offset
887
888
889def value_noise(p: np.ndarray) -> np.ndarray:
890 """Smooth 3D value noise in [-1, 1] at integer lattice spacing."""
891 i = np.floor(p).astype(np.int64)
892 f = p - i
893 u = f * f * (3 - 2 * f)
894 out = np.zeros(len(p))
895 for dx in (0, 1):
896 for dy in (0, 1):
897 for dz in (0, 1):
898 h = ((i[:, 0] + dx) * 73856093) ^ ((i[:, 1] + dy) * 19349663) ^ ((i[:, 2] + dz) * 83492791)
899 h = (h * 0x2545F491) & 0xFFFFFF
900 w = (u[:, 0] if dx else 1 - u[:, 0]) * (u[:, 1] if dy else 1 - u[:, 1]) \
901 * (u[:, 2] if dz else 1 - u[:, 2])
902 out += w * (h / 0x7FFFFF - 1.0)
903 return out
904
905
906def cloud_detail(a: np.ndarray, lat: np.ndarray, lon: np.ndarray, q: np.ndarray) -> np.ndarray:
907 """Turn the soft, upsampled cloud edges into crisp fractal ones.
908
909 Octaves (in km on the ground) are only added where they span at least
910 two output px and lie below the data resolution, so the far, small parts
911 of the map stay as they are.
912 """
913 out = a.copy()
914 flat_a, flat_q = a.ravel(), q.ravel()
915 todo = np.nonzero((flat_a > 0.02) & (flat_a < 0.98) & (flat_q < 2 * CLOUD_DATA_KM))[0]
916 octaves = [(0.4, 0.35), (0.8, 0.5), (1.6, 0.75), (3.2, 1.0)] # wavelength km, amplitude
917 total = sum(amp for _, amp in octaves)
918 flat_out = out.ravel()
919 for c0 in range(0, len(todo), 1_000_000):
920 idx = todo[c0 : c0 + 1_000_000]
921 la, lo = np.radians(lat.ravel()[idx]), np.radians(lon.ravel()[idx])
922 p = R_EARTH * np.stack([np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)], -1)
923 qi, ai = flat_q[idx], flat_a[idx]
924 n = np.zeros(len(idx))
925 weight = np.zeros(len(idx))
926 for wavelength, amp in octaves:
927 w = smoothstep(2.0, 4.0, wavelength / qi) * amp
928 n += w * value_noise(p / wavelength)
929 weight += w
930 detail = np.clip(weight / total * 1.5, 0, 1) * smoothstep(2 * CLOUD_DATA_KM, 0.5 * CLOUD_DATA_KM, qi)
931 nn = n / np.maximum(weight, 1e-6) # fBm in [-1, 1]
932 # modulate rather than threshold: thin parts get the most texture,
933 # dense cores only some, and translucency survives
934 textured = np.clip(ai * (1 + 0.9 * nn * (1 - 0.6 * ai)), 0, 1)
935 textured = textured + (smoothstep(0.15, 0.85, textured) - textured) * 0.5
936 flat_out[idx] = ai + (textured - ai) * detail
937 return out
938
939
940def box_blur(a: np.ndarray, r: int) -> np.ndarray:
941 """Approximate Gaussian blur (three box passes per axis)."""
942 if r < 1:
943 return a
944 out = a.astype(np.float32)
945 for axis in (0, 1):
946 for _ in range(3):
947 pad = [(0, 0), (0, 0)]
948 pad[axis] = (r + 1, r)
949 c = np.cumsum(np.pad(out, pad, mode="edge"), axis=axis, dtype=np.float64)
950 hi = np.take(c, np.arange(2 * r + 1, c.shape[axis]), axis=axis)
951 lo = np.take(c, np.arange(0, c.shape[axis] - 2 * r - 1), axis=axis)
952 out = ((hi - lo) / (2 * r + 1)).astype(np.float32)
953 return out
954
955
956def gather(a: np.ndarray, x: np.ndarray, y: np.ndarray) -> np.ndarray:
957 """Bilinear lookup of a 2D array at fractional pixel positions."""
958 h, w = a.shape
959 x = np.clip(x, 0, w - 1.001)
960 y = np.clip(y, 0, h - 1.001)
961 x0, y0 = np.floor(x).astype(np.int64), np.floor(y).astype(np.int64)
962 fx, fy = x - x0, y - y0
963 top = a[y0, x0] * (1 - fx) + a[y0, x0 + 1] * fx
964 bottom = a[y0 + 1, x0] * (1 - fx) + a[y0 + 1, x0 + 1] * fx
965 return top * (1 - fy) + bottom * fy
966
967
968def render_base(proj: LensMap, lay: dict[str, Layer], sun: tuple[float, float]) -> np.ndarray:
969 W, H = proj.width, proj.height
970 ground = np.zeros((H, W, 3), np.float32)
971 cloud = np.zeros((H, W), np.float32)
972 cosz_f = np.zeros((H, W), np.float32)
973 lights = np.zeros((H, W), np.float32)
974 lat_f = np.zeros((H, W), np.float32)
975 lon_f = np.zeros((H, W), np.float32)
976 q_f = np.zeros((H, W), np.float32)
977 chunk = max(1, 1_000_000 // W)
978 xs = np.arange(W) + 0.5
979 snow_rgb = np.array([0.87, 0.87, 0.87])
980 ice_rgb = np.array([0.90, 0.94, 1.00])
981 fire_rgb = np.array([0.94, 0.52, 0.05])
982 # chlorophyll: turquoise for moderately productive water, green for blooms
983 plankton_rgb = np.array([0.06, 0.32, 0.40])
984 bloom_rgb = np.array([0.14, 0.40, 0.16])
985 # the night side is the (high resolution) day map in faint moonlight
986 moonlight = np.array([0.20, 0.24, 0.36])
987 for y0 in range(0, H, chunk):
988 ys = np.arange(y0, min(H, y0 + chunk)) + 0.5
989 rows = slice(y0, y0 + len(ys))
990 X, Y = np.meshgrid(xs, ys)
991 lat, lon, q = proj.inverse(X.ravel(), Y.ravel())
992 logq = np.log(q)
993
994 def get(name: str):
995 return sample(lay[name], lat, lon, logq) if name in lay else None
996
997 day, cov = get("day")
998 day /= np.maximum(cov, 1e-6)[:, None]
999
1000 if (s := get("chlorophyll")) is not None:
1001 # the value is log-scaled 0.01..20 mg/m³: 0.30 is ~0.1 mg/m³
1002 # (productive shelf water), 0.75 is ~3 mg/m³ (a bloom)
1003 v, _ = s
1004 present = v[:, 1]
1005 value = v[:, 0] / np.maximum(present, 1e-6)
1006 tint = lerp(np.tile(plankton_rgb, (len(value), 1)), bloom_rgb,
1007 smoothstep(0.55, 0.85, value))
1008 day = lerp(day, tint, present * smoothstep(0.30, 0.75, value) * 0.7)
1009 if (s := get("ice")) is not None:
1010 day = lerp(day, ice_rgb, s[0][:, 0] * 0.9)
1011 if (s := get("snow")) is not None:
1012 day = lerp(day, snow_rgb, s[0][:, 0] * 0.85)
1013
1014 cosz = subsolar_cos(lat, lon, *sun)
1015 light = smoothstep(-0.10, 0.08, cosz)
1016 day_lit = day * (0.55 + 0.45 * np.sqrt(np.clip(cosz, 0, 1)))[:, None]
1017 color = lerp(day * moonlight, day_lit, light)
1018 if (s := get("fire")) is not None:
1019 color = np.maximum(color, fire_rgb * 0.8 * s[0][:, 0][:, None])
1020 if (s := get("clouds")) is not None:
1021 cloud[rows] = np.clip(s[0][:, 0], 0, 1).reshape(len(ys), W)
1022 if (s := get("lights")) is not None:
1023 lights[rows] = (s[0][:, 0] * (1 - light)).reshape(len(ys), W)
1024
1025 ground[rows] = color.reshape(len(ys), W, 3)
1026 cosz_f[rows] = cosz.reshape(len(ys), W)
1027 lat_f[rows] = lat.reshape(len(ys), W)
1028 lon_f[rows] = lon.reshape(len(ys), W)
1029 q_f[rows] = q.reshape(len(ys), W)
1030
1031 scale = H / 1800
1032 light = smoothstep(-0.10, 0.08, cosz_f)
1033 if "clouds" in lay:
1034 cloud = cloud_detail(cloud, lat_f, lon_f, q_f)
1035 del lat_f, lon_f
1036
1037 # sun direction on screen: cos(zenith) grows towards the sun
1038 gy, gx = np.gradient(cosz_f)
1039 norm = np.hypot(gx, gy)
1040 lx = np.where(norm > 1e-9, gx / np.maximum(norm, 1e-9), 0.0).astype(np.float32)
1041 ly = np.where(norm > 1e-9, gy / np.maximum(norm, 1e-9), 0.0).astype(np.float32)
1042
1043 # shadows: the cloud seen from the ground towards the sun, displaced by
1044 # height / tan(elevation), only where the sun is up
1045 elev = np.arcsin(np.clip(cosz_f, 0.0, 1.0))
1046 offset = np.clip(CLOUD_HEIGHT_KM / np.maximum(np.tan(elev), 0.1) / q_f, 0, 40 * scale)
1047 yy, xx = np.mgrid[0:H, 0:W].astype(np.float32)
1048 soft = box_blur(cloud, max(1, round(1.5 * scale)))
1049 shadow = gather(soft, xx + lx * offset, yy + ly * offset)
1050 del xx, yy
1051 ground *= (1 - 0.45 * shadow * light * smoothstep(0.0, 0.15, cosz_f))[..., None]
1052
1053 # relief: treat the (blurred) opacity as cloud height and light its
1054 # slopes; a slope facing the sun has its height falling towards the sun
1055 height = box_blur(cloud, max(1, round(2 * scale)))
1056 hy, hx = np.gradient(height)
1057 relief = np.clip(1 - 6.0 * (hx * lx + hy * ly) / scale, 0.6, 1.12)
1058
1059 sunlit = light * (0.6 + 0.4 * np.sqrt(np.clip(cosz_f, 0, 1))) * relief
1060 cloud_rgb = (
1061 np.array([0.11, 0.12, 0.16], np.float32) * (1 - light)[..., None]
1062 + np.array([0.88, 0.89, 0.90], np.float32) * sunlit[..., None]
1063 )
1064 out = ground + (cloud_rgb - ground) * cloud[..., None]
1065
1066 # city lights as a glow drawn in screen space: smooth at any zoom, so the
1067 # ~1 km squares of the source never show; clouds dim them
1068 lights *= 1 - 0.75 * cloud
1069 img = Image.fromarray(to_u8(lights), "L")
1070 core = np.asarray(img.filter(ImageFilter.GaussianBlur(0.9 * scale)), np.float32) / 255
1071 halo = np.asarray(img.filter(ImageFilter.GaussianBlur(6 * scale)), np.float32) / 255
1072 glow = np.minimum(core * 1.1 + halo * 0.8, 1.5)[..., None]
1073 out += glow * np.array([1.0, 0.80, 0.50], np.float32)
1074 return to_u8(out)
1075
1076
1077# --------------------------------------------------------------------------
1078# sky: sun, moon, planets and constellations at their sub-points
1079
1080
1081PLANETS = {
1082 "mercury": "☿",
1083 "venus": "♀",
1084 "mars": "♂",
1085 "jupiter": "♃",
1086 "saturn": "♄",
1087 "uranus": "♅",
1088 "neptune": "♆",
1089}
1090
1091
1092@dataclass
1093class Sky:
1094 bodies: dict[str, tuple[float, float]]
1095 moon_elongation: float
1096 lines: list[tuple[tuple[float, float], tuple[float, float]]]
1097
1098
1099def wrap180(x):
1100 return (np.asarray(x) + 180.0) % 360.0 - 180.0
1101
1102
1103def compute_sky(when: dt.datetime, constellations: Path | None) -> Sky:
1104 warnings.filterwarnings("ignore", module="astropy")
1105 warnings.filterwarnings("ignore", module="erfa")
1106 from astropy.utils import iers
1107
1108 # offline: bundled IERS tables, UT1-UTC drift is far below a pixel
1109 iers.conf.auto_download = False
1110 iers.conf.auto_max_age = None
1111 iers.conf.iers_degraded_accuracy = "warn"
1112 import astropy.units as u
1113 from astropy.coordinates import TETE, SkyCoord, get_body
1114 from astropy.time import Time
1115
1116 t = Time(when)
1117 gast = t.sidereal_time("apparent", "greenwich").deg
1118 frame = TETE(obstime=t)
1119
1120 def subpoints(coord):
1121 c = coord.transform_to(frame)
1122 return np.asarray(c.dec.deg), wrap180(np.asarray(c.ra.deg) - gast)
1123
1124 bodies = {}
1125 raw = {}
1126 for name in ["sun", "moon", *PLANETS]:
1127 raw[name] = get_body(name, t)
1128 la, lo = subpoints(raw[name])
1129 bodies[name] = (float(la), float(lo))
1130 elong = float(raw["sun"].separation(raw["moon"]).deg)
1131
1132 lines = []
1133 if constellations is not None:
1134 feats = json.loads(constellations.read_text())["features"]
1135 strips = [
1136 line for f in feats for line in f["geometry"]["coordinates"] if len(line) > 1
1137 ]
1138 pts = np.array([p for line in strips for p in line], dtype=np.float64)
1139 la, lo = subpoints(SkyCoord(ra=pts[:, 0] * u.deg, dec=pts[:, 1] * u.deg, frame="icrs"))
1140 i = 0
1141 for line in strips:
1142 for k in range(len(line) - 1):
1143 lines.append(((la[i + k], lo[i + k]), (la[i + k + 1], lo[i + k + 1])))
1144 i += len(line)
1145 return Sky(bodies, elong, lines)
1146
1147
1148def to_vec(lat: float, lon: float) -> np.ndarray:
1149 la, lo = math.radians(lat), math.radians(lon)
1150 return np.array([math.cos(la) * math.cos(lo), math.cos(la) * math.sin(lo), math.sin(la)])
1151
1152
1153def great_circle(a, b, proj: LensMap, max_px: float = 3.0) -> list[list[tuple[float, float]]]:
1154 """Project the great circle arc a->b, split where it wraps around the map."""
1155 va, vb = to_vec(*a), to_vec(*b)
1156 omega = math.acos(max(-1.0, min(1.0, float(va @ vb))))
1157 # sample densely, then drop points that add nothing on screen
1158 n = 2
1159 while n < 4096:
1160 ts = np.linspace(0, 1, n)
1161 if omega < 1e-9:
1162 vs = np.outer(1 - ts, va) + np.outer(ts, vb)
1163 else:
1164 vs = (np.outer(np.sin((1 - ts) * omega), va) + np.outer(np.sin(ts * omega), vb)) / math.sin(omega)
1165 lat = np.degrees(np.arcsin(np.clip(vs[:, 2], -1, 1)))
1166 lon = np.degrees(np.arctan2(vs[:, 1], vs[:, 0]))
1167 x, y = proj.forward(lat, lon)
1168 step = np.hypot(np.diff(x), np.diff(y))
1169 if step.max() <= max_px:
1170 break
1171 n *= 2
1172 pieces, cur = [], [(x[0], y[0])]
1173 for k, s in enumerate(step):
1174 if s > proj.width * 0.25:
1175 if len(cur) > 1:
1176 pieces.append(cur)
1177 cur = []
1178 cur.append((x[k + 1], y[k + 1]))
1179 if len(cur) > 1:
1180 pieces.append(cur)
1181 return pieces
1182
1183
1184def sun_icon(path: Path | None, size: int) -> Image.Image:
1185 if path is not None:
1186 try:
1187 img = Image.open(path).convert("RGB").resize((size, size), Image.LANCZOS)
1188 yy, xx = np.mgrid[:size, :size]
1189 r = np.hypot(xx - size / 2 + 0.5, yy - size / 2 + 0.5) / (size / 2)
1190 alpha = np.clip((0.86 - r) / 0.06, 0, 1)
1191 rgba = np.dstack([np.asarray(img), (alpha * 255).astype(np.uint8)])
1192 return Image.fromarray(rgba, "RGBA")
1193 except Exception as e: # noqa: BLE001
1194 log(f"warning: sun image: {e}")
1195 img = Image.new("RGBA", (size, size))
1196 ImageDraw.Draw(img).ellipse((size * 0.2, size * 0.2, size * 0.8, size * 0.8), fill=(255, 210, 80, 255))
1197 return img
1198
1199
1200def moon_icon(size: int, elongation: float, sun_dir: tuple[float, float]) -> Image.Image:
1201 """Moon disc lit from screen direction sun_dir with the real phase."""
1202 phase = math.radians(180.0 - elongation) # sun-moon-earth angle
1203 ex, ey = sun_dir
1204 yy, xx = np.mgrid[:size, :size]
1205 x = (xx + 0.5 - size / 2) / (size / 2)
1206 y = -(yy + 0.5 - size / 2) / (size / 2)
1207 r2 = x * x + y * y
1208 z = np.sqrt(np.clip(1 - r2, 0, 1))
1209 lit = (x * ex + y * ey) * math.sin(phase) + z * math.cos(phase)
1210 shade = np.clip(lit * 4, 0, 1) * 0.85 + 0.12
1211 alpha = np.clip((1 - np.sqrt(r2)) * size / 2, 0, 1)
1212 g = (shade * 235).astype(np.uint8)
1213 rgba = np.dstack([g, g, (shade * 225).astype(np.uint8), (alpha * 255).astype(np.uint8)])
1214 return Image.fromarray(rgba, "RGBA")
1215
1216
1217def paste_centered(canvas: Image.Image, icon: Image.Image, x: float, y: float) -> None:
1218 canvas.alpha_composite(icon, (int(round(x - icon.width / 2)), int(round(y - icon.height / 2))))
1219
1220
1221def draw_lines(canvas: Image.Image, proj: LensMap, sky: Sky, scale: float) -> None:
1222 layer = Image.new("RGBA", canvas.size)
1223 d = ImageDraw.Draw(layer)
1224 width = max(1, round(1.2 * scale))
1225 for a, b in sky.lines:
1226 for piece in great_circle(a, b, proj):
1227 d.line(piece, fill=(255, 255, 255, 110), width=width, joint="curve")
1228 canvas.alpha_composite(layer)
1229
1230
1231def draw_bodies(canvas, proj, sky: Sky, scale: float, font_path, sun_image) -> None:
1232 sx, sy = proj.forward(*sky.bodies["sun"])
1233 mx, my = proj.forward(*sky.bodies["moon"])
1234 font = ImageFont.truetype(str(font_path), max(8, round(34 * scale))) if font_path else None
1235 d = ImageDraw.Draw(canvas)
1236 for name, glyph in PLANETS.items():
1237 x, y = proj.forward(*sky.bodies[name])
1238 if font is not None:
1239 d.text((float(x), float(y)), glyph, font=font, anchor="mm", fill=(254, 128, 25, 255),
1240 stroke_width=max(1, round(2 * scale)), stroke_fill=(0, 0, 0, 200))
1241 else:
1242 r = 6 * scale
1243 d.ellipse((x - r, y - r, x + r, y + r), fill=(254, 128, 25, 255))
1244 dx, dy = float(sx - mx), float(my - sy)
1245 norm = math.hypot(dx, dy) or 1.0
1246 paste_centered(canvas, moon_icon(max(8, round(50 * scale)), sky.moon_elongation, (dx / norm, dy / norm)),
1247 float(mx), float(my))
1248 paste_centered(canvas, sun_icon(sun_image, max(8, round(100 * scale))), float(sx), float(sy))
1249
1250
1251def draw_markers(canvas, proj, markers: list[tuple[float, float]], icon_path, scale: float) -> None:
1252 if not markers or icon_path is None:
1253 return
1254 size = max(8, round(20 * scale))
1255 icon = Image.open(icon_path).convert("RGBA").resize((size, size), Image.LANCZOS)
1256 for lat, lon in markers:
1257 x, y = proj.forward(lat, lon)
1258 paste_centered(canvas, icon, float(x), float(y))
1259
1260
1261def load_markers(path: Path | None) -> list[tuple[float, float]]:
1262 if path is None:
1263 return []
1264 try:
1265 data = json.loads(path.read_text())
1266 except Exception as e: # noqa: BLE001
1267 log(f"warning: markers: {e}")
1268 return []
1269 entries = data.values() if isinstance(data, dict) else data
1270 out = []
1271 for v in entries:
1272 if isinstance(v, dict) and v.get("latitude") is not None and v.get("longitude") is not None:
1273 out.append((float(v["latitude"]), float(v["longitude"])))
1274 return out
1275
1276
1277def save(img: Image.Image, path: Path, **kw) -> None:
1278 tmp = path.with_name(f".{path.name}.tmp")
1279 img.save(tmp, format=path.suffix.lstrip(".").replace("jpg", "jpeg").upper(), **kw)
1280 os.replace(tmp, path)
1281
1282
1283# --------------------------------------------------------------------------
1284
1285
1286def main() -> None:
1287 env = os.environ.get
1288 p = argparse.ArgumentParser(prog="realwallpaper", description=__doc__.splitlines()[0])
1289 p.add_argument("--dir", type=Path, default=Path("."),
1290 help="output and cache directory (default: .)")
1291 p.add_argument("--lat", type=float, default=52.520008,
1292 help="home latitude in degrees, south negative (default Berlin)")
1293 p.add_argument("--lon", type=float, default=13.404954,
1294 help="home longitude in degrees, west negative (default Berlin)")
1295 p.add_argument("--width", type=int, default=3840)
1296 p.add_argument("--height", type=int, default=2160)
1297 p.add_argument("--magnification", type=float, default=6.0,
1298 help="zoom at home relative to the plain map (default 6)")
1299 p.add_argument("--lens", type=float, default=1.0,
1300 help="lens size relative to the screen; 1 reaches the edges (default 1)")
1301 p.add_argument("--lens-shape", choices=["circle", "ellipse", "rect"], default="ellipse",
1302 help="circle; ellipse or rect follow the screen's aspect (default ellipse)")
1303 p.add_argument("--focus", type=float, default=0.0,
1304 help="share of the lens zoomed undistorted before the transition (default 0)")
1305 p.add_argument("--falloff", type=float, default=0.65,
1306 help="where across the lens (0..1) the distortion peaks; higher keeps "
1307 "the middle closer to a plain zoom (default 0.65)")
1308 p.add_argument("--projection", choices=["lens", "plain", "mercator"], default="lens",
1309 help="lens: centred on home with the zoom lens; plain: centred on home "
1310 "without distortion; mercator: classic north-up map (default lens)")
1311 p.add_argument("--marker-url", default=None,
1312 help="JSON with {name: {latitude, longitude}} for krebs markers")
1313 p.add_argument("--time", type=dt.datetime.fromisoformat, default=None,
1314 help="render for this UTC time instead of now (ISO 8601)")
1315 p.add_argument("--no-archive", action="store_true", help="do not keep a jpg in archive/")
1316 args = p.parse_args()
1317
1318 when = args.time or dt.datetime.now(dt.timezone.utc)
1319 if when.tzinfo is None:
1320 when = when.replace(tzinfo=dt.timezone.utc)
1321 when = when.astimezone(dt.timezone.utc)
1322
1323 out_dir = args.dir.resolve()
1324 out_dir.mkdir(parents=True, exist_ok=True)
1325 cache = Cache(out_dir / "cache")
1326 proj = LensMap(args.lat, args.lon, args.width, args.height,
1327 args.magnification, args.lens, args.lens_shape, args.focus, args.falloff,
1328 args.projection)
1329 home_x, home_y = proj.forward(args.lat, args.lon)
1330 home_q = proj.inverse(np.array([home_x]), np.array([home_y]))[2][0]
1331 log(f"projection: {args.projection}, home at {home_q * 1000:.0f} m/px, "
1332 f"finest {proj.q0 * 1000:.0f} m/px")
1333
1334 font = env("REALWALLPAPER_FONT")
1335 constellations = env("REALWALLPAPER_CONSTELLATIONS")
1336 krebs_icon = env("REALWALLPAPER_KREBS_ICON")
1337
1338 t0 = time.monotonic()
1339 sun_image = cache.get_url("sun.jpg", "https://sdo.gsfc.nasa.gov/assets/img/latest/latest_512_0171.jpg", 3600)
1340 marker_file = cache.get_url("markers.json", args.marker_url, 0) if args.marker_url else None
1341 lay = load_layers(proj, cache, when)
1342 log(f"data ready in {time.monotonic() - t0:.1f}s")
1343
1344 try:
1345 sky = compute_sky(when, Path(constellations) if constellations else None)
1346 except Exception as e: # noqa: BLE001 - the map is still useful without the sky
1347 log(f"warning: sky: {e}")
1348 sky = None
1349 sun = sky.bodies["sun"] if sky else (0.0, 0.0)
1350
1351 t0 = time.monotonic()
1352 base = Image.fromarray(render_base(proj, lay, sun)).convert("RGBA")
1353 log(f"rendered in {time.monotonic() - t0:.1f}s")
1354
1355 scale = args.height / 1800
1356 markers = load_markers(marker_file)
1357
1358 plain = base
1359 marker = plain.copy()
1360 if sky:
1361 draw_bodies(marker, proj, sky, scale, font, sun_image)
1362 krebs = marker.copy()
1363 draw_markers(krebs, proj, markers, krebs_icon, scale)
1364 stars = plain.copy()
1365 if sky:
1366 draw_lines(stars, proj, sky, scale)
1367 draw_bodies(stars, proj, sky, scale, font, sun_image)
1368 draw_markers(stars, proj, markers, krebs_icon, scale)
1369
1370 for name, img in [
1371 ("realwallpaper.png", plain),
1372 ("realwallpaper-marker.png", marker),
1373 ("realwallpaper-krebs.png", krebs),
1374 ("realwallpaper-krebs-stars.png", stars),
1375 ]:
1376 save(img.convert("RGB"), out_dir / name, compress_level=3)
1377 if not args.no_archive:
1378 (out_dir / "archive").mkdir(exist_ok=True)
1379 save(krebs.convert("RGB"), out_dir / "archive" / f"{when.strftime('%Y-%m-%dT%H:%M:%SZ')}.jpg",
1380 quality=88)
1381 cache.cleanup()
1382 log("done")
1383
1384
1385if __name__ == "__main__":
1386 main()