The Moon's shadow from Meteosat — the 12 Aug 2026 total solar eclipse (MTG-FCI)¶
On 12 August 2026 the Moon's shadow swept the far north — the Arctic, eastern Greenland, western Iceland, the North Atlantic and northern Spain — with greatest eclipse at 17:47 UTC. Unlike GOES-East (which sees this reach only at its limb), Meteosat-12 (MTG-I1) sits at 0°E, so its FCI full-disk imager looks straight at the eclipse.
We build true-colour imagery from FCI Level-1c granules on the EUMETSAT Data Store and render the shadow two ways, both from-space perspectives on the same full disk: a vertical near-side view (looking straight down on the North Atlantic) and a tilted, curved-limb view (the low-angle shot a crew in orbit would see, shadow riding the northern limb).
About Meteosat-12 / MTG-I1¶
Meteosat is EUMETSAT's family of European geostationary weather satellites. MTG-I1 (Meteosat Third Generation – Imaging 1, operationally Meteosat-12) carries the Flexible Combined Imager (FCI), a 16-band imager that scans the full Earth disk every 10 minutes. Like GOES it is geostationary — ~35,786 km over the equator, turning with Earth so it stares at one fixed hemisphere.
The crucial difference for this eclipse is where it is parked: MTG-I1 sits at 0° longitude (the prime meridian), directly over the Gulf of Guinea. That places Europe and the North Atlantic — where the 12 Aug 2026 eclipse tracked — near the centre of its disk, so it looks almost head-on at the shadow instead of catching it at the limb the way GOES-East does. That is why this notebook's imagery shows the eclipse so much more dramatically.
In [ ]:
Copied!
import gc
import hashlib
import shutil
import warnings
import zipfile
from pathlib import Path
import cleopatra.glyphs.gridded.array_glyph as _array_glyph
import matplotlib.patheffects as pe
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from cleopatra.glyphs.gridded.array_glyph import Animation, FrameLabel
from dotenv import find_dotenv, load_dotenv
from IPython.display import Image, display
from loguru import logger
from matplotlib.patches import Circle, FancyBboxPatch, Wedge
from pyproj import CRS, Transformer
from pyramids.dataset import Dataset, GeoReference
from pyramids.dataset.collection import DatasetCollection
from pyramids_eo.composites import (
night_ir,
rayleigh_correct,
relative_azimuth,
satellite_zenith_azimuth,
solar_zenith_azimuth,
static_image,
sunz_correct,
sunz_reduce,
true_color,
true_color_with_night_ir,
)
from pyramids_eo.enhance import stretch
from pyramids_eo.resample import to_area
from pyramids_eo.sensors.readers import read_fci_l1c
from scipy.ndimage import binary_erosion, gaussian_filter
from earthlens.core import EarthLens, cache_dir
warnings.filterwarnings("ignore")
logger.remove()
plt.rcParams["figure.dpi"] = 80
import gc
import hashlib
import shutil
import warnings
import zipfile
from pathlib import Path
import cleopatra.glyphs.gridded.array_glyph as _array_glyph
import matplotlib.patheffects as pe
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from cleopatra.glyphs.gridded.array_glyph import Animation, FrameLabel
from dotenv import find_dotenv, load_dotenv
from IPython.display import Image, display
from loguru import logger
from matplotlib.patches import Circle, FancyBboxPatch, Wedge
from pyproj import CRS, Transformer
from pyramids.dataset import Dataset, GeoReference
from pyramids.dataset.collection import DatasetCollection
from pyramids_eo.composites import (
night_ir,
rayleigh_correct,
relative_azimuth,
satellite_zenith_azimuth,
solar_zenith_azimuth,
static_image,
sunz_correct,
sunz_reduce,
true_color,
true_color_with_night_ir,
)
from pyramids_eo.enhance import stretch
from pyramids_eo.resample import to_area
from pyramids_eo.sensors.readers import read_fci_l1c
from scipy.ndimage import binary_erosion, gaussian_filter
from earthlens.core import EarthLens, cache_dir
warnings.filterwarnings("ignore")
logger.remove()
plt.rcParams["figure.dpi"] = 80
In [ ]:
Copied!
# A quick schematic of the geometry described above (illustrative; the orbit ring is to scale
# relative to Earth's radius, the satellite marker is not).
R, D = (
1.0,
6.61,
) # Earth radius; geostationary orbit radius = 42,164 km / 6,371 km Earth radii
half = np.degrees(
np.arccos(R / D)
) # half-angle of the visible cap (~81 deg -> nearly a full hemisphere)
fig, ax = plt.subplots(figsize=(9, 6.2))
ax.set_aspect("equal")
ax.axis("off")
ax.add_patch(
Circle((0, 0), R, facecolor="#1b3a63", edgecolor="#0d2440", lw=1.5, zorder=3)
)
ax.add_patch(Wedge((0, 0), R, -half, half, facecolor="#2e6fb0", alpha=0.9, zorder=4))
ax.add_patch(Circle((0, 0), R, facecolor="none", edgecolor="#0d2440", lw=1.5, zorder=6))
ax.plot([-R, R], [0, 0], color="#9fb8d6", lw=0.8, ls=(0, (4, 3)), zorder=5) # equator
ax.text(
0,
R + 0.12,
"Earth",
ha="center",
va="bottom",
fontsize=11,
color="#1b3a63",
weight="bold",
)
theta = np.linspace(0, 2 * np.pi, 400)
ax.plot(
D * np.cos(theta),
D * np.sin(theta),
color="#8a8a8a",
lw=1.0,
ls=(0, (6, 4)),
zorder=1,
)
ax.text(
0,
-D - 0.35,
"geostationary orbit",
ha="center",
va="top",
fontsize=9,
color="#6a6a6a",
style="italic",
)
sx = D
ax.plot([sx], [0], marker="s", ms=11, color="#7a3ea8", zorder=7)
ax.text(
sx + 0.15,
0,
"Meteosat-12\n(MTG-I1)",
ha="left",
va="center",
fontsize=10,
weight="bold",
color="#5c2e82",
)
tx, ty = R * np.cos(np.radians(half)), R * np.sin(np.radians(half))
ax.plot([sx, tx], [0, ty], color="#7a3ea8", lw=1.1, alpha=0.75, zorder=2)
ax.plot([sx, tx], [0, -ty], color="#7a3ea8", lw=1.1, alpha=0.75, zorder=2)
ax.fill([sx, tx, tx], [0, ty, -ty], color="#7a3ea8", alpha=0.06, zorder=1)
ax.plot([R], [0], marker="o", ms=7, color="#f2c200", mec="#7a6400", zorder=8)
ax.annotate(
"sub-satellite point\n0\u00b0N, 0\u00b0E\n(Gulf of Guinea)",
xy=(R, 0),
xytext=(-1.9, 2.6),
ha="left",
va="bottom",
fontsize=9,
color="#5a4a00",
arrowprops={"arrowstyle": "->", "color": "#7a6400", "lw": 1},
)
ax.annotate(
"",
xy=(R, -0.5),
xytext=(D, -0.5),
arrowprops={"arrowstyle": "<->", "color": "#444", "lw": 1.1},
)
ax.text(
(R + D) / 2,
-0.72,
"35,786 km altitude",
ha="center",
va="top",
fontsize=9,
color="#333",
)
ax.text(
sx - 0.3,
ty + 1.15,
"sees one full hemisphere\n(Europe, Africa, the Atlantic)",
ha="center",
va="bottom",
fontsize=9,
color="#5c2e82",
style="italic",
)
ax.text(
0.0,
-R - 1.7,
"The eclipse (far-north Atlantic / Europe) lies near this central meridian,\n"
"so Meteosat sees it close to disk-centre and head-on \u2014 where GOES-East\n"
"saw the same eclipse out at its limb.",
ha="center",
va="top",
fontsize=8.5,
color="#444",
)
ax.set_xlim(-2.0, D + 2.6)
ax.set_ylim(-D - 1.8, D + 0.6)
ax.set_title(
"Geostationary geometry \u2014 Meteosat-12 (MTG-I1) hangs fixed at 0\u00b0 and looks head-on at Europe",
fontsize=11.5,
color="#222",
pad=10,
)
plt.show()
# A quick schematic of the geometry described above (illustrative; the orbit ring is to scale
# relative to Earth's radius, the satellite marker is not).
R, D = (
1.0,
6.61,
) # Earth radius; geostationary orbit radius = 42,164 km / 6,371 km Earth radii
half = np.degrees(
np.arccos(R / D)
) # half-angle of the visible cap (~81 deg -> nearly a full hemisphere)
fig, ax = plt.subplots(figsize=(9, 6.2))
ax.set_aspect("equal")
ax.axis("off")
ax.add_patch(
Circle((0, 0), R, facecolor="#1b3a63", edgecolor="#0d2440", lw=1.5, zorder=3)
)
ax.add_patch(Wedge((0, 0), R, -half, half, facecolor="#2e6fb0", alpha=0.9, zorder=4))
ax.add_patch(Circle((0, 0), R, facecolor="none", edgecolor="#0d2440", lw=1.5, zorder=6))
ax.plot([-R, R], [0, 0], color="#9fb8d6", lw=0.8, ls=(0, (4, 3)), zorder=5) # equator
ax.text(
0,
R + 0.12,
"Earth",
ha="center",
va="bottom",
fontsize=11,
color="#1b3a63",
weight="bold",
)
theta = np.linspace(0, 2 * np.pi, 400)
ax.plot(
D * np.cos(theta),
D * np.sin(theta),
color="#8a8a8a",
lw=1.0,
ls=(0, (6, 4)),
zorder=1,
)
ax.text(
0,
-D - 0.35,
"geostationary orbit",
ha="center",
va="top",
fontsize=9,
color="#6a6a6a",
style="italic",
)
sx = D
ax.plot([sx], [0], marker="s", ms=11, color="#7a3ea8", zorder=7)
ax.text(
sx + 0.15,
0,
"Meteosat-12\n(MTG-I1)",
ha="left",
va="center",
fontsize=10,
weight="bold",
color="#5c2e82",
)
tx, ty = R * np.cos(np.radians(half)), R * np.sin(np.radians(half))
ax.plot([sx, tx], [0, ty], color="#7a3ea8", lw=1.1, alpha=0.75, zorder=2)
ax.plot([sx, tx], [0, -ty], color="#7a3ea8", lw=1.1, alpha=0.75, zorder=2)
ax.fill([sx, tx, tx], [0, ty, -ty], color="#7a3ea8", alpha=0.06, zorder=1)
ax.plot([R], [0], marker="o", ms=7, color="#f2c200", mec="#7a6400", zorder=8)
ax.annotate(
"sub-satellite point\n0\u00b0N, 0\u00b0E\n(Gulf of Guinea)",
xy=(R, 0),
xytext=(-1.9, 2.6),
ha="left",
va="bottom",
fontsize=9,
color="#5a4a00",
arrowprops={"arrowstyle": "->", "color": "#7a6400", "lw": 1},
)
ax.annotate(
"",
xy=(R, -0.5),
xytext=(D, -0.5),
arrowprops={"arrowstyle": "<->", "color": "#444", "lw": 1.1},
)
ax.text(
(R + D) / 2,
-0.72,
"35,786 km altitude",
ha="center",
va="top",
fontsize=9,
color="#333",
)
ax.text(
sx - 0.3,
ty + 1.15,
"sees one full hemisphere\n(Europe, Africa, the Atlantic)",
ha="center",
va="bottom",
fontsize=9,
color="#5c2e82",
style="italic",
)
ax.text(
0.0,
-R - 1.7,
"The eclipse (far-north Atlantic / Europe) lies near this central meridian,\n"
"so Meteosat sees it close to disk-centre and head-on \u2014 where GOES-East\n"
"saw the same eclipse out at its limb.",
ha="center",
va="top",
fontsize=8.5,
color="#444",
)
ax.set_xlim(-2.0, D + 2.6)
ax.set_ylim(-D - 1.8, D + 0.6)
ax.set_title(
"Geostationary geometry \u2014 Meteosat-12 (MTG-I1) hangs fixed at 0\u00b0 and looks head-on at Europe",
fontsize=11.5,
color="#222",
pad=10,
)
plt.show()
How the frames are built¶
earthlens reaches EUMETSAT through the eumdac Data Store client (the same client the
earthlens.eumetsat backend wraps). FCI L1C has a 10-minute repeat cycle, so we query the Data Store
for each 10-minute slot from 16:30 to 19:00 UTC and download that one full-disk cycle (~770 MB, ~40
chunk files). pyramids-eo's read_fci_l1c calibrates and assembles the channels, and its compositing
primitives build the day/night frame: true_color by day (Rayleigh-corrected in-package), grey thermal
cloud tops over NASA Black Marble city lights by night. We resample once per projection, then clean up the
frame edges. The raw granules stay in the cache dir and the resampled frames are cached too, so a re-run
costs no downloads and no resampling.
In [ ]:
Copied!
SLOTS = pd.date_range("2026-08-12 16:30", "2026-08-12 19:00", freq="10min") # 16 cycles
DATASET = "mtg-fci-l1c" # earthlens catalog key for EO:EUM:DAT:0662, 10-min full disk
DISC = {"lat_lim": [-79.0, 79.0], "lon_lim": [-79.0, 79.0]} # the FCI full disk
OUT = Path("out") / "eclipse_fci"
# Raw granules are staged under the earthlens cache dir (set_cache_dir() /
# EARTHLENS_CACHE, else the per-platform user cache) and KEPT, so re-running -- or
# re-framing to another projection later -- never re-downloads the ~770 MB/cycle.
RAW_DIR = cache_dir() / "eclipse_fci" / "raw"
RAW_DIR.mkdir(parents=True, exist_ok=True)
# Resampling is the expensive step, so its output is cached in `frames_*_raw` and the
# edge treatment below writes its own copies into `frames_*`. Keeping the two apart
# means the cheap pass can be retuned and re-run without resampling anything again.
NSPER_RAW = OUT / "frames_nsper_raw"
TPERS_RAW = OUT / "frames_tpers_raw"
NSPER_DIR = OUT / "frames_nsper"
TPERS_DIR = OUT / "frames_tpers"
for d in (NSPER_RAW, TPERS_RAW, NSPER_DIR, TPERS_DIR):
d.mkdir(parents=True, exist_ok=True)
# Two "from space" perspectives on the eclipse track, both rendered from the full disk:
# 1) NSPER -- vertical near-side perspective, looking straight down on the N Atlantic.
# 2) TPERS -- *tilted* perspective: the low-angle, curved-limb view a crew in orbit would
# see, with Iberia in the foreground and the shadow riding the northern limb.
NSPER = {
"crs": (
"+proj=nsper +lat_0=60 +lon_0=-10 +h=35785831 "
"+a=6378137 +b=6356752.314 +units=m +no_defs"
),
"extent": (-4.2e6, -4.2e6, 4.2e6, 4.2e6),
"width": 1600,
"height": 1600,
}
TPERS = {
"crs": (
"+proj=tpers +lat_0=45 +lon_0=0 +tilt=28 +azi=0 +h=35785831 "
"+a=6378137 +b=6356752.314 +units=m +no_defs"
),
"extent": (-5.0e6, -3.2e6, 5.0e6, 3.4e6),
"width": 1680,
"height": 1120,
}
# The five channels the composite needs: FCI's blue / green / red / veggie bands for true
# colour, and the 10.5 um window channel for the night side.
SOLAR_CHANNELS = ["vis_04", "vis_05", "vis_06", "vis_08"]
IR_CHANNEL = "ir_105"
# The threshold between the solar bands' native 1 km grid and the IR channel's native
# 2 km one, used below to tell which pixel size a just-read channel came in at.
NATIVE_1KM_THRESHOLD_M = 1500
# EUMETSAT credentials come from the repo-root .env (EUMETSAT_CONSUMER_KEY / _SECRET).
# The earthlens backend mints its own token from those same variables, so nothing else
# is needed here.
load_dotenv(find_dotenv(usecwd=True))
# --- Day/night city-lights ------------------------------------------------------
# The night side is painted with NASA Black Marble (VIIRS city lights). `static_image`
# downloads it once, caches it, and warps it onto whatever grid it is given, so there is
# no composite configuration to write and no dead bundled URL to work around.
# Black Marble and the cached lon/lat grid are large, derived caches rather than
# notebook outputs, so they live beside the raw granules under the cache dir.
AUX = cache_dir() / "eclipse_fci" / "aux"
AUX.mkdir(parents=True, exist_ok=True)
BLACK_MARBLE = (
"https://eoimages.gsfc.nasa.gov/images/imagerecords/144000/144898/"
"BlackMarble_2016_3km_geo.tif"
)
# Band centre wavelengths for the Rayleigh correction. Molecular scattering goes as
# wavelength^-4, so the blue band needs several times the correction the red band does;
# without it the disc carries a blue veil, measured at +37 DN on the blue channel.
RAYLEIGH_UM = {"vis_04": 0.444, "vis_05": 0.510, "vis_06": 0.640}
SLOTS = pd.date_range("2026-08-12 16:30", "2026-08-12 19:00", freq="10min") # 16 cycles
DATASET = "mtg-fci-l1c" # earthlens catalog key for EO:EUM:DAT:0662, 10-min full disk
DISC = {"lat_lim": [-79.0, 79.0], "lon_lim": [-79.0, 79.0]} # the FCI full disk
OUT = Path("out") / "eclipse_fci"
# Raw granules are staged under the earthlens cache dir (set_cache_dir() /
# EARTHLENS_CACHE, else the per-platform user cache) and KEPT, so re-running -- or
# re-framing to another projection later -- never re-downloads the ~770 MB/cycle.
RAW_DIR = cache_dir() / "eclipse_fci" / "raw"
RAW_DIR.mkdir(parents=True, exist_ok=True)
# Resampling is the expensive step, so its output is cached in `frames_*_raw` and the
# edge treatment below writes its own copies into `frames_*`. Keeping the two apart
# means the cheap pass can be retuned and re-run without resampling anything again.
NSPER_RAW = OUT / "frames_nsper_raw"
TPERS_RAW = OUT / "frames_tpers_raw"
NSPER_DIR = OUT / "frames_nsper"
TPERS_DIR = OUT / "frames_tpers"
for d in (NSPER_RAW, TPERS_RAW, NSPER_DIR, TPERS_DIR):
d.mkdir(parents=True, exist_ok=True)
# Two "from space" perspectives on the eclipse track, both rendered from the full disk:
# 1) NSPER -- vertical near-side perspective, looking straight down on the N Atlantic.
# 2) TPERS -- *tilted* perspective: the low-angle, curved-limb view a crew in orbit would
# see, with Iberia in the foreground and the shadow riding the northern limb.
NSPER = {
"crs": (
"+proj=nsper +lat_0=60 +lon_0=-10 +h=35785831 "
"+a=6378137 +b=6356752.314 +units=m +no_defs"
),
"extent": (-4.2e6, -4.2e6, 4.2e6, 4.2e6),
"width": 1600,
"height": 1600,
}
TPERS = {
"crs": (
"+proj=tpers +lat_0=45 +lon_0=0 +tilt=28 +azi=0 +h=35785831 "
"+a=6378137 +b=6356752.314 +units=m +no_defs"
),
"extent": (-5.0e6, -3.2e6, 5.0e6, 3.4e6),
"width": 1680,
"height": 1120,
}
# The five channels the composite needs: FCI's blue / green / red / veggie bands for true
# colour, and the 10.5 um window channel for the night side.
SOLAR_CHANNELS = ["vis_04", "vis_05", "vis_06", "vis_08"]
IR_CHANNEL = "ir_105"
# The threshold between the solar bands' native 1 km grid and the IR channel's native
# 2 km one, used below to tell which pixel size a just-read channel came in at.
NATIVE_1KM_THRESHOLD_M = 1500
# EUMETSAT credentials come from the repo-root .env (EUMETSAT_CONSUMER_KEY / _SECRET).
# The earthlens backend mints its own token from those same variables, so nothing else
# is needed here.
load_dotenv(find_dotenv(usecwd=True))
# --- Day/night city-lights ------------------------------------------------------
# The night side is painted with NASA Black Marble (VIIRS city lights). `static_image`
# downloads it once, caches it, and warps it onto whatever grid it is given, so there is
# no composite configuration to write and no dead bundled URL to work around.
# Black Marble and the cached lon/lat grid are large, derived caches rather than
# notebook outputs, so they live beside the raw granules under the cache dir.
AUX = cache_dir() / "eclipse_fci" / "aux"
AUX.mkdir(parents=True, exist_ok=True)
BLACK_MARBLE = (
"https://eoimages.gsfc.nasa.gov/images/imagerecords/144000/144898/"
"BlackMarble_2016_3km_geo.tif"
)
# Band centre wavelengths for the Rayleigh correction. Molecular scattering goes as
# wavelength^-4, so the blue band needs several times the correction the red band does;
# without it the disc carries a blue veil, measured at +37 DN on the blue channel.
RAYLEIGH_UM = {"vis_04": 0.444, "vis_05": 0.510, "vis_06": 0.640}
Geometry for the atmospheric correction¶
Rayleigh correction needs to know where the Sun and the satellite are for every pixel. The solar zenith comes from pyramids-eo; the satellite angles and the solar azimuth are plain geostationary geometry, computed once for the fixed disc grid.
In [ ]:
Copied!
# --- The disc grid ----------------------------------------------------------------
# pyramids-eo supplies the solar and viewing geometry itself, so all this notebook
# still needs is the longitude / latitude of the fixed disc grid to feed it.
def disc_lonlat(geo, crs, shape):
"""Longitude / latitude of every cell centre on the geostationary grid.
The grid never moves, so this is computed once and reused for every cycle.
Args:
geo: The dataset geotransform.
crs: The dataset CRS, as WKT.
shape: The `(rows, columns)` of the grid.
Returns:
A `(lon, lat)` pair of arrays; off-disc cells are not finite.
"""
# The cache key folds in the grid's shape and geotransform (and a hash of the CRS
# WKT) so a changed grid -- a different channel, a resampled resolution, an
# upstream geotransform fix -- gets a fresh cache entry instead of silently
# reusing coordinates for a grid it was never computed on.
rows, cols = shape
# sha256, not sha1: a cache-key derivation, not a security use of hashing, but
# sha1 still trips SonarCloud's weak-hash rule (S4790).
key = hashlib.sha256(
f"{rows}x{cols}:{geo}:{hashlib.sha256(crs.encode()).hexdigest()}".encode()
).hexdigest()[:16]
cache = AUX / f"disc_lonlat_{key}.npz"
if cache.exists():
# The hash key already changes with the grid, so this is not the guard
# against a changed grid -- that never reaches here under a different key.
# It only catches a hash collision or a corrupted/manually-placed file at
# this exact key, which is why it falls through to recomputing rather than
# trusting a mismatched cache blindly.
cached = np.load(cache)
if cached["lon"].shape == shape:
return cached["lon"], cached["lat"]
xs = geo[0] + (np.arange(cols) + 0.5) * geo[1]
ys = geo[3] + (np.arange(rows) + 0.5) * geo[5]
grid_x, grid_y = np.meshgrid(xs, ys)
transformer = Transformer.from_crs(
CRS.from_wkt(crs), CRS.from_epsg(4326), always_xy=True
)
lon, lat = transformer.transform(grid_x, grid_y)
lon, lat = lon.astype("float32"), lat.astype("float32")
np.savez_compressed(cache, lon=lon, lat=lat)
return lon, lat
# --- The disc grid ----------------------------------------------------------------
# pyramids-eo supplies the solar and viewing geometry itself, so all this notebook
# still needs is the longitude / latitude of the fixed disc grid to feed it.
def disc_lonlat(geo, crs, shape):
"""Longitude / latitude of every cell centre on the geostationary grid.
The grid never moves, so this is computed once and reused for every cycle.
Args:
geo: The dataset geotransform.
crs: The dataset CRS, as WKT.
shape: The `(rows, columns)` of the grid.
Returns:
A `(lon, lat)` pair of arrays; off-disc cells are not finite.
"""
# The cache key folds in the grid's shape and geotransform (and a hash of the CRS
# WKT) so a changed grid -- a different channel, a resampled resolution, an
# upstream geotransform fix -- gets a fresh cache entry instead of silently
# reusing coordinates for a grid it was never computed on.
rows, cols = shape
# sha256, not sha1: a cache-key derivation, not a security use of hashing, but
# sha1 still trips SonarCloud's weak-hash rule (S4790).
key = hashlib.sha256(
f"{rows}x{cols}:{geo}:{hashlib.sha256(crs.encode()).hexdigest()}".encode()
).hexdigest()[:16]
cache = AUX / f"disc_lonlat_{key}.npz"
if cache.exists():
# The hash key already changes with the grid, so this is not the guard
# against a changed grid -- that never reaches here under a different key.
# It only catches a hash collision or a corrupted/manually-placed file at
# this exact key, which is why it falls through to recomputing rather than
# trusting a mismatched cache blindly.
cached = np.load(cache)
if cached["lon"].shape == shape:
return cached["lon"], cached["lat"]
xs = geo[0] + (np.arange(cols) + 0.5) * geo[1]
ys = geo[3] + (np.arange(rows) + 0.5) * geo[5]
grid_x, grid_y = np.meshgrid(xs, ys)
transformer = Transformer.from_crs(
CRS.from_wkt(crs), CRS.from_epsg(4326), always_xy=True
)
lon, lat = transformer.transform(grid_x, grid_y)
lon, lat = lon.astype("float32"), lat.astype("float32")
np.savez_compressed(cache, lon=lon, lat=lat)
return lon, lat
Cleaning up the frame edges¶
Two artefacts survive the resample, and both sit exactly where the eye goes first.
The first is a teal ridge along the day/night seam. An IR false colour with cold cloud tops coming out
cyan would be the obvious suspect on the night side — but build_composite below already renders the night
clouds in neutral grey (night_ir(0.92 * grey, ...)), and the ridge survives that. It is on the day
side of the blend: at the terminator the true-colour Rayleigh correction is working at a grazing solar
angle and over- or under-corrects green and blue depending on how closely the correction matches the real
atmosphere at that angle. declip_cyan below removes it by comparing the green/blue excess against its own
locally-blurred background, so a thin bright ridge is pulled down while broad real colour is left alone.
The second is the hard cut where FCI's scan coverage stops: true_color_with_night_ir(keep_alpha=True)
derives its coverage band as a binary 0/1 mask (0 off-disk), so the imagery ends in a razor edge against
black. soften ramps the brightness down over the last few pixels and adds the atmospheric halo that real
full-disk imagery carries, because at the horizon the line of sight grazes the atmosphere. Neither pass
invents imagery — both only attenuate or tint pixels that already carry data.
In [ ]:
Copied!
def declip_cyan(rgb, tol=10.0, sigma=22.0, strength=0.9):
"""Tame a localised cyan/teal excess along the day/night seam, if present.
Written against satpy's Rayleigh correction, which over-corrected at grazing solar
angles and left a sharp ridge a few pixels wide. pyramids-eo's correction is weaker
(measured in the PR that ported this notebook off satpy: it removes roughly a third
less path radiance), which shows up as a broader, disc-wide blue veil rather than a
sharp localised ridge -- and measured directly on the rendered frames, this pass now
changes very little: the count of pixels with an anomalous green/blue excess is
essentially unchanged before and after. It is kept as a low-risk residual safeguard
against any localised artifact rather than as the fix for the disc-wide bias, which
would need retuning verified against a fresh render. Clamping colour globally would
also wash out real vegetation and shallow water, so the excess is compared against
its own locally-blurred background instead: a broad region matches its background
and is left alone, while a localised spike stands far above it and is pulled down.
Args:
rgb: `(H, W, 3)` uint8 image.
tol: How far above the local background the excess may sit before clamping.
sigma: Blur radius, in pixels, defining "local background".
strength: How much of the anomaly to remove, 0..1.
Returns:
`(H, W, 3)` uint8 image with the seam neutralised.
"""
out = rgb.astype(np.float32).copy()
red = out[..., 0]
for chan in (1, 2):
excess = out[..., chan] - red
ridge = excess - gaussian_filter(excess, sigma=sigma)
out[..., chan] -= strength * np.clip(ridge - tol, 0.0, None)
return np.clip(out, 0, 255).astype(np.uint8)
def soften(rgb, alpha, feather=14.0, glow=0.55, glow_rgb=(0.60, 0.78, 1.0), trim=4):
"""Soften the hard edge where FCI's scan coverage stops, and add a limb glow.
pyramids-eo's true_color_with_night_ir(keep_alpha=True) derives a binary 0/1 coverage
band, so the imagery ends in a razor cut against black. Two passes fix it: the feather ramps brightness down over the last
few pixels of real data, and the glow adds the atmospheric halo that real full-disk
imagery carries, because the line of sight grazes the atmosphere at the horizon.
Neither pass invents imagery -- both only attenuate or tint pixels that already carry
data.
Args:
rgb: `(H, W, 3)` uint8 image.
alpha: `(H, W)` binary coverage mask (0 outside the scan, 255 inside).
feather: Gaussian sigma, in pixels, of the fade-out ramp.
glow: Strength of the atmospheric halo; 0 disables it.
glow_rgb: Colour of the halo.
trim: Pixels of the outermost, ragged rim to drop before feathering.
Returns:
`(H, W, 3)` uint8 image with a soft boundary.
"""
data = alpha > 0
if trim:
data = binary_erosion(data, iterations=int(trim))
data = data.astype(np.float32)
ramp = np.clip(gaussian_filter(data, sigma=feather), 0.0, 1.0)
out = rgb.astype(np.float32) * ramp[..., None]
if glow > 0:
limb = (
np.clip(
gaussian_filter(data, sigma=feather * 0.55)
- gaussian_filter(data, sigma=feather * 1.9),
0.0,
1.0,
)
* data
)
peak = float(limb.max())
if peak > 1e-6:
limb = np.clip(limb / peak, 0.0, 1.0) ** 1.3 # tighten the rim
out += (
glow * 255.0 * limb[..., None] * np.asarray(glow_rgb, dtype=np.float32)
)
return np.clip(out, 0, 255).astype(np.uint8)
def soften_file(path, **kwargs):
"""Read an RGBA GeoTIFF frame, clean its seam and boundary, write it back in place."""
ds = Dataset.read_file(path, read_only=False)
bands = ds.read_array()
# read_array returns (bands, rows, cols), so `bands[:3]` slices the band
# axis; a single-band read comes back 2-D and would slice rows instead.
if bands.ndim != 3:
raise ValueError(f"expected an RGBA frame, got shape {bands.shape}")
rgb = np.dstack(bands[:3])
alpha = (
bands[3] if ds.band_count > 3 else (rgb.max(axis=2) > 0).astype(np.uint8) * 255
)
out = soften(declip_cyan(rgb), alpha, **kwargs)
for i in range(3):
ds.write_array(out[..., i], band=i)
ds.close() # flush the in-place writes before the next frame reads them
return path
def declip_cyan(rgb, tol=10.0, sigma=22.0, strength=0.9):
"""Tame a localised cyan/teal excess along the day/night seam, if present.
Written against satpy's Rayleigh correction, which over-corrected at grazing solar
angles and left a sharp ridge a few pixels wide. pyramids-eo's correction is weaker
(measured in the PR that ported this notebook off satpy: it removes roughly a third
less path radiance), which shows up as a broader, disc-wide blue veil rather than a
sharp localised ridge -- and measured directly on the rendered frames, this pass now
changes very little: the count of pixels with an anomalous green/blue excess is
essentially unchanged before and after. It is kept as a low-risk residual safeguard
against any localised artifact rather than as the fix for the disc-wide bias, which
would need retuning verified against a fresh render. Clamping colour globally would
also wash out real vegetation and shallow water, so the excess is compared against
its own locally-blurred background instead: a broad region matches its background
and is left alone, while a localised spike stands far above it and is pulled down.
Args:
rgb: `(H, W, 3)` uint8 image.
tol: How far above the local background the excess may sit before clamping.
sigma: Blur radius, in pixels, defining "local background".
strength: How much of the anomaly to remove, 0..1.
Returns:
`(H, W, 3)` uint8 image with the seam neutralised.
"""
out = rgb.astype(np.float32).copy()
red = out[..., 0]
for chan in (1, 2):
excess = out[..., chan] - red
ridge = excess - gaussian_filter(excess, sigma=sigma)
out[..., chan] -= strength * np.clip(ridge - tol, 0.0, None)
return np.clip(out, 0, 255).astype(np.uint8)
def soften(rgb, alpha, feather=14.0, glow=0.55, glow_rgb=(0.60, 0.78, 1.0), trim=4):
"""Soften the hard edge where FCI's scan coverage stops, and add a limb glow.
pyramids-eo's true_color_with_night_ir(keep_alpha=True) derives a binary 0/1 coverage
band, so the imagery ends in a razor cut against black. Two passes fix it: the feather ramps brightness down over the last
few pixels of real data, and the glow adds the atmospheric halo that real full-disk
imagery carries, because the line of sight grazes the atmosphere at the horizon.
Neither pass invents imagery -- both only attenuate or tint pixels that already carry
data.
Args:
rgb: `(H, W, 3)` uint8 image.
alpha: `(H, W)` binary coverage mask (0 outside the scan, 255 inside).
feather: Gaussian sigma, in pixels, of the fade-out ramp.
glow: Strength of the atmospheric halo; 0 disables it.
glow_rgb: Colour of the halo.
trim: Pixels of the outermost, ragged rim to drop before feathering.
Returns:
`(H, W, 3)` uint8 image with a soft boundary.
"""
data = alpha > 0
if trim:
data = binary_erosion(data, iterations=int(trim))
data = data.astype(np.float32)
ramp = np.clip(gaussian_filter(data, sigma=feather), 0.0, 1.0)
out = rgb.astype(np.float32) * ramp[..., None]
if glow > 0:
limb = (
np.clip(
gaussian_filter(data, sigma=feather * 0.55)
- gaussian_filter(data, sigma=feather * 1.9),
0.0,
1.0,
)
* data
)
peak = float(limb.max())
if peak > 1e-6:
limb = np.clip(limb / peak, 0.0, 1.0) ** 1.3 # tighten the rim
out += (
glow * 255.0 * limb[..., None] * np.asarray(glow_rgb, dtype=np.float32)
)
return np.clip(out, 0, 255).astype(np.uint8)
def soften_file(path, **kwargs):
"""Read an RGBA GeoTIFF frame, clean its seam and boundary, write it back in place."""
ds = Dataset.read_file(path, read_only=False)
bands = ds.read_array()
# read_array returns (bands, rows, cols), so `bands[:3]` slices the band
# axis; a single-band read comes back 2-D and would slice rows instead.
if bands.ndim != 3:
raise ValueError(f"expected an RGBA frame, got shape {bands.shape}")
rgb = np.dstack(bands[:3])
alpha = (
bands[3] if ds.band_count > 3 else (rgb.max(axis=2) > 0).astype(np.uint8) * 255
)
out = soften(declip_cyan(rgb), alpha, **kwargs)
for i in range(3):
ds.write_array(out[..., i], band=i)
ds.close() # flush the in-place writes before the next frame reads them
return path
In [ ]:
Copied!
def rendered(path):
"""True when `path` holds a complete frame; an interrupted write leaves a stub."""
return path.exists() and path.stat().st_size > (1 << 20)
def load_channel(chunks, channel):
"""Read one FCI channel onto the common 2 km grid.
Reads one channel at a time with `read_fci_l1c(chunks, channel)` rather than
all five via `channels=` in a single call: each channel's `Dataset` is
materialised, downsampled and freed (`del` + `gc.collect()`) before the next
one is read, bounding peak memory to one channel at a time -- the four 1 km
solar bands held simultaneously would roughly quadruple it. `channels=`
trades that for fewer chunk-decode passes; not taken here deliberately.
Args:
chunks: The cycle's chunk files.
channel: Channel name, e.g. `"vis_06"`.
Returns:
An `(array, geotransform, crs)` triple on the 2 km grid.
"""
dataset = read_fci_l1c(chunks, channel)
array = np.asarray(dataset.read_array(), dtype="float32")
geo, crs = list(dataset.geotransform), str(dataset.crs)
del dataset
gc.collect()
if (
abs(geo[1]) < NATIVE_1KM_THRESHOLD_M
): # average the 1 km solar bands down to 2 km
rows, cols = array.shape[0] // 2 * 2, array.shape[1] // 2 * 2
with np.errstate(invalid="ignore"):
array = np.nanmean(
array[:rows, :cols].reshape(rows // 2, 2, cols // 2, 2), axis=(1, 3)
).astype("float32")
geo[1] *= 2
geo[5] *= 2
return array, tuple(geo), crs
def build_composite(chunks, when):
"""Build the day/night composite for one cycle, on the geostationary grid.
Sun-zenith corrected, Rayleigh corrected true colour by day; grey cloud tops over
Black Marble city lights by night; the two cross-faded by solar zenith angle.
Args:
chunks: The cycle's chunk files.
when: The cycle's nominal time, in UTC.
Returns:
A `(frame, geotransform, crs)` triple, where `frame` is `(4, H, W)` RGBA in 0..1.
"""
# geo / crs are captured per channel and taken from IR_CHANNEL specifically --
# not "whichever channel the loop reads last" -- so reordering SOLAR_CHANNELS can
# never silently change the grid every downstream computation is built on.
bands, geos, crses = {}, {}, {}
for channel in SOLAR_CHANNELS + [IR_CHANNEL]:
bands[channel], geos[channel], crses[channel] = load_channel(chunks, channel)
geo, crs = geos[IR_CHANNEL], crses[IR_CHANNEL]
ir_shape = bands[IR_CHANNEL].shape
for name, arr in bands.items():
assert arr.shape == ir_shape, (
f"{name} grid {arr.shape} does not match {IR_CHANNEL} grid {ir_shape}"
)
lon, lat = disc_lonlat(geo, crs, ir_shape)
# Solar and viewing geometry, both from pyramids-eo.
with np.errstate(invalid="ignore"):
sza, sun_azimuth = solar_zenith_azimuth(when, lat=lat, lon=lon)
view_zenith, sat_azimuth = satellite_zenith_azimuth(lat=lat, lon=lon, sat_lon=0.0)
azidiff = relative_azimuth(sun_azimuth, sat_azimuth)
del sun_azimuth, sat_azimuth
gc.collect()
# Sun-zenith handling, the pair satpy applies as `sunz_corrected` then
# `sunz_reduced`: divide out cos(sza) so the disc does not darken toward the
# terminator, then taper the signal back down where illumination is marginal, so
# deep shadow -- the eclipse umbra above all -- reads dark instead of grey.
for channel in SOLAR_CHANNELS:
corrected = sunz_correct(bands[channel], sza)
bands[channel] = np.asarray(sunz_reduce(corrected, sza), dtype="float32")
# Rayleigh correction. Molecular scattering goes as lambda^-4, so the blue band is
# corrected several times as hard as the red one; the NIR is left alone, as satpy
# also leaves it.
for channel, wavelength in RAYLEIGH_UM.items():
bands[channel] = np.asarray(
rayleigh_correct(
bands[channel],
wavelength_um=wavelength,
sza=sza,
vza=view_zenith,
azidiff=azidiff,
),
dtype="float32",
)
del view_zenith, azidiff
gc.collect()
# FCI's 0.51 um band deliberately misses chlorophyll, so vegetation reads brown unless
# some NIR is mixed back in -- hence the NDVI hybrid green rather than the raw band.
day = true_color(
bands["vis_06"],
bands["vis_04"],
bands["vis_08"],
green=bands["vis_05"],
green_mode="ndvi_hybrid",
clip=True,
)
# Grey cloud tops rather than a three-channel IR false colour, whose cold tops come
# out cyan; the alpha lets the city lights show through the gaps.
brightness_temperature = bands[IR_CHANNEL]
opacity = np.clip((273.0 - brightness_temperature) / 70.0, 0.0, 1.0)
opacity = np.where(np.isfinite(brightness_temperature), opacity, 0.0)
grey = np.clip((290.0 - brightness_temperature) / 90.0, 0.0, 1.0)
clouds = night_ir(0.92 * grey, 0.92 * grey, 0.92 * grey, alpha=opacity)
like = Dataset.from_array(
brightness_temperature,
no_data_value=np.nan,
geo_ref=GeoReference(geo=geo, epsg=crs),
)
city_lights = np.nan_to_num(
np.asarray(
static_image(BLACK_MARBLE, like=like, cache_dir=AUX).read_array(),
dtype="float32",
),
nan=0.0,
)
if city_lights.max() > 1.5: # Black Marble arrives as 0-255 bytes
city_lights = city_lights / 255.0
frame = np.asarray(
true_color_with_night_ir(
np.asarray(day),
np.asarray(clouds),
city_lights,
sza,
lim_low=78,
lim_high=88,
keep_alpha=True,
)
)
del bands, day, clouds, city_lights, like
gc.collect()
return frame, geo, crs
def render_cycle(slot):
"""Fetch one FCI cycle (kept in the cache dir) and write its nsper + tilted frames.
Resampling is cached in the `frames_*_raw` directories; the edge treatment always
re-runs, writing its own copies into `frames_nsper` / `frames_tpers`.
"""
tag = slot.strftime("%Y%m%d%H%M")
raw = ((NSPER, NSPER_RAW / f"{tag}.tif"), (TPERS, TPERS_RAW / f"{tag}.tif"))
final = (NSPER_DIR / f"{tag}.tif", TPERS_DIR / f"{tag}.tif")
if not all(rendered(path) for _, path in raw):
work = RAW_DIR / tag
if not list(work.rglob("*.nc")): # not cached yet -> fetch once
work.mkdir(parents=True, exist_ok=True)
got = EarthLens(
data_source="eumetsat",
variables={DATASET: []}, # whole product, no band subset
start=slot.strftime("%Y-%m-%d %H:%M"),
end=(slot + pd.Timedelta(minutes=9)).strftime("%Y-%m-%d %H:%M"),
fmt="%Y-%m-%d %H:%M",
path=work,
**DISC,
).download(progress_bar=False)
if not got:
return None # occasional missing cycle
for archive in got: # FCI ships as a zipped product
with zipfile.ZipFile(archive) as bundle:
bundle.extractall(work)
archive.unlink()
chunks = sorted(str(p) for p in work.rglob("*.nc"))
frame, geo, crs = build_composite(chunks, slot.to_pydatetime())
composite = Dataset.from_array(
frame, no_data_value=np.nan, geo_ref=GeoReference(geo=geo, epsg=crs)
)
for area, path in raw:
# Warp in float and tone-map afterwards: one quantisation instead of two.
framed = to_area(
composite,
area["crs"],
area["extent"],
area["width"],
area["height"],
method="bilinear",
)
eight_bit = np.asarray(
stretch(
np.asarray(framed.read_array(), dtype="float32"),
# Tuned against the satpy frames: max_stretch sets overall
# brightness, gamma the lift on the low end. A high gamma turned
# the near-zero signal inside the umbra into a visible green cast,
# so the lift comes down and brightness is recovered via max_stretch.
kind="crude",
min_stretch=0.0,
max_stretch=0.75,
gamma=1.4,
preserve_alpha=True,
)
)
Dataset.from_array(
eight_bit,
# A scalar 0 across all four RGBA bands would flag any
# genuinely
# covered, genuinely dark pixel as missing -- measured on this
# pipeline's own frames: 5% of on-disc pixels have R == 0 during the
# umbra's deepest cycles. Passing None for R/G/B does not leave them
# nodata-free: pyramids-gis's SetNoDataValue(None) call fails on a
# Byte band and falls back to NaN, so they end up tagged
# a NaN `GDAL_NODATA` instead. That is harmless here -- no Byte pixel is
# ever exactly NaN, so it behaves like "no nodata" for every real
# value -- but it is real metadata a stricter reader would see,
# not an omitted one. Alpha still carries the real coverage signal
# (0 off-disc, ~255 on) for anything that wants it.
no_data_value=[None, None, None, 0],
geo_ref=GeoReference(geo=framed.geotransform, epsg=str(framed.crs)),
).to_file(path)
del frame, composite
gc.collect()
for (_, src), dst in zip(raw, final):
shutil.copyfile(src, dst)
soften_file(dst)
return final
def rendered(path):
"""True when `path` holds a complete frame; an interrupted write leaves a stub."""
return path.exists() and path.stat().st_size > (1 << 20)
def load_channel(chunks, channel):
"""Read one FCI channel onto the common 2 km grid.
Reads one channel at a time with `read_fci_l1c(chunks, channel)` rather than
all five via `channels=` in a single call: each channel's `Dataset` is
materialised, downsampled and freed (`del` + `gc.collect()`) before the next
one is read, bounding peak memory to one channel at a time -- the four 1 km
solar bands held simultaneously would roughly quadruple it. `channels=`
trades that for fewer chunk-decode passes; not taken here deliberately.
Args:
chunks: The cycle's chunk files.
channel: Channel name, e.g. `"vis_06"`.
Returns:
An `(array, geotransform, crs)` triple on the 2 km grid.
"""
dataset = read_fci_l1c(chunks, channel)
array = np.asarray(dataset.read_array(), dtype="float32")
geo, crs = list(dataset.geotransform), str(dataset.crs)
del dataset
gc.collect()
if (
abs(geo[1]) < NATIVE_1KM_THRESHOLD_M
): # average the 1 km solar bands down to 2 km
rows, cols = array.shape[0] // 2 * 2, array.shape[1] // 2 * 2
with np.errstate(invalid="ignore"):
array = np.nanmean(
array[:rows, :cols].reshape(rows // 2, 2, cols // 2, 2), axis=(1, 3)
).astype("float32")
geo[1] *= 2
geo[5] *= 2
return array, tuple(geo), crs
def build_composite(chunks, when):
"""Build the day/night composite for one cycle, on the geostationary grid.
Sun-zenith corrected, Rayleigh corrected true colour by day; grey cloud tops over
Black Marble city lights by night; the two cross-faded by solar zenith angle.
Args:
chunks: The cycle's chunk files.
when: The cycle's nominal time, in UTC.
Returns:
A `(frame, geotransform, crs)` triple, where `frame` is `(4, H, W)` RGBA in 0..1.
"""
# geo / crs are captured per channel and taken from IR_CHANNEL specifically --
# not "whichever channel the loop reads last" -- so reordering SOLAR_CHANNELS can
# never silently change the grid every downstream computation is built on.
bands, geos, crses = {}, {}, {}
for channel in SOLAR_CHANNELS + [IR_CHANNEL]:
bands[channel], geos[channel], crses[channel] = load_channel(chunks, channel)
geo, crs = geos[IR_CHANNEL], crses[IR_CHANNEL]
ir_shape = bands[IR_CHANNEL].shape
for name, arr in bands.items():
assert arr.shape == ir_shape, (
f"{name} grid {arr.shape} does not match {IR_CHANNEL} grid {ir_shape}"
)
lon, lat = disc_lonlat(geo, crs, ir_shape)
# Solar and viewing geometry, both from pyramids-eo.
with np.errstate(invalid="ignore"):
sza, sun_azimuth = solar_zenith_azimuth(when, lat=lat, lon=lon)
view_zenith, sat_azimuth = satellite_zenith_azimuth(lat=lat, lon=lon, sat_lon=0.0)
azidiff = relative_azimuth(sun_azimuth, sat_azimuth)
del sun_azimuth, sat_azimuth
gc.collect()
# Sun-zenith handling, the pair satpy applies as `sunz_corrected` then
# `sunz_reduced`: divide out cos(sza) so the disc does not darken toward the
# terminator, then taper the signal back down where illumination is marginal, so
# deep shadow -- the eclipse umbra above all -- reads dark instead of grey.
for channel in SOLAR_CHANNELS:
corrected = sunz_correct(bands[channel], sza)
bands[channel] = np.asarray(sunz_reduce(corrected, sza), dtype="float32")
# Rayleigh correction. Molecular scattering goes as lambda^-4, so the blue band is
# corrected several times as hard as the red one; the NIR is left alone, as satpy
# also leaves it.
for channel, wavelength in RAYLEIGH_UM.items():
bands[channel] = np.asarray(
rayleigh_correct(
bands[channel],
wavelength_um=wavelength,
sza=sza,
vza=view_zenith,
azidiff=azidiff,
),
dtype="float32",
)
del view_zenith, azidiff
gc.collect()
# FCI's 0.51 um band deliberately misses chlorophyll, so vegetation reads brown unless
# some NIR is mixed back in -- hence the NDVI hybrid green rather than the raw band.
day = true_color(
bands["vis_06"],
bands["vis_04"],
bands["vis_08"],
green=bands["vis_05"],
green_mode="ndvi_hybrid",
clip=True,
)
# Grey cloud tops rather than a three-channel IR false colour, whose cold tops come
# out cyan; the alpha lets the city lights show through the gaps.
brightness_temperature = bands[IR_CHANNEL]
opacity = np.clip((273.0 - brightness_temperature) / 70.0, 0.0, 1.0)
opacity = np.where(np.isfinite(brightness_temperature), opacity, 0.0)
grey = np.clip((290.0 - brightness_temperature) / 90.0, 0.0, 1.0)
clouds = night_ir(0.92 * grey, 0.92 * grey, 0.92 * grey, alpha=opacity)
like = Dataset.from_array(
brightness_temperature,
no_data_value=np.nan,
geo_ref=GeoReference(geo=geo, epsg=crs),
)
city_lights = np.nan_to_num(
np.asarray(
static_image(BLACK_MARBLE, like=like, cache_dir=AUX).read_array(),
dtype="float32",
),
nan=0.0,
)
if city_lights.max() > 1.5: # Black Marble arrives as 0-255 bytes
city_lights = city_lights / 255.0
frame = np.asarray(
true_color_with_night_ir(
np.asarray(day),
np.asarray(clouds),
city_lights,
sza,
lim_low=78,
lim_high=88,
keep_alpha=True,
)
)
del bands, day, clouds, city_lights, like
gc.collect()
return frame, geo, crs
def render_cycle(slot):
"""Fetch one FCI cycle (kept in the cache dir) and write its nsper + tilted frames.
Resampling is cached in the `frames_*_raw` directories; the edge treatment always
re-runs, writing its own copies into `frames_nsper` / `frames_tpers`.
"""
tag = slot.strftime("%Y%m%d%H%M")
raw = ((NSPER, NSPER_RAW / f"{tag}.tif"), (TPERS, TPERS_RAW / f"{tag}.tif"))
final = (NSPER_DIR / f"{tag}.tif", TPERS_DIR / f"{tag}.tif")
if not all(rendered(path) for _, path in raw):
work = RAW_DIR / tag
if not list(work.rglob("*.nc")): # not cached yet -> fetch once
work.mkdir(parents=True, exist_ok=True)
got = EarthLens(
data_source="eumetsat",
variables={DATASET: []}, # whole product, no band subset
start=slot.strftime("%Y-%m-%d %H:%M"),
end=(slot + pd.Timedelta(minutes=9)).strftime("%Y-%m-%d %H:%M"),
fmt="%Y-%m-%d %H:%M",
path=work,
**DISC,
).download(progress_bar=False)
if not got:
return None # occasional missing cycle
for archive in got: # FCI ships as a zipped product
with zipfile.ZipFile(archive) as bundle:
bundle.extractall(work)
archive.unlink()
chunks = sorted(str(p) for p in work.rglob("*.nc"))
frame, geo, crs = build_composite(chunks, slot.to_pydatetime())
composite = Dataset.from_array(
frame, no_data_value=np.nan, geo_ref=GeoReference(geo=geo, epsg=crs)
)
for area, path in raw:
# Warp in float and tone-map afterwards: one quantisation instead of two.
framed = to_area(
composite,
area["crs"],
area["extent"],
area["width"],
area["height"],
method="bilinear",
)
eight_bit = np.asarray(
stretch(
np.asarray(framed.read_array(), dtype="float32"),
# Tuned against the satpy frames: max_stretch sets overall
# brightness, gamma the lift on the low end. A high gamma turned
# the near-zero signal inside the umbra into a visible green cast,
# so the lift comes down and brightness is recovered via max_stretch.
kind="crude",
min_stretch=0.0,
max_stretch=0.75,
gamma=1.4,
preserve_alpha=True,
)
)
Dataset.from_array(
eight_bit,
# A scalar 0 across all four RGBA bands would flag any
# genuinely
# covered, genuinely dark pixel as missing -- measured on this
# pipeline's own frames: 5% of on-disc pixels have R == 0 during the
# umbra's deepest cycles. Passing None for R/G/B does not leave them
# nodata-free: pyramids-gis's SetNoDataValue(None) call fails on a
# Byte band and falls back to NaN, so they end up tagged
# a NaN `GDAL_NODATA` instead. That is harmless here -- no Byte pixel is
# ever exactly NaN, so it behaves like "no nodata" for every real
# value -- but it is real metadata a stricter reader would see,
# not an omitted one. Alpha still carries the real coverage signal
# (0 off-disc, ~255 on) for anything that wants it.
no_data_value=[None, None, None, 0],
geo_ref=GeoReference(geo=framed.geotransform, epsg=str(framed.crs)),
).to_file(path)
del frame, composite
gc.collect()
for (_, src), dst in zip(raw, final):
shutil.copyfile(src, dst)
soften_file(dst)
return final
Fetch each 10-minute cycle and render both frames¶
One cycle (~770 MB) is downloaded and both frames are written. Only the zip archive is removed — the extracted chunks stay in the cache so a re-run costs no downloads. Cache-aware: a slot whose frames exist is skipped. Across the 16 cycles this leaves roughly 12 GB under the earthlens cache directory; delete <cache dir>/eclipse_fci when you are done with it.
In [ ]:
Copied!
nsper_frames, tpers_frames = [], []
for slot in SLOTS:
result = render_cycle(slot)
if result is None:
continue # occasional missing cycle -> skip the frame
nsper_frames.append(result[0])
tpers_frames.append(result[1])
print(f"{len(nsper_frames)} frames: {SLOTS[0]:%H:%M} -> {SLOTS[-1]:%H:%M} UTC")
nsper_frames, tpers_frames = [], []
for slot in SLOTS:
result = render_cycle(slot)
if result is None:
continue # occasional missing cycle -> skip the frame
nsper_frames.append(result[0])
tpers_frames.append(result[1])
print(f"{len(nsper_frames)} frames: {SLOTS[0]:%H:%M} -> {SLOTS[-1]:%H:%M} UTC")
The cinematic — the Moon's shadow from space¶
Near-side perspective centred on the North Atlantic, one frame every 10 minutes (fps 3). The penumbra darkens the middle of the disk around greatest eclipse, then lifts.
Branding the frames¶
The satellite name is the fixed identity so it takes the upper line, with the per-frame clock beneath it; top-left is the cleanest text zone in these clips, though not clean enough for bare white text, so both lines get an outline and a shared scrim. The earthlens mark goes bottom-right at 11% of the frame width -- the largest size that stays clear of the Sun in the geometry animation.
In [ ]:
Copied!
BRAND_DIR = Path("../../../_images/branding/earthlens-brand-kit/logo")
BRAND_MARK = BRAND_DIR / "earthlens-lockup-stacked-overlay-full.png"
# the plated variant: it carries its own navy scrim, which the transparent one needs over
# busy imagery -- under the Meteosat mark the background runs mean 99, std 84, p90 254
LOGO_FRAC = 0.18 # of figure width
LOGO_MARGIN_Y = 0.0 # of figure height -- tucked low, clear of the imagery
LOGO_MARGIN = 0.025
SCRIM = {"facecolor": "#091c30", "alpha": 0.55, "edgecolor": "none"}
def stamp_titles(
fig,
label,
satellite,
sizes=(16, 14),
x=0.02,
y=0.965,
step=0.048,
clock_sample="Date = 00:00 UTC",
):
"""Put the satellite name above the frame label, both on a shared scrim.
The name is the fixed identity so it takes the upper line; the clock changes every frame
and sits beneath it. Top-left is the cleanest text zone in these clips, but not clean
enough for bare white text, so both lines get an outline and a scrim behind them.
Everything is added to the axes rather than the figure. The clock is an axes-level artist,
and matplotlib draws figure-level artists after every axes -- so a figure-level scrim would
paint over the clock whatever zorder it carried, dimming it to the scrim's alpha.
Args:
fig: The figure being animated.
label: The frame-label `Text` artist `animate` created.
satellite: Identity line, e.g. `Meteosat-12 - MTG-I1 / FCI`.
sizes: Font sizes of the identity line and the clock, in points.
x: Left edge, in figure coordinates.
y: Top edge of the identity line, in figure coordinates.
step: Line spacing, in figure coordinates.
clock_sample: Stand-in used to size the scrim; `animate` leaves the clock empty until
a frame is drawn, so measuring it directly would size the scrim for one line.
Returns:
Text: The identity-line artist.
"""
ax = label.axes
stroke = [pe.withStroke(linewidth=3, foreground="black")]
top = ax.transAxes.inverted().transform(fig.transFigure.transform((x, y)))
below = ax.transAxes.inverted().transform(fig.transFigure.transform((x, y - step)))
name = ax.text(
top[0],
top[1],
satellite,
transform=ax.transAxes,
ha="left",
va="top",
color="white",
fontsize=sizes[0],
zorder=5,
path_effects=stroke,
)
label.set_position(tuple(below))
label.set_ha("left")
label.set_va("top")
label.set_size(sizes[1])
label.set_path_effects(stroke)
label.set_zorder(5)
was = label.get_text()
label.set_text(clock_sample)
fig.canvas.draw() # text extents are only known once the figure has been laid out
boxes = [t.get_window_extent() for t in (name, label)]
label.set_text(was)
pad = 0.30 * min(b.height for b in boxes)
x0, y0 = min(b.x0 for b in boxes) - pad, min(b.y0 for b in boxes) - pad
x1, y1 = max(b.x1 for b in boxes) + pad, max(b.y1 for b in boxes) + pad
(ax0, ay0), (ax1, ay1) = ax.transAxes.inverted().transform([(x0, y0), (x1, y1)])
ax.add_patch(
FancyBboxPatch(
(ax0, ay0),
ax1 - ax0,
ay1 - ay0,
boxstyle="round,pad=0.003,rounding_size=0.010",
transform=ax.transAxes,
zorder=4,
clip_on=False,
**SCRIM,
)
)
return name
# cleopatra 0.31 counts the in-domain cells of every ArrayGlyph by materialising one Python
# tuple per cell, and for a 4-D (RGB) stack it does that over the whole stack rather than one
# frame. A 25-frame full-disk animation needs ~15 GB of tuples for a number the library never
# reads, so it dies with MemoryError. Report the count directly instead -- same value, no list.
# Upstream: https://github.com/serapeum-org/cleopatra/issues/304
def _domain_cell_count(arr, mask=None):
"""Stand-in for `get_indices2` whose length is the domain-cell count.
`mask` is accepted only for signature compatibility with the real
`get_indices2` -- cleopatra's own call sites always pass `[np.nan]`,
which is exactly what this counts, so it is safe to ignore rather than
honour a caller-supplied sentinel list. `arr` is one frame, `(bands, H,
W)` for a multi-band array or `(H, W)` for a single-band one; band 0
stands in for the whole frame, matching what cleopatra passes here.
"""
frame = arr[0] if arr.ndim >= 3 else arr
return range(int(np.count_nonzero(~np.isnan(frame))))
# Installed for the rest of this kernel, not scoped to one cell -- the animation cell
# below also builds a DatasetCollection and needs the same patch. Not restored.
_array_glyph.get_indices2 = _domain_cell_count
BRAND_DIR = Path("../../../_images/branding/earthlens-brand-kit/logo")
BRAND_MARK = BRAND_DIR / "earthlens-lockup-stacked-overlay-full.png"
# the plated variant: it carries its own navy scrim, which the transparent one needs over
# busy imagery -- under the Meteosat mark the background runs mean 99, std 84, p90 254
LOGO_FRAC = 0.18 # of figure width
LOGO_MARGIN_Y = 0.0 # of figure height -- tucked low, clear of the imagery
LOGO_MARGIN = 0.025
SCRIM = {"facecolor": "#091c30", "alpha": 0.55, "edgecolor": "none"}
def stamp_titles(
fig,
label,
satellite,
sizes=(16, 14),
x=0.02,
y=0.965,
step=0.048,
clock_sample="Date = 00:00 UTC",
):
"""Put the satellite name above the frame label, both on a shared scrim.
The name is the fixed identity so it takes the upper line; the clock changes every frame
and sits beneath it. Top-left is the cleanest text zone in these clips, but not clean
enough for bare white text, so both lines get an outline and a scrim behind them.
Everything is added to the axes rather than the figure. The clock is an axes-level artist,
and matplotlib draws figure-level artists after every axes -- so a figure-level scrim would
paint over the clock whatever zorder it carried, dimming it to the scrim's alpha.
Args:
fig: The figure being animated.
label: The frame-label `Text` artist `animate` created.
satellite: Identity line, e.g. `Meteosat-12 - MTG-I1 / FCI`.
sizes: Font sizes of the identity line and the clock, in points.
x: Left edge, in figure coordinates.
y: Top edge of the identity line, in figure coordinates.
step: Line spacing, in figure coordinates.
clock_sample: Stand-in used to size the scrim; `animate` leaves the clock empty until
a frame is drawn, so measuring it directly would size the scrim for one line.
Returns:
Text: The identity-line artist.
"""
ax = label.axes
stroke = [pe.withStroke(linewidth=3, foreground="black")]
top = ax.transAxes.inverted().transform(fig.transFigure.transform((x, y)))
below = ax.transAxes.inverted().transform(fig.transFigure.transform((x, y - step)))
name = ax.text(
top[0],
top[1],
satellite,
transform=ax.transAxes,
ha="left",
va="top",
color="white",
fontsize=sizes[0],
zorder=5,
path_effects=stroke,
)
label.set_position(tuple(below))
label.set_ha("left")
label.set_va("top")
label.set_size(sizes[1])
label.set_path_effects(stroke)
label.set_zorder(5)
was = label.get_text()
label.set_text(clock_sample)
fig.canvas.draw() # text extents are only known once the figure has been laid out
boxes = [t.get_window_extent() for t in (name, label)]
label.set_text(was)
pad = 0.30 * min(b.height for b in boxes)
x0, y0 = min(b.x0 for b in boxes) - pad, min(b.y0 for b in boxes) - pad
x1, y1 = max(b.x1 for b in boxes) + pad, max(b.y1 for b in boxes) + pad
(ax0, ay0), (ax1, ay1) = ax.transAxes.inverted().transform([(x0, y0), (x1, y1)])
ax.add_patch(
FancyBboxPatch(
(ax0, ay0),
ax1 - ax0,
ay1 - ay0,
boxstyle="round,pad=0.003,rounding_size=0.010",
transform=ax.transAxes,
zorder=4,
clip_on=False,
**SCRIM,
)
)
return name
# cleopatra 0.31 counts the in-domain cells of every ArrayGlyph by materialising one Python
# tuple per cell, and for a 4-D (RGB) stack it does that over the whole stack rather than one
# frame. A 25-frame full-disk animation needs ~15 GB of tuples for a number the library never
# reads, so it dies with MemoryError. Report the count directly instead -- same value, no list.
# Upstream: https://github.com/serapeum-org/cleopatra/issues/304
def _domain_cell_count(arr, mask=None):
"""Stand-in for `get_indices2` whose length is the domain-cell count.
`mask` is accepted only for signature compatibility with the real
`get_indices2` -- cleopatra's own call sites always pass `[np.nan]`,
which is exactly what this counts, so it is safe to ignore rather than
honour a caller-supplied sentinel list. `arr` is one frame, `(bands, H,
W)` for a multi-band array or `(H, W)` for a single-band one; band 0
stands in for the whole frame, matching what cleopatra passes here.
"""
frame = arr[0] if arr.ndim >= 3 else arr
return range(int(np.count_nonzero(~np.isnan(frame))))
# Installed for the rest of this kernel, not scoped to one cell -- the animation cell
# below also builds a DatasetCollection and needs the same patch. Not restored.
_array_glyph.get_indices2 = _domain_cell_count
In [ ]:
Copied!
dc = DatasetCollection.from_files(
[str(p) for p in nsper_frames], date_format="%Y%m%d%H%M", date_regex=r"\d{12}"
)
labels = [t.strftime("%H:%M UTC") for t in dc.time]
glyph = dc.plot(
rgb_options={"rgb": [0, 1, 2], "surface_reflectance": 255},
figsize=(10, 10),
full_bleed=True,
)
glyph.animate(
labels,
playback=Animation(interval=350, frame_label=FrameLabel(color="white", size=20)),
full_bleed=True,
)
# `full_bleed` draws no title and no colorbar, so the frame label is the only text on the axes.
(label,) = glyph.fig.axes[0].texts
stamp_titles(glyph.fig, label, "Meteosat-12 · MTG-I1 / FCI")
glyph.stamp_mark(BRAND_MARK, frac=LOGO_FRAC, margin=(LOGO_MARGIN, LOGO_MARGIN_Y))
mp4 = OUT / "eclipse_fci_nsper.mp4"
glyph.save_animation(str(mp4), fps=3, dpi=150)
gif = str(OUT / "eclipse_fci_nsper.gif")
glyph.fig.set_dpi(90)
glyph.save_animation(gif, fps=3)
plt.close("all")
del dc, glyph
gc.collect()
display(Image(filename=gif))
dc = DatasetCollection.from_files(
[str(p) for p in nsper_frames], date_format="%Y%m%d%H%M", date_regex=r"\d{12}"
)
labels = [t.strftime("%H:%M UTC") for t in dc.time]
glyph = dc.plot(
rgb_options={"rgb": [0, 1, 2], "surface_reflectance": 255},
figsize=(10, 10),
full_bleed=True,
)
glyph.animate(
labels,
playback=Animation(interval=350, frame_label=FrameLabel(color="white", size=20)),
full_bleed=True,
)
# `full_bleed` draws no title and no colorbar, so the frame label is the only text on the axes.
(label,) = glyph.fig.axes[0].texts
stamp_titles(glyph.fig, label, "Meteosat-12 · MTG-I1 / FCI")
glyph.stamp_mark(BRAND_MARK, frac=LOGO_FRAC, margin=(LOGO_MARGIN, LOGO_MARGIN_Y))
mp4 = OUT / "eclipse_fci_nsper.mp4"
glyph.save_animation(str(mp4), fps=3, dpi=150)
gif = str(OUT / "eclipse_fci_nsper.gif")
glyph.fig.set_dpi(90)
glyph.save_animation(gif, fps=3)
plt.close("all")
del dc, glyph
gc.collect()
display(Image(filename=gif))
The tilted view — the shadow riding the limb¶
The same cycles in a tilted near-side perspective (PROJ tpers): a low-angle, curved-limb view like the
one a crew in orbit would get, with Iberia and the Sahara in the foreground and the Moon's shadow sweeping
the northern limb. Because it is rendered from the full disk, the limb and the edge of the atmosphere
are real, not a crop.
In [ ]:
Copied!
dc2 = DatasetCollection.from_files(
[str(p) for p in tpers_frames], date_format="%Y%m%d%H%M", date_regex=r"\d{12}"
)
labels2 = [t.strftime("%H:%M UTC") for t in dc2.time]
glyph2 = dc2.plot(
rgb_options={"rgb": [0, 1, 2], "surface_reflectance": 255},
figsize=(12, 8),
full_bleed=True,
)
glyph2.animate(
labels2,
playback=Animation(interval=350, frame_label=FrameLabel(color="white", size=18)),
full_bleed=True,
)
# `full_bleed` draws no title and no colorbar, so the frame label is the only text on the axes.
(label2,) = glyph2.fig.axes[0].texts
stamp_titles(glyph2.fig, label2, "Meteosat-12 · MTG-I1 / FCI")
glyph2.stamp_mark(BRAND_MARK, frac=LOGO_FRAC, margin=(LOGO_MARGIN, LOGO_MARGIN_Y))
mp4b = OUT / "eclipse_fci_tpers.mp4"
glyph2.save_animation(str(mp4b), fps=3, dpi=150)
gif2 = str(OUT / "eclipse_fci_tpers.gif")
glyph2.fig.set_dpi(90)
glyph2.save_animation(gif2, fps=3)
plt.close("all")
del dc2, glyph2
gc.collect()
display(Image(filename=gif2))
dc2 = DatasetCollection.from_files(
[str(p) for p in tpers_frames], date_format="%Y%m%d%H%M", date_regex=r"\d{12}"
)
labels2 = [t.strftime("%H:%M UTC") for t in dc2.time]
glyph2 = dc2.plot(
rgb_options={"rgb": [0, 1, 2], "surface_reflectance": 255},
figsize=(12, 8),
full_bleed=True,
)
glyph2.animate(
labels2,
playback=Animation(interval=350, frame_label=FrameLabel(color="white", size=18)),
full_bleed=True,
)
# `full_bleed` draws no title and no colorbar, so the frame label is the only text on the axes.
(label2,) = glyph2.fig.axes[0].texts
stamp_titles(glyph2.fig, label2, "Meteosat-12 · MTG-I1 / FCI")
glyph2.stamp_mark(BRAND_MARK, frac=LOGO_FRAC, margin=(LOGO_MARGIN, LOGO_MARGIN_Y))
mp4b = OUT / "eclipse_fci_tpers.mp4"
glyph2.save_animation(str(mp4b), fps=3, dpi=150)
gif2 = str(OUT / "eclipse_fci_tpers.gif")
glyph2.fig.set_dpi(90)
glyph2.save_animation(gif2, fps=3)
plt.close("all")
del dc2, glyph2
gc.collect()
display(Image(filename=gif2))
Both MP4s are written under out/eclipse_fci/ (eclipse_fci_nsper.mp4, eclipse_fci_tpers.mp4) for sharing, and the raw granules stay in the earthlens cache dir so a re-run costs no downloads.