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