forked from
lass.pds.lassul.us/superconfig
This repository has no description
tools.realwallpaper: real colour imagery on the day side
The day side now shows what the geostationary satellites see, like
zoom.earth, instead of the base map under white clouds synthesised from
infrared:
- Europe/Africa/Atlantic: MTG (Meteosat-12) rgb_truecolour from EUMETSAT.
rgb_geocolour looks the same by day but renders ~4x slower on the
server and timed out (502) for full 2048 px tiles.
- Americas/Pacific: GOES-East and GOES-West GeoColor from GIBS.
- Asia/Australia: GIBS has no Himawari colour product. Himawari band 3
(red, 10 min) is turned into white cloud over the base map, with the
cloud opacity being the reflectance beyond the ground's own albedo.
The daily VIIRS true colour mosaics were the alternative, but their
clouds are up to a day old, have swath seams of different times and
would not match the IR clouds next to them.
Each source is weighted by how well its satellite sees the point: full
trust up to 62° view zenith, none from 78° (the limb is smeared), and
the more nadir satellite wins in overlaps with cos^8 of the zenith
angle, so seams are soft. Trust also needs the sun not to be in front
of the satellite: looking towards a low sun both MTG (mornings over the
Indian Ocean) and Himawari (evenings) are washed out to a pale band
with straight edges, while a long path with the sun behind the
satellite is fine (Berlin at winter noon). A plain path-length limit
would have dropped the latter. Near the terminator (cosz < 0.2)
GeoColor fades into its purple night blend, so the live imagery fades
out there and the IR clouds over the base map take over as before.
Where it is shown, the synthetic clouds and so their relief shading and
shadows are scaled down, so the real image is not shaded twice.
All live layers (including the IR clouds) now request an explicit TIME,
the render time floored to 10 min: without TIME GIBS served tiles of
one layer from different frames (sun glint cut at a tile edge), and it
makes --time render the imagery of that time, not only the sun. The
TIME is capped at the newest frame the server lists (GIBS
DescribeDomains, EUMETSAT per-layer capabilities), since GIBS answers
an unpublished time with an empty image and EUMETSAT with an error;
MTG true colour stops producing frames while its disk is dark.
GIBS also now and then answers with an empty 1 kB image for a frame it
has (the same request succeeds when repeated), which left a quarter of
the Himawari disk missing, so empty tiles are asked for up to 3 times.
author
Hermes
committer
lassulus
date
(Oct 17, 2026, 1:32 AM CEST)
commit
02965f48
02965f489aea0742d049c65ba09c25db33bdc180
parent
f044be6c
f044be6cf9befe576fda3ab1b6a722a7b11d17ea
+141
-9
+141
-9
tools/realwallpaper/realwallpaper.py
+141
-9
tools/realwallpaper/realwallpaper.py
···
7
7
8
8
Imagery is fetched over WMS at the resolution each pixel needs: a few local
9
9
levels around home, one global level, and downsampled global mipmaps. The
10
-
night side is the day map darkened, plus city lights as a glow. Every fetch
10
+
day side shows live colour imagery of the geostationary satellites where
11
+
they see it well, elsewhere the base map under infrared clouds. The night
12
+
side is the day map darkened, plus city lights as a glow. Every fetch
11
13
is cached with a TTL; a failing source falls back to its last good copy or
12
14
is skipped, so one dead upstream no longer kills the whole wallpaper.
13
15
"""
···
38
40
39
41
R_EARTH = 6371.0088 # km
40
42
KM_PER_DEG = math.pi * R_EARTH / 180
43
+
GEO_RADIUS = 42164.0 / R_EARTH # geostationary orbit in earth radii
41
44
GLOBAL_WIDTH = 8192 # px of the global equirectangular level
42
45
# px per WMS request side. The servers render on the fly at a roughly fixed
43
46
# px/s per request (EOX: ~0.35 Mpx/s) but serve requests in parallel, so
···
312
315
oversample: int = 1
313
316
# per-channel gain to match the colours of the other day source
314
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
315
321
316
322
def times(self, when: dt.datetime) -> list[str | None]:
317
323
if self.days_back:
···
486
492
return np.concatenate([rgb.astype(np.uint8), rgba[..., 3:4]], axis=-1)
487
493
488
494
495
+
def prep_grey(layer: "Layer", rgba: np.ndarray, cmap: Colormap | None) -> np.ndarray:
496
+
return prep_rgb(layer, rgba, cmap)[..., [0, 3]]
497
+
498
+
489
499
def prep_lights(layer: "Layer", rgba: np.ndarray, cmap: Colormap | None) -> np.ndarray:
490
500
"""Black Marble -> city light intensity; drops its dim blue land texture."""
491
501
red = rgba[..., 0].astype(np.float32) / 255
···
608
618
w, e = max(-180.0, w + shift), min(180.0, e + shift)
609
619
n = level.north - r0 * px_lat
610
620
s = level.north - r1 * px_lat
611
-
tile = fetch_image(src.url((s, w, n, e), c1 - c0, r1 - r0, t))
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
612
627
if tile.shape[:2] != (r1 - r0, c1 - c0):
613
628
tile = np.asarray(Image.fromarray(tile).resize((c1 - c0, r1 - r0), Image.BILINEAR))
614
629
img[r0:r1, c0:c1] = tile
···
630
645
combine: Callable[[list[np.ndarray]], np.ndarray] = combine_max
631
646
colormap: str | None = None # GIBS colormap name
632
647
transform: Callable[[np.ndarray], np.ndarray] | None = None
648
+
sat_lon: float | None = None # sub-satellite longitude of geostationary day imagery
633
649
required: bool = False
634
650
cmap: Colormap | None = None
635
651
levels: list[Level] = field(default_factory=list)
···
750
766
return smoothstep(0.42, 0.78, ir) ** 1.2
751
767
752
768
753
-
def layers() -> list[Layer]:
769
+
def layers(when: dt.datetime) -> list[Layer]:
754
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)
755
775
return [
756
776
Layer(
757
777
"day",
···
771
791
),
772
792
Layer(
773
793
"clouds",
774
-
Wms("mtg-ir105", EUMETSAT, "mtg_fd:ir105_hrfi", 1.0, 600),
775
-
Wms("georing-ir108", EUMETSAT, "mumi:worldcloudmap_ir108", 3.0, 3600),
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),
776
796
prep_cloud,
777
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,
778
826
),
779
827
Layer(
780
828
"snow",
···
806
854
]
807
855
808
856
857
+
def 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
+
809
886
def load_layers(proj: LensMap, cache: Cache, when: dt.datetime) -> dict[str, Layer]:
810
887
all_layers = []
811
-
for layer in 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))
812
894
if layer.colormap:
813
895
path = cache.get_url(f"colormap-{layer.colormap}.xml",
814
896
f"https://gibs.earthdata.nasa.gov/colormaps/v1.3/{layer.colormap}.xml",
···
974
1056
return np.sin(la) * math.sin(sl) + np.cos(la) * math.cos(sl) * np.cos(lo - so)
975
1057
976
1058
1059
+
def 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
+
977
1104
# --------------------------------------------------------------------------
978
1105
# cloud post-processing: the best live cloud data is ~1 km, far coarser than
979
1106
# the zoomed-in pixels around home. Synthesise the missing detail (noise
···
1112
1239
1113
1240
cosz = subsolar_cos(lat, lon, *sun)
1114
1241
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)
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)
1117
1249
if (s := get("clouds")) is not None:
1118
-
cloud[rows] = np.clip(s[0][:, 0], 0, 1).reshape(len(ys), W)
1250
+
cloud[rows] = (np.clip(s[0][:, 0], 0, 1) * (1 - sunny)).reshape(len(ys), W)
1119
1251
if (s := get("lights")) is not None:
1120
1252
lights[rows] = (s[0][:, 0] * (1 - light)).reshape(len(ys), W)
1121
1253