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