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