Every satellite earthlens can reach, in orbit¶
earthlens talks to 61 data providers. A good number of them are not satellites at all — ground station networks, reanalysis models, vector basemaps — but behind the rest sits a fleet of spacecraft, and this notebook draws it.
No Sun, no Moon: just Earth and the satellites, each on its own orbit, moving at its own speed.
The point is the structure. Two populations dominate, and they could hardly be less alike:
- a tight, fast shell of low-Earth orbiters a few hundred kilometres up, crossing the poles roughly every 100 minutes, sweeping a new strip of ground on every pass;
- a sparse, distant ring of geostationary satellites almost 36,000 km out, each one parked over a fixed longitude, turning at exactly Earth's own rate so it never moves in the frame.
Watch the geostationary ring: it looks frozen relative to the surface below it. That is not an artefact — it is what "geostationary" means. The altitude is what does it: Kepler's third law turns 35,786 km into a period of one sidereal day, which is precisely the rate Earth turns at.
VIIRS Earth at Night¶
Linked rather than embedded. The clip is committed under docs/_images/animation/, which the
documentation site serves, so the picture survives this repo's pre-commit output stripping without
carrying a base64 copy of itself in the notebook.
NASA Blue Marble¶
Linked rather than embedded. The clip is committed under docs/_images/animation/, which the
documentation site serves, so the picture survives this repo's pre-commit output stripping without
carrying a base64 copy of itself in the notebook.
What is accurate here, and what is not¶
Every orbit is drawn from published orbital elements — altitude and inclination per satellite — with the period derived from the radius by Kepler's third law rather than assumed. So the relative sizes of the orbits and the tilt of each orbital plane are physically right, and every period is a real period rather than a chosen one.
Radii are compressed, and this is the one place the drawing departs from the numbers. At true scale the
entire low-orbit fleet sits in a band 1% of the frame wide, pinned to Earth's limb, while geostationary is 6.6
radii out — accurate, and unreadable. Altitudes above the surface are therefore raised to a power (RADIAL_GAMMA,
0.34), which lifts the low-orbit fleet clear of the surface and pulls the two populations together, while
preserving their order and the fact that a gap exists:
| true radius | drawn at | |
|---|---|---|
| GPM Core (lowest) | 1.064 | 1.48 |
| Sentinel-1 | 1.109 | 1.58 |
| NOAA-20 | 1.129 | 1.61 |
| Sentinel-6 | 1.210 | 1.72 |
| geostationary | 6.617 | 3.20 |
Read the radial axis as ordering, not as distance: geostationary really is about six times further out than the low-orbit shell, not twice. Radii are the single presentational liberty; every rate, period and inclination is the real one.
That liberty reaches one derived thing as well. Each sensor beam hangs from its satellite at the drawn radius while its ground edge comes from the true swath, so the beam lands on the right strip but opens at the wrong angle -- VIIRS draws about 20° from the spacecraft where the real half-angle is 56°. It is a footprint indicator, not a measured field of view. Widening the beams to their true angle at the compressed radius would have them overlap into a solid sheet across the shell, which is the picture the compression exists to avoid.
Rates are not touched at all, and the window is short instead. A low orbiter at its true rate laps the frame about fifteen times a day, which is far too fast to follow. The obvious fix -- draw the low orbits at a fraction of their real speed -- is the wrong one, because it silently breaks the ground tracks below. A ground track is the product of two motions, the satellite going round and Earth turning underneath, and its shape encodes the ratio between them. Slow one side only and every swath drifts west four times too fast.
So the clip covers 1.5 hours rather than six, and everything inside it runs at its true rate. That yields the same on-screen speed -- 0.25 x 6.0 h and 1.0 x 1.5 h both put a low orbiter a little under one circuit around its ring -- while keeping the physics self-consistent. The cost is that Earth turns 22.6 degrees across the clip rather than 90, so the geostationary ring sweeps a sixteenth of a turn instead of a quarter. It still hangs over the same patch of ground, which is the part worth seeing.
What is not modelled: orbital precession and drift, the equation of time, atmospheric drag, eccentricity (every orbit here is drawn circular, which for these missions is close), and the true right ascension of each ascending node. Where a satellite's plane orientation is not pinned by its mission, the notebook spreads the planes evenly so the shell reads clearly instead of collapsing into one line — the inclinations are real, the phasing between planes is presentational.
Satellite bodies are drawn as glyphs, not to scale. At true scale a satellite would be far smaller than one pixel.
The strip each satellite actually sees¶
A satellite in low orbit does not photograph the planet. It sweeps a strip, and only that strip, then relies on Earth turning beneath it to bring a new one round on the next pass. The trailing bands show that directly: each one is the ground its satellite has crossed over the last quarter of an hour, fading out behind it, painted onto the surface and carried round with the rotation rather than sliding along with the spacecraft.
The widths are the real instrument swaths, and they differ by more than two orders of magnitude. VIIRS on Suomi-NPP and the two NOAA satellites cuts a band about 3,000 km across, wide enough to cover the planet in a day. Landsat's imager manages 185 km. OCO-2 and ICESat-2 draw lines thinner than the coastlines beneath them -- which is not a rendering failure but the honest answer to how much of Earth a spectrometer or a laser altimeter sees at once.
In [ ]:
Copied!
import hashlib
import os
import subprocess
import tempfile
import warnings
from pathlib import Path
import imageio_ffmpeg
import matplotlib
import matplotlib.animation as manim
import matplotlib.patheffects as pe
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import Image, display
from matplotlib.colors import to_rgba
from matplotlib.path import Path as MPath # the glyph path; `Path` is pathlib's
from mpl_toolkits.mplot3d.art3d import Line3DCollection, Poly3DCollection
from scipy.ndimage import gaussian_filter
# Two known-noisy sources, filtered by name rather than by silencing the
# kernel: matplotlib re-warns about the 3D axes on every frame, and the tile
# fetch emits an urllib3 chunking warning per connection. Anything else --
# a degenerate normalisation, a numpy invalid value -- should still be heard.
warnings.filterwarnings("ignore", category=DeprecationWarning, module="matplotlib")
warnings.filterwarnings("ignore", category=UserWarning, module="matplotlib")
# render_master writes the MP4 through ffmpeg; imageio_ffmpeg ships the binary
# so nothing has to be installed system-wide.
matplotlib.rcParams["animation.ffmpeg_path"] = imageio_ffmpeg.get_ffmpeg_exe()
# Earth is the unit: every distance below is in Earth radii. Radii are
# compressed for drawing (see draw_radius); the physics uses the true values.
R_EARTH_KM = 6371.0
MU = 398600.4418 # Earth's gravitational parameter, km^3/s^2
# The drawing works in mean radii, but the 35,786 km geostationary altitude is
# defined against the equatorial radius. The check below only lands on the
# sidereal day when it is computed with the radius the altitude refers to --
# with the mean radius it comes out 22 s short, which is the difference between
# the two radii, not anything about the orbit.
R_EQUATORIAL_KM = 6378.137
TILT = np.deg2rad(23.44) # axial tilt
# The globe is lit evenly rather than from a direction: no Sun, no terminator,
# no night side. The gain that multiplies the texture is not fixed here -- it
# belongs to the basemap, since imagery that is dark by design needs lifting and
# imagery that is already bright does not. See BASEMAPS below.
#: Textures already built this session, keyed by every setting that shapes
#: one. A bare single-slot memo would hand back the previous basemap's
#: pixels after BASEMAP is changed and the cells below are re-run, which
#: publishes a file named for one basemap containing another's imagery.
_BASEMAP: dict[tuple, np.ndarray] = {}
#: The spin-independent half of the globe, cached the same way.
_EARTH_FACE: dict[tuple, tuple] = {}
import hashlib
import os
import subprocess
import tempfile
import warnings
from pathlib import Path
import imageio_ffmpeg
import matplotlib
import matplotlib.animation as manim
import matplotlib.patheffects as pe
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import Image, display
from matplotlib.colors import to_rgba
from matplotlib.path import Path as MPath # the glyph path; `Path` is pathlib's
from mpl_toolkits.mplot3d.art3d import Line3DCollection, Poly3DCollection
from scipy.ndimage import gaussian_filter
# Two known-noisy sources, filtered by name rather than by silencing the
# kernel: matplotlib re-warns about the 3D axes on every frame, and the tile
# fetch emits an urllib3 chunking warning per connection. Anything else --
# a degenerate normalisation, a numpy invalid value -- should still be heard.
warnings.filterwarnings("ignore", category=DeprecationWarning, module="matplotlib")
warnings.filterwarnings("ignore", category=UserWarning, module="matplotlib")
# render_master writes the MP4 through ffmpeg; imageio_ffmpeg ships the binary
# so nothing has to be installed system-wide.
matplotlib.rcParams["animation.ffmpeg_path"] = imageio_ffmpeg.get_ffmpeg_exe()
# Earth is the unit: every distance below is in Earth radii. Radii are
# compressed for drawing (see draw_radius); the physics uses the true values.
R_EARTH_KM = 6371.0
MU = 398600.4418 # Earth's gravitational parameter, km^3/s^2
# The drawing works in mean radii, but the 35,786 km geostationary altitude is
# defined against the equatorial radius. The check below only lands on the
# sidereal day when it is computed with the radius the altitude refers to --
# with the mean radius it comes out 22 s short, which is the difference between
# the two radii, not anything about the orbit.
R_EQUATORIAL_KM = 6378.137
TILT = np.deg2rad(23.44) # axial tilt
# The globe is lit evenly rather than from a direction: no Sun, no terminator,
# no night side. The gain that multiplies the texture is not fixed here -- it
# belongs to the basemap, since imagery that is dark by design needs lifting and
# imagery that is already bright does not. See BASEMAPS below.
#: Textures already built this session, keyed by every setting that shapes
#: one. A bare single-slot memo would hand back the previous basemap's
#: pixels after BASEMAP is changed and the cells below are re-run, which
#: publishes a file named for one basemap containing another's imagery.
_BASEMAP: dict[tuple, np.ndarray] = {}
#: The spin-independent half of the globe, cached the same way.
_EARTH_FACE: dict[tuple, tuple] = {}
In [ ]:
Copied!
def tilt_y(x, z):
"""Apply Earth's 23.44° axial tilt, a rotation about the y-axis.
Args:
x: The x coordinate, or an array of them.
z: The z coordinate, or an array of them, matching `x`.
Returns:
tuple: The rotated `(x, z)`; `y` is the rotation axis and is unchanged.
"""
return x * np.cos(TILT) + z * np.sin(TILT), -x * np.sin(TILT) + z * np.cos(TILT)
def sphere(c, r, n=40):
"""Surface mesh of a sphere.
Args:
c: Centre as an `(x, y, z)` sequence.
r: Radius, in the same units as `c`.
n: Samples per angular axis; the mesh is `n` by `n`.
Returns:
tuple: The `(x, y, z)` meshes, each of shape `(n, n)`.
"""
u = np.linspace(0, 2 * np.pi, n)
v = np.linspace(0, np.pi, n)
return (
c[0] + r * np.outer(np.cos(u), np.sin(v)),
c[1] + r * np.outer(np.sin(u), np.sin(v)),
c[2] + r * np.outer(np.ones_like(u), np.cos(v)),
)
def ll_to_xyz(lon, lat, spin):
"""Longitude and latitude to a tilted, spun point on Earth's surface.
Args:
lon: Longitude in degrees, or an array of them.
lat: Latitude in degrees, matching `lon`.
spin: Earth's rotation angle at this moment, in radians.
Returns:
tuple: The `(x, y, z)` coordinates in Earth radii, on the unit sphere.
"""
lon = np.deg2rad(lon) + spin
lat = np.deg2rad(lat)
x = np.cos(lat) * np.cos(lon)
y = np.cos(lat) * np.sin(lon)
z = np.sin(lat)
xr, zr = tilt_y(x, z)
return xr, y, zr
#: Every basemap this notebook publishes a clip from, and the handling each
#: one's imagery needs. To render a different one, change BASEMAP below and run
#: the notebook again -- nothing else has to be touched.
#:
#: - `two_tone` flattens the basemap onto OCEAN_RGB / LAND_RGB. Right for a
#: cartographic source like DarkMatter, which carries no imagery worth
#: keeping; wrong for real imagery, where the picture itself is the point.
#: - `black_point` subtracts the imagery's own colour floor. NASA's night lights
#: paint the unlit surface deep blue rather than black, and a gain then
#: amplifies that into a blue planet; subtracting the per-channel floor puts
#: the background at true black and barely touches the city lights.
#: - `gain` multiplies the texture, above 1 for imagery that is dark by design.
BASEMAPS = {
"night": {
"provider": "NASAGIBS.ViirsEarthAtNight2012",
"two_tone": False,
"black_point": True,
"gain": 2.3,
},
"bluemarble": {
"provider": "NASAGIBS.BlueMarble",
"two_tone": False,
"black_point": False,
"gain": 1.0,
},
"darkmatter": {
"provider": "CartoDB.DarkMatterNoLabels",
"two_tone": True,
"black_point": False,
"gain": 1.0,
},
}
#: The basemaps this notebook publishes a clip from. The render cell walks this
#: tuple, so one top-to-bottom run produces the whole published set from one
#: state of the code -- rather than refreshing whichever half BASEMAP happened
#: to name and leaving the other stale against it.
PUBLISHED_BASEMAPS = ("night", "bluemarble")
def select_basemap(name):
"""Point every basemap-derived constant at one entry in `BASEMAPS`.
The constants below are read all over the notebook, so switching basemap
means reassigning them together. Doing it here keeps the render loop from
setting five globals by hand and getting one of them wrong.
Args:
name: A key of `BASEMAPS`.
Returns:
dict: That basemap's entry.
"""
global BASEMAP, BASEMAP_PROVIDER, BASEMAP_TWO_TONE, BASEMAP_BLACK_POINT
global BRIGHTNESS
spec = BASEMAPS[name]
BASEMAP = name
BASEMAP_PROVIDER = spec["provider"]
BASEMAP_TWO_TONE = spec["two_tone"]
BASEMAP_BLACK_POINT = spec["black_point"]
BRIGHTNESS = spec["gain"]
return spec
#: Which basemap the preview frame further down is drawn with. The render cell
#: overrides this per clip, so it only affects what you see while working.
BASEMAP = "night"
select_basemap(BASEMAP)
#: Zoom of the tile pyramid to fetch: the whole world is `4 ** zoom` tiles, so 5
#: is 1024. Ample for the texture below, where 6 would be 4096 for no visible gain.
BASEMAP_ZOOM = 5
#: Web Mercator's half-extent in metres and the latitude it stops at. The
#: projection never reaches the poles, so the caps repeat the last real row.
MERC_MAX = 20037508.342789244
MERC_LAT_MAX = 85.05112877980659
#: The two tones a `two_tone` basemap is reduced to, and the luminance below
#: which a pixel counts as land. DarkMatter draws land darker than water, so the
#: test reads inverted from what it looks like it should be.
OCEAN_RGB = np.array([0.045, 0.062, 0.095])
LAND_RGB = np.array([0.235, 0.255, 0.225])
LAND_LUMA_CUT = 0.09
#: Percentile of each channel taken as the imagery's own colour floor when
#: `black_point` is on. Low enough to sit in the unlit background rather
#: than in anything lit.
BLACK_POINT_PCT = 2.0
def basemap_texture(n_lon=2880, n_lat=1440):
"""The basemap as an `(n_lat, n_lon, 3)` equirectangular array in 0..1.
XYZ tiles arrive in Web Mercator, but a sphere is textured in longitude and
latitude, so the stitched mosaic is resampled onto a lon/lat grid here. The
result is cached under the earthlens cache directory, so the tiles are
fetched once per provider rather than on every run.
Args:
n_lon: Width of the returned texture, spanning -180..180 degrees.
n_lat: Height of the returned texture, spanning 90..-90 degrees.
Returns:
numpy.ndarray: The `(n_lat, n_lon, 3)` texture, values in 0..1.
"""
# Everything that changes a pixel, in one tuple: it is both the in-session
# memo key and the on-disk cache key. Tuning a constant that is not in here
# would silently serve the texture built before the change -- the failure is
# invisible, and the honest alternative costs a refetch of a thousand tiles.
settings = (
BASEMAP,
BASEMAP_PROVIDER,
BASEMAP_ZOOM,
n_lon,
n_lat,
BASEMAP_TWO_TONE,
BASEMAP_BLACK_POINT,
BLACK_POINT_PCT,
tuple(OCEAN_RGB),
tuple(LAND_RGB),
LAND_LUMA_CUT,
)
if settings in _BASEMAP:
return _BASEMAP[settings]
from earthlens.core import cache_dir
cache = Path(cache_dir()) / "basemap-textures"
cache.mkdir(parents=True, exist_ok=True)
digest = hashlib.blake2s(repr(settings).encode(), digest_size=8).hexdigest()
key = cache / f"{BASEMAP_PROVIDER.replace('.', '_')}_z{BASEMAP_ZOOM}_{digest}.npy"
if key.exists():
_BASEMAP[settings] = np.load(key)
return _BASEMAP[settings]
from cleopatra.basemap.tiles import Tile, fetch_tiles, get_provider, stitch_tiles
provider = get_provider(BASEMAP_PROVIDER)
side = 2**BASEMAP_ZOOM
tiles = [Tile(x, y, BASEMAP_ZOOM) for x in range(side) for y in range(side)]
mosaic, (west, south, east, north) = stitch_tiles(
fetch_tiles(tiles, provider, max_workers=8, timeout=20, retries=3),
tiles,
BASEMAP_ZOOM,
)
mosaic = np.asarray(mosaic)[..., :3]
h, w = mosaic.shape[:2]
# Every output cell averages the whole contiguous block of source pixels that
# falls inside it, found from the cell's *edges* rather than its centre.
#
# Taking a fixed number of point samples per cell instead was inverted where
# it mattered: Mercator stretches towards the poles, so three samples landed
# on about one source row at the equator (an average of nothing) and skipped
# seven rows out of every eight above 80 degrees, which is precisely where
# the stretch makes aliasing worst. Half the fetched mosaic was never read.
#
# The mosaic's own extent is used, not an assumed full-world one. They agree
# while the tile list is the complete grid -- but a reprojection should read
# its bounds from the thing it is reprojecting.
lat_edges = np.clip(
np.linspace(90.0, -90.0, n_lat + 1), -MERC_LAT_MAX, MERC_LAT_MAX
)
y_edges = np.log(np.tan(np.pi / 4 + np.deg2rad(lat_edges) / 2)) * MERC_MAX / np.pi
row_edges = np.clip(
np.round((north - y_edges) / (north - south) * h).astype(int), 0, h
)
x_edges = np.linspace(-180.0, 180.0, n_lon + 1) / 180.0 * MERC_MAX
col_edges = np.clip(
np.round((x_edges - west) / (east - west) * w).astype(int), 0, w
)
# Where a cell is thinner than one source pixel -- the polar rows, which the
# Mercator clamp collapses onto a single row -- reduceat returns that row on
# its own, which is the right fallback. Dividing by max(width, 1) keeps it.
row_start = np.minimum(row_edges[:-1], h - 1)
col_start = np.minimum(col_edges[:-1], w - 1)
# float32 divisors: dividing a float32 band by integer counts would promote
# the result to float64, doubling both the resident texture and the cached
# file for precision a 0..1 colour has no use for.
row_w = np.maximum(np.diff(row_edges), 1)[:, None].astype("float32")
col_w = np.maximum(np.diff(col_edges), 1)[None, :].astype("float32")
# One band at a time: the float32 copy of an 8192-square mosaic is 256 MB per
# band, and converting all three at once put three quarters of a gigabyte
# live at the peak for a 47 MB result.
bands = []
for b in range(3):
band = mosaic[..., b].astype("float32")
if band.max() > 1.5:
band /= 255.0
band = np.add.reduceat(band, row_start, axis=0) / row_w
bands.append(np.add.reduceat(band, col_start, axis=1) / col_w)
grid = np.stack(bands, axis=-1)
# DarkMatter paints land near-black and water a lighter grey, the opposite of
# how a picture of Earth reads. The tiles are effectively two tones, so they
# are remapped onto a chosen pair: a deep ocean and a lighter land. The water
# stays dark enough that a swath laid over it is the brightest thing in its
# area, while the continents stay legible.
if BASEMAP_TWO_TONE:
lum = grid @ np.array([0.2126, 0.7152, 0.0722])
out = np.empty_like(grid)
out[...] = OCEAN_RGB
out[lum < LAND_LUMA_CUT] = LAND_RGB
else:
out = grid
if BASEMAP_BLACK_POINT:
floor = np.percentile(out.reshape(-1, 3), BLACK_POINT_PCT, axis=0)
out = out - floor
texture = np.clip(out, 0, 1)
# Written aside and renamed, because np.save straight onto the final path
# leaves a truncated array at the exact name the next run trusts if the
# write is interrupted -- and the remedy is knowing to delete a file buried
# in the platform cache directory.
tmp = key.with_name(key.name + ".tmp")
with open(tmp, "wb") as fh:
np.save(fh, texture)
os.replace(tmp, key)
_BASEMAP[settings] = texture
return texture
def sample_texture(tex, lon_g, lat_g):
"""Bilinear sample of an equirectangular texture on a longitude/latitude grid.
Nearest-neighbour sampling reads one texel per mesh cell, so at any mesh coarser
than the texture most of the texture is simply skipped and coastlines break into
steps. Interpolating between the four surrounding texels uses the detail that is
there, which is what makes a finer mesh worth paying for.
Args:
tex: `(h, w, 3)` texture in 0..1, spanning -180..180 longitude and 90..-90 latitude.
lon_g: Longitudes to sample, in degrees.
lat_g: Latitudes to sample, in degrees, matching the shape of `lon_g`.
Returns:
An `(..., 3)` array of sampled colours in 0..1.
"""
h, w = tex.shape[:2]
fx = (lon_g + 180.0) / 360.0 * (w - 1)
fy = (90.0 - lat_g) / 180.0 * (h - 1)
x0 = np.clip(np.floor(fx).astype(int), 0, w - 1)
y0 = np.clip(np.floor(fy).astype(int), 0, h - 1)
x1, y1 = np.clip(x0 + 1, 0, w - 1), np.clip(y0 + 1, 0, h - 1)
tx, ty = (fx - x0)[..., None], (fy - y0)[..., None]
top = tex[y0, x0] * (1 - tx) + tex[y0, x1] * tx
bot = tex[y1, x0] * (1 - tx) + tex[y1, x1] * tx
return top * (1 - ty) + bot * ty
#: How finely the globe is meshed. This is the render's cost driver and it is
#: quadratic in the numbers below, so it dominates everything else: at 1920 px
#: wide a 1440x720 mesh works out at roughly one screen pixel per cell, which is
#: the point past which a finer mesh cannot be seen. It must be defined before
#: `textured_earth`, which takes it as a default argument.
MESH_LON, MESH_LAT = 1440, 720
def _earth_face(n_lon, n_lat):
"""The lon/lat mesh and its lit colours, which are the same in every frame.
Neither the mesh nor the sampled texture depends on `spin` -- rotating the
globe moves where those colours are drawn, not what they are. Sampling is a
million bilinear gathers, so doing it inside the draw loop repeated the
identical work on all 480 frames.
Args:
n_lon: Mesh columns, spanning -180..180 degrees.
n_lat: Mesh rows, spanning 90..-90 degrees.
Returns:
tuple: `(lon_grid, lat_grid, face_colours)`, the last already lit and
clipped to 0..1.
"""
key = (n_lon, n_lat, BASEMAP, BRIGHTNESS)
if key not in _EARTH_FACE:
lon_g, lat_g = np.meshgrid(
np.linspace(-180.0, 180.0, n_lon), np.linspace(90.0, -90.0, n_lat)
)
rgb = sample_texture(basemap_texture(), lon_g, lat_g)
# Evenly lit: no light direction, so no terminator and no night side.
# Every longitude reads the same, which keeps a satellite passing over
# the far limb as legible as one crossing the middle.
_EARTH_FACE[key] = (lon_g, lat_g, np.clip(rgb * BRIGHTNESS, 0, 1))
return _EARTH_FACE[key]
def textured_earth(ax, spin, n_lon=MESH_LON, n_lat=MESH_LAT):
"""Earth as an evenly lit, texture-mapped sphere -- real continents, tilted and spun.
The mesh is sampled at 1440 x 720 against a 2880 x 1440 basemap. A coarser mesh
was the source of the blockiness: the texture was always this detailed, the
sphere just was not sampling it.
Args:
ax: The 3D axis to draw onto.
spin: Earth's rotation angle for this frame, in radians.
n_lon: Mesh columns spanning -180..180 degrees.
n_lat: Mesh rows spanning 90..-90 degrees.
"""
lon_g, lat_g, face = _earth_face(n_lon, n_lat)
lo, la = np.deg2rad(lon_g) + spin, np.deg2rad(lat_g)
x, y, z = np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)
xr, zr = tilt_y(x, z)
ax.plot_surface(
xr,
y,
zr,
facecolors=face,
rstride=1,
cstride=1,
linewidth=0,
antialiased=False,
shade=False,
zorder=Z_GLOBE,
)
def render_master(fig, render, n_frames, stem, fps=24, bitrate=12000, dpi=160):
"""Draw every frame once, into the master mp4.
This is the only expensive step: `render` is called `n_frames` times and each call
draws a full frame. Every published size is then derived from the file this writes,
by `publish_clip`, which transcodes rather than redrawing.
Args:
fig: The figure `render` draws onto.
render: Callable taking a frame index and drawing that frame onto `fig`.
n_frames: How many frames to draw.
stem: Output path without a suffix; `.mp4` is written beside it.
fps: Frame rate of the master.
bitrate: Encoder bitrate in kbit/s. An mp4's size is bitrate times duration
regardless of resolution, so this has to rise with the frame size or the
detail the render spent hours on is simply compressed away.
dpi: Dots per inch, which together with the figure size sets the pixel
dimensions of the master.
Returns:
str: Path to the mp4 that was written.
"""
writer = manim.FFMpegWriter(fps=fps, bitrate=bitrate)
mp4 = f"{stem}.mp4"
with writer.saving(fig, mp4, dpi):
for i in range(n_frames):
render(i)
writer.grab_frame()
return mp4
def publish_clip(
master, web_width=1280, gif_width=880, gif_fps=12, webp_width=1080, webp_fps=15
):
"""Derive every published size from the rendered master.
Drawing the frames is the expensive part and it happens once, at the largest size
any target needs; every other file is a transcode of those same frames.
Three outputs, because the three places this clip appears do not accept the same
thing. The documentation site takes an mp4, which it can embed in a `<video>`.
GitHub strips video tags from README markdown, so the README needs an image --
and an animated webp beats a gif on every axis here, carrying truecolour at a
wider size for fewer bytes. The gif remains as the conservative fallback.
The gif settings are not defaults. A 256-colour palette measured 1.8 dB better
than 128 on this material, `stats_mode=full` suits a composition that is mostly
static, and dithering is off because error diffusion speckles a dark, smooth
frame -- measured better *and* smaller than the diffused version. 880 px matches
GitHub's README column, so nothing is scaled in either direction.
Args:
master: Path to the rendered master mp4.
web_width: Width of the mp4 the documentation site embeds.
gif_width: Width of the gif. Below about 700 px the satellite labels stop
being legible, so this is a floor rather than a preference.
gif_fps: Frame rate of the gif; halving it halves the file with no visible
judder, since the clip's duration is unchanged.
webp_width: Width of the animated webp used in the README.
webp_fps: Frame rate of the webp, which can afford more than the gif.
Returns:
tuple: The `(web_mp4, gif, webp)` paths that were written.
"""
master = Path(master)
exe = imageio_ffmpeg.get_ffmpeg_exe()
stem = str(master.with_suffix("")).rsplit("-1920", 1)[0]
web, gif, webp = f"{stem}-{web_width}.mp4", f"{stem}.gif", f"{stem}.webp"
# The palette is scratch, so it goes to the system temp directory rather
# than beside the published clips: the output directory is tracked, and a
# failed transcode used to strand a stray PNG in it that nothing ignores.
palette = os.path.join(tempfile.gettempdir(), f"{master.stem}_palette.png")
subprocess.run(
[
exe,
"-y",
"-v",
"error",
"-i",
str(master),
"-vf",
f"scale={web_width}:-2:flags=lanczos",
"-c:v",
"libx264",
"-b:v",
"6000k",
"-pix_fmt",
"yuv420p",
web,
],
check=True,
)
subprocess.run(
[
exe,
"-y",
"-v",
"error",
"-i",
str(master),
"-vf",
f"fps={webp_fps},scale={webp_width}:-1:flags=lanczos",
"-c:v",
"libwebp",
"-lossless",
"0",
"-q:v",
"72",
"-loop",
"0",
"-an",
webp,
],
check=True,
)
scale = f"fps={gif_fps},scale={gif_width}:-1:flags=lanczos"
try:
subprocess.run(
[
exe,
"-y",
"-v",
"error",
"-i",
str(master),
"-vf",
f"{scale},palettegen=max_colors=256:stats_mode=full",
palette,
],
check=True,
)
subprocess.run(
[
exe,
"-y",
"-v",
"error",
"-i",
str(master),
"-i",
palette,
"-lavfi",
f"{scale} [x]; [x][1:v] paletteuse=dither=none:diff_mode=rectangle",
gif,
],
check=True,
)
finally:
if os.path.exists(palette):
os.remove(palette)
return web, gif, webp
def _glyph_rect(x0, y0, x1, y1):
"""Vertices and codes for one closed rectangle of a compound marker path.
Args:
x0: Left edge, in marker units.
y0: Bottom edge, in marker units.
x1: Right edge, in marker units.
y1: Top edge, in marker units.
Returns:
tuple: The `(vertices, codes)` lists to concatenate into a `Path`.
"""
return (
[(x0, y0), (x1, y0), (x1, y1), (x0, y1), (x0, y0)],
[MPath.MOVETO, MPath.LINETO, MPath.LINETO, MPath.LINETO, MPath.CLOSEPOLY],
)
def _satellite_marker():
"""A satellite silhouette: two solar wings flanking a small bus on short booms.
Drawn as a marker rather than as 3D geometry on purpose. A modelled satellite would
be two or three pixels across in the wide framing and vanish; a marker is sized in
points, so it stays legible at every zoom.
Returns:
matplotlib.path.Path: The compound path, centred on the origin and
scaled to roughly one marker unit across.
"""
verts, codes = [], []
for box in (
(-0.17, -0.30, 0.17, 0.30), # bus
(-0.40, -0.035, -0.17, 0.035), # left boom
(0.17, -0.035, 0.40, 0.035), # right boom
(-1.00, -0.62, -0.40, 0.62), # left wing
(0.40, -0.62, 1.00, 0.62),
): # right wing
v, c = _glyph_rect(*box)
verts += v
codes += c
return MPath(verts, codes)
#: The glyph every satellite is drawn with.
SATELLITE = _satellite_marker()
def tilt_y(x, z):
"""Apply Earth's 23.44° axial tilt, a rotation about the y-axis.
Args:
x: The x coordinate, or an array of them.
z: The z coordinate, or an array of them, matching `x`.
Returns:
tuple: The rotated `(x, z)`; `y` is the rotation axis and is unchanged.
"""
return x * np.cos(TILT) + z * np.sin(TILT), -x * np.sin(TILT) + z * np.cos(TILT)
def sphere(c, r, n=40):
"""Surface mesh of a sphere.
Args:
c: Centre as an `(x, y, z)` sequence.
r: Radius, in the same units as `c`.
n: Samples per angular axis; the mesh is `n` by `n`.
Returns:
tuple: The `(x, y, z)` meshes, each of shape `(n, n)`.
"""
u = np.linspace(0, 2 * np.pi, n)
v = np.linspace(0, np.pi, n)
return (
c[0] + r * np.outer(np.cos(u), np.sin(v)),
c[1] + r * np.outer(np.sin(u), np.sin(v)),
c[2] + r * np.outer(np.ones_like(u), np.cos(v)),
)
def ll_to_xyz(lon, lat, spin):
"""Longitude and latitude to a tilted, spun point on Earth's surface.
Args:
lon: Longitude in degrees, or an array of them.
lat: Latitude in degrees, matching `lon`.
spin: Earth's rotation angle at this moment, in radians.
Returns:
tuple: The `(x, y, z)` coordinates in Earth radii, on the unit sphere.
"""
lon = np.deg2rad(lon) + spin
lat = np.deg2rad(lat)
x = np.cos(lat) * np.cos(lon)
y = np.cos(lat) * np.sin(lon)
z = np.sin(lat)
xr, zr = tilt_y(x, z)
return xr, y, zr
#: Every basemap this notebook publishes a clip from, and the handling each
#: one's imagery needs. To render a different one, change BASEMAP below and run
#: the notebook again -- nothing else has to be touched.
#:
#: - `two_tone` flattens the basemap onto OCEAN_RGB / LAND_RGB. Right for a
#: cartographic source like DarkMatter, which carries no imagery worth
#: keeping; wrong for real imagery, where the picture itself is the point.
#: - `black_point` subtracts the imagery's own colour floor. NASA's night lights
#: paint the unlit surface deep blue rather than black, and a gain then
#: amplifies that into a blue planet; subtracting the per-channel floor puts
#: the background at true black and barely touches the city lights.
#: - `gain` multiplies the texture, above 1 for imagery that is dark by design.
BASEMAPS = {
"night": {
"provider": "NASAGIBS.ViirsEarthAtNight2012",
"two_tone": False,
"black_point": True,
"gain": 2.3,
},
"bluemarble": {
"provider": "NASAGIBS.BlueMarble",
"two_tone": False,
"black_point": False,
"gain": 1.0,
},
"darkmatter": {
"provider": "CartoDB.DarkMatterNoLabels",
"two_tone": True,
"black_point": False,
"gain": 1.0,
},
}
#: The basemaps this notebook publishes a clip from. The render cell walks this
#: tuple, so one top-to-bottom run produces the whole published set from one
#: state of the code -- rather than refreshing whichever half BASEMAP happened
#: to name and leaving the other stale against it.
PUBLISHED_BASEMAPS = ("night", "bluemarble")
def select_basemap(name):
"""Point every basemap-derived constant at one entry in `BASEMAPS`.
The constants below are read all over the notebook, so switching basemap
means reassigning them together. Doing it here keeps the render loop from
setting five globals by hand and getting one of them wrong.
Args:
name: A key of `BASEMAPS`.
Returns:
dict: That basemap's entry.
"""
global BASEMAP, BASEMAP_PROVIDER, BASEMAP_TWO_TONE, BASEMAP_BLACK_POINT
global BRIGHTNESS
spec = BASEMAPS[name]
BASEMAP = name
BASEMAP_PROVIDER = spec["provider"]
BASEMAP_TWO_TONE = spec["two_tone"]
BASEMAP_BLACK_POINT = spec["black_point"]
BRIGHTNESS = spec["gain"]
return spec
#: Which basemap the preview frame further down is drawn with. The render cell
#: overrides this per clip, so it only affects what you see while working.
BASEMAP = "night"
select_basemap(BASEMAP)
#: Zoom of the tile pyramid to fetch: the whole world is `4 ** zoom` tiles, so 5
#: is 1024. Ample for the texture below, where 6 would be 4096 for no visible gain.
BASEMAP_ZOOM = 5
#: Web Mercator's half-extent in metres and the latitude it stops at. The
#: projection never reaches the poles, so the caps repeat the last real row.
MERC_MAX = 20037508.342789244
MERC_LAT_MAX = 85.05112877980659
#: The two tones a `two_tone` basemap is reduced to, and the luminance below
#: which a pixel counts as land. DarkMatter draws land darker than water, so the
#: test reads inverted from what it looks like it should be.
OCEAN_RGB = np.array([0.045, 0.062, 0.095])
LAND_RGB = np.array([0.235, 0.255, 0.225])
LAND_LUMA_CUT = 0.09
#: Percentile of each channel taken as the imagery's own colour floor when
#: `black_point` is on. Low enough to sit in the unlit background rather
#: than in anything lit.
BLACK_POINT_PCT = 2.0
def basemap_texture(n_lon=2880, n_lat=1440):
"""The basemap as an `(n_lat, n_lon, 3)` equirectangular array in 0..1.
XYZ tiles arrive in Web Mercator, but a sphere is textured in longitude and
latitude, so the stitched mosaic is resampled onto a lon/lat grid here. The
result is cached under the earthlens cache directory, so the tiles are
fetched once per provider rather than on every run.
Args:
n_lon: Width of the returned texture, spanning -180..180 degrees.
n_lat: Height of the returned texture, spanning 90..-90 degrees.
Returns:
numpy.ndarray: The `(n_lat, n_lon, 3)` texture, values in 0..1.
"""
# Everything that changes a pixel, in one tuple: it is both the in-session
# memo key and the on-disk cache key. Tuning a constant that is not in here
# would silently serve the texture built before the change -- the failure is
# invisible, and the honest alternative costs a refetch of a thousand tiles.
settings = (
BASEMAP,
BASEMAP_PROVIDER,
BASEMAP_ZOOM,
n_lon,
n_lat,
BASEMAP_TWO_TONE,
BASEMAP_BLACK_POINT,
BLACK_POINT_PCT,
tuple(OCEAN_RGB),
tuple(LAND_RGB),
LAND_LUMA_CUT,
)
if settings in _BASEMAP:
return _BASEMAP[settings]
from earthlens.core import cache_dir
cache = Path(cache_dir()) / "basemap-textures"
cache.mkdir(parents=True, exist_ok=True)
digest = hashlib.blake2s(repr(settings).encode(), digest_size=8).hexdigest()
key = cache / f"{BASEMAP_PROVIDER.replace('.', '_')}_z{BASEMAP_ZOOM}_{digest}.npy"
if key.exists():
_BASEMAP[settings] = np.load(key)
return _BASEMAP[settings]
from cleopatra.basemap.tiles import Tile, fetch_tiles, get_provider, stitch_tiles
provider = get_provider(BASEMAP_PROVIDER)
side = 2**BASEMAP_ZOOM
tiles = [Tile(x, y, BASEMAP_ZOOM) for x in range(side) for y in range(side)]
mosaic, (west, south, east, north) = stitch_tiles(
fetch_tiles(tiles, provider, max_workers=8, timeout=20, retries=3),
tiles,
BASEMAP_ZOOM,
)
mosaic = np.asarray(mosaic)[..., :3]
h, w = mosaic.shape[:2]
# Every output cell averages the whole contiguous block of source pixels that
# falls inside it, found from the cell's *edges* rather than its centre.
#
# Taking a fixed number of point samples per cell instead was inverted where
# it mattered: Mercator stretches towards the poles, so three samples landed
# on about one source row at the equator (an average of nothing) and skipped
# seven rows out of every eight above 80 degrees, which is precisely where
# the stretch makes aliasing worst. Half the fetched mosaic was never read.
#
# The mosaic's own extent is used, not an assumed full-world one. They agree
# while the tile list is the complete grid -- but a reprojection should read
# its bounds from the thing it is reprojecting.
lat_edges = np.clip(
np.linspace(90.0, -90.0, n_lat + 1), -MERC_LAT_MAX, MERC_LAT_MAX
)
y_edges = np.log(np.tan(np.pi / 4 + np.deg2rad(lat_edges) / 2)) * MERC_MAX / np.pi
row_edges = np.clip(
np.round((north - y_edges) / (north - south) * h).astype(int), 0, h
)
x_edges = np.linspace(-180.0, 180.0, n_lon + 1) / 180.0 * MERC_MAX
col_edges = np.clip(
np.round((x_edges - west) / (east - west) * w).astype(int), 0, w
)
# Where a cell is thinner than one source pixel -- the polar rows, which the
# Mercator clamp collapses onto a single row -- reduceat returns that row on
# its own, which is the right fallback. Dividing by max(width, 1) keeps it.
row_start = np.minimum(row_edges[:-1], h - 1)
col_start = np.minimum(col_edges[:-1], w - 1)
# float32 divisors: dividing a float32 band by integer counts would promote
# the result to float64, doubling both the resident texture and the cached
# file for precision a 0..1 colour has no use for.
row_w = np.maximum(np.diff(row_edges), 1)[:, None].astype("float32")
col_w = np.maximum(np.diff(col_edges), 1)[None, :].astype("float32")
# One band at a time: the float32 copy of an 8192-square mosaic is 256 MB per
# band, and converting all three at once put three quarters of a gigabyte
# live at the peak for a 47 MB result.
bands = []
for b in range(3):
band = mosaic[..., b].astype("float32")
if band.max() > 1.5:
band /= 255.0
band = np.add.reduceat(band, row_start, axis=0) / row_w
bands.append(np.add.reduceat(band, col_start, axis=1) / col_w)
grid = np.stack(bands, axis=-1)
# DarkMatter paints land near-black and water a lighter grey, the opposite of
# how a picture of Earth reads. The tiles are effectively two tones, so they
# are remapped onto a chosen pair: a deep ocean and a lighter land. The water
# stays dark enough that a swath laid over it is the brightest thing in its
# area, while the continents stay legible.
if BASEMAP_TWO_TONE:
lum = grid @ np.array([0.2126, 0.7152, 0.0722])
out = np.empty_like(grid)
out[...] = OCEAN_RGB
out[lum < LAND_LUMA_CUT] = LAND_RGB
else:
out = grid
if BASEMAP_BLACK_POINT:
floor = np.percentile(out.reshape(-1, 3), BLACK_POINT_PCT, axis=0)
out = out - floor
texture = np.clip(out, 0, 1)
# Written aside and renamed, because np.save straight onto the final path
# leaves a truncated array at the exact name the next run trusts if the
# write is interrupted -- and the remedy is knowing to delete a file buried
# in the platform cache directory.
tmp = key.with_name(key.name + ".tmp")
with open(tmp, "wb") as fh:
np.save(fh, texture)
os.replace(tmp, key)
_BASEMAP[settings] = texture
return texture
def sample_texture(tex, lon_g, lat_g):
"""Bilinear sample of an equirectangular texture on a longitude/latitude grid.
Nearest-neighbour sampling reads one texel per mesh cell, so at any mesh coarser
than the texture most of the texture is simply skipped and coastlines break into
steps. Interpolating between the four surrounding texels uses the detail that is
there, which is what makes a finer mesh worth paying for.
Args:
tex: `(h, w, 3)` texture in 0..1, spanning -180..180 longitude and 90..-90 latitude.
lon_g: Longitudes to sample, in degrees.
lat_g: Latitudes to sample, in degrees, matching the shape of `lon_g`.
Returns:
An `(..., 3)` array of sampled colours in 0..1.
"""
h, w = tex.shape[:2]
fx = (lon_g + 180.0) / 360.0 * (w - 1)
fy = (90.0 - lat_g) / 180.0 * (h - 1)
x0 = np.clip(np.floor(fx).astype(int), 0, w - 1)
y0 = np.clip(np.floor(fy).astype(int), 0, h - 1)
x1, y1 = np.clip(x0 + 1, 0, w - 1), np.clip(y0 + 1, 0, h - 1)
tx, ty = (fx - x0)[..., None], (fy - y0)[..., None]
top = tex[y0, x0] * (1 - tx) + tex[y0, x1] * tx
bot = tex[y1, x0] * (1 - tx) + tex[y1, x1] * tx
return top * (1 - ty) + bot * ty
#: How finely the globe is meshed. This is the render's cost driver and it is
#: quadratic in the numbers below, so it dominates everything else: at 1920 px
#: wide a 1440x720 mesh works out at roughly one screen pixel per cell, which is
#: the point past which a finer mesh cannot be seen. It must be defined before
#: `textured_earth`, which takes it as a default argument.
MESH_LON, MESH_LAT = 1440, 720
def _earth_face(n_lon, n_lat):
"""The lon/lat mesh and its lit colours, which are the same in every frame.
Neither the mesh nor the sampled texture depends on `spin` -- rotating the
globe moves where those colours are drawn, not what they are. Sampling is a
million bilinear gathers, so doing it inside the draw loop repeated the
identical work on all 480 frames.
Args:
n_lon: Mesh columns, spanning -180..180 degrees.
n_lat: Mesh rows, spanning 90..-90 degrees.
Returns:
tuple: `(lon_grid, lat_grid, face_colours)`, the last already lit and
clipped to 0..1.
"""
key = (n_lon, n_lat, BASEMAP, BRIGHTNESS)
if key not in _EARTH_FACE:
lon_g, lat_g = np.meshgrid(
np.linspace(-180.0, 180.0, n_lon), np.linspace(90.0, -90.0, n_lat)
)
rgb = sample_texture(basemap_texture(), lon_g, lat_g)
# Evenly lit: no light direction, so no terminator and no night side.
# Every longitude reads the same, which keeps a satellite passing over
# the far limb as legible as one crossing the middle.
_EARTH_FACE[key] = (lon_g, lat_g, np.clip(rgb * BRIGHTNESS, 0, 1))
return _EARTH_FACE[key]
def textured_earth(ax, spin, n_lon=MESH_LON, n_lat=MESH_LAT):
"""Earth as an evenly lit, texture-mapped sphere -- real continents, tilted and spun.
The mesh is sampled at 1440 x 720 against a 2880 x 1440 basemap. A coarser mesh
was the source of the blockiness: the texture was always this detailed, the
sphere just was not sampling it.
Args:
ax: The 3D axis to draw onto.
spin: Earth's rotation angle for this frame, in radians.
n_lon: Mesh columns spanning -180..180 degrees.
n_lat: Mesh rows spanning 90..-90 degrees.
"""
lon_g, lat_g, face = _earth_face(n_lon, n_lat)
lo, la = np.deg2rad(lon_g) + spin, np.deg2rad(lat_g)
x, y, z = np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)
xr, zr = tilt_y(x, z)
ax.plot_surface(
xr,
y,
zr,
facecolors=face,
rstride=1,
cstride=1,
linewidth=0,
antialiased=False,
shade=False,
zorder=Z_GLOBE,
)
def render_master(fig, render, n_frames, stem, fps=24, bitrate=12000, dpi=160):
"""Draw every frame once, into the master mp4.
This is the only expensive step: `render` is called `n_frames` times and each call
draws a full frame. Every published size is then derived from the file this writes,
by `publish_clip`, which transcodes rather than redrawing.
Args:
fig: The figure `render` draws onto.
render: Callable taking a frame index and drawing that frame onto `fig`.
n_frames: How many frames to draw.
stem: Output path without a suffix; `.mp4` is written beside it.
fps: Frame rate of the master.
bitrate: Encoder bitrate in kbit/s. An mp4's size is bitrate times duration
regardless of resolution, so this has to rise with the frame size or the
detail the render spent hours on is simply compressed away.
dpi: Dots per inch, which together with the figure size sets the pixel
dimensions of the master.
Returns:
str: Path to the mp4 that was written.
"""
writer = manim.FFMpegWriter(fps=fps, bitrate=bitrate)
mp4 = f"{stem}.mp4"
with writer.saving(fig, mp4, dpi):
for i in range(n_frames):
render(i)
writer.grab_frame()
return mp4
def publish_clip(
master, web_width=1280, gif_width=880, gif_fps=12, webp_width=1080, webp_fps=15
):
"""Derive every published size from the rendered master.
Drawing the frames is the expensive part and it happens once, at the largest size
any target needs; every other file is a transcode of those same frames.
Three outputs, because the three places this clip appears do not accept the same
thing. The documentation site takes an mp4, which it can embed in a `<video>`.
GitHub strips video tags from README markdown, so the README needs an image --
and an animated webp beats a gif on every axis here, carrying truecolour at a
wider size for fewer bytes. The gif remains as the conservative fallback.
The gif settings are not defaults. A 256-colour palette measured 1.8 dB better
than 128 on this material, `stats_mode=full` suits a composition that is mostly
static, and dithering is off because error diffusion speckles a dark, smooth
frame -- measured better *and* smaller than the diffused version. 880 px matches
GitHub's README column, so nothing is scaled in either direction.
Args:
master: Path to the rendered master mp4.
web_width: Width of the mp4 the documentation site embeds.
gif_width: Width of the gif. Below about 700 px the satellite labels stop
being legible, so this is a floor rather than a preference.
gif_fps: Frame rate of the gif; halving it halves the file with no visible
judder, since the clip's duration is unchanged.
webp_width: Width of the animated webp used in the README.
webp_fps: Frame rate of the webp, which can afford more than the gif.
Returns:
tuple: The `(web_mp4, gif, webp)` paths that were written.
"""
master = Path(master)
exe = imageio_ffmpeg.get_ffmpeg_exe()
stem = str(master.with_suffix("")).rsplit("-1920", 1)[0]
web, gif, webp = f"{stem}-{web_width}.mp4", f"{stem}.gif", f"{stem}.webp"
# The palette is scratch, so it goes to the system temp directory rather
# than beside the published clips: the output directory is tracked, and a
# failed transcode used to strand a stray PNG in it that nothing ignores.
palette = os.path.join(tempfile.gettempdir(), f"{master.stem}_palette.png")
subprocess.run(
[
exe,
"-y",
"-v",
"error",
"-i",
str(master),
"-vf",
f"scale={web_width}:-2:flags=lanczos",
"-c:v",
"libx264",
"-b:v",
"6000k",
"-pix_fmt",
"yuv420p",
web,
],
check=True,
)
subprocess.run(
[
exe,
"-y",
"-v",
"error",
"-i",
str(master),
"-vf",
f"fps={webp_fps},scale={webp_width}:-1:flags=lanczos",
"-c:v",
"libwebp",
"-lossless",
"0",
"-q:v",
"72",
"-loop",
"0",
"-an",
webp,
],
check=True,
)
scale = f"fps={gif_fps},scale={gif_width}:-1:flags=lanczos"
try:
subprocess.run(
[
exe,
"-y",
"-v",
"error",
"-i",
str(master),
"-vf",
f"{scale},palettegen=max_colors=256:stats_mode=full",
palette,
],
check=True,
)
subprocess.run(
[
exe,
"-y",
"-v",
"error",
"-i",
str(master),
"-i",
palette,
"-lavfi",
f"{scale} [x]; [x][1:v] paletteuse=dither=none:diff_mode=rectangle",
gif,
],
check=True,
)
finally:
if os.path.exists(palette):
os.remove(palette)
return web, gif, webp
def _glyph_rect(x0, y0, x1, y1):
"""Vertices and codes for one closed rectangle of a compound marker path.
Args:
x0: Left edge, in marker units.
y0: Bottom edge, in marker units.
x1: Right edge, in marker units.
y1: Top edge, in marker units.
Returns:
tuple: The `(vertices, codes)` lists to concatenate into a `Path`.
"""
return (
[(x0, y0), (x1, y0), (x1, y1), (x0, y1), (x0, y0)],
[MPath.MOVETO, MPath.LINETO, MPath.LINETO, MPath.LINETO, MPath.CLOSEPOLY],
)
def _satellite_marker():
"""A satellite silhouette: two solar wings flanking a small bus on short booms.
Drawn as a marker rather than as 3D geometry on purpose. A modelled satellite would
be two or three pixels across in the wide framing and vanish; a marker is sized in
points, so it stays legible at every zoom.
Returns:
matplotlib.path.Path: The compound path, centred on the origin and
scaled to roughly one marker unit across.
"""
verts, codes = [], []
for box in (
(-0.17, -0.30, 0.17, 0.30), # bus
(-0.40, -0.035, -0.17, 0.035), # left boom
(0.17, -0.035, 0.40, 0.035), # right boom
(-1.00, -0.62, -0.40, 0.62), # left wing
(0.40, -0.62, 1.00, 0.62),
): # right wing
v, c = _glyph_rect(*box)
verts += v
codes += c
return MPath(verts, codes)
#: The glyph every satellite is drawn with.
SATELLITE = _satellite_marker()
The fleet¶
Each row is a real spacecraft earthlens can pull data from, with the backend that reaches it. Altitudes and inclinations are the published mission values; the sources are listed at the foot of this notebook.
kind drives how the orbit is drawn:
geo— geostationary, parked over a fixed longitude at 35,786 km.sso— sun-synchronous, near-polar, retrograde (inclination just past 90°).leo— everything else in low orbit, inclined but not sun-synchronous. ICESat-2 belongs here despite its near-polar 92° inclination: that value is chosen precisely so the orbit is not sun-synchronous.
In [ ]:
Copied!
#: One row per spacecraft earthlens can reach. `alt_km` and `incl_deg` are the
#: published mission values (sources at the foot of the notebook); `lon_deg` is
#: the parked longitude for geostationary satellites and `None` otherwise.
SATELLITES = [
# --- geostationary: parked over a fixed longitude at 35,786 km -----------
("GOES-19", "goes", "geo", 35786, 0.0, -75.2),
("GOES-18", "goes", "geo", 35786, 0.0, -137.0),
("GOES-16", "goes", "geo", 35786, 0.0, -104.7),
("Meteosat-12", "eumetsat", "geo", 35786, 0.0, 0.0),
("Meteosat-11", "eumetsat", "geo", 35786, 0.0, 9.5),
("Meteosat-9", "eumetsat", "geo", 35786, 0.0, 45.5),
("Himawari-9", "jaxa", "geo", 35786, 0.0, 140.7),
# --- sun-synchronous: near-polar, retrograde, fixed local crossing time --
("Sentinel-1A", "asf", "sso", 693, 98.18, None),
("Sentinel-1C", "asf", "sso", 693, 98.18, None),
("Sentinel-2A", "stac", "sso", 786, 98.62, None),
("Sentinel-2B", "stac", "sso", 786, 98.62, None),
("Sentinel-2C", "stac", "sso", 786, 98.62, None),
("Sentinel-3A", "eumetsat", "sso", 815, 98.60, None),
("Sentinel-3B", "eumetsat", "sso", 815, 98.60, None),
("Sentinel-5P", "eumetsat", "sso", 824, 98.74, None),
("Landsat 8", "stac", "sso", 705, 98.20, None),
("Landsat 9", "stac", "sso", 705, 98.20, None),
("Terra", "earthdata", "sso", 694, 98.20, None),
("Aqua", "earthdata", "sso", 705, 98.20, None),
("Suomi-NPP", "earthdata", "sso", 824, 98.74, None),
("NOAA-20", "firms", "sso", 824, 98.79, None),
("NOAA-21", "firms", "sso", 833, 98.80, None),
("Metop-B", "eumetsat", "sso", 817, 98.70, None),
("Metop-C", "eumetsat", "sso", 817, 98.70, None),
("GCOM-C", "jaxa", "sso", 798, 98.60, None),
("GCOM-W", "jaxa", "sso", 700, 98.20, None),
("ALOS-2", "asf", "sso", 628, 97.90, None),
("GOSAT-2", "jaxa", "sso", 613, 97.84, None),
("SMAP", "asf", "sso", 685, 98.00, None),
("PACE", "earthdata", "sso", 677, 98.00, None),
("OCO-2", "earthdata", "sso", 705, 98.20, None),
# ICESat-2 sits in this block because the rows are ordered by mission,
# not by kind, and the index sets each orbital plane. Its 92 deg
# inclination is deliberately NOT sun-synchronous: the orbit precesses
# so ATLAS samples the poles across every local time rather than one.
("ICESat-2", "earthdata", "leo", 496, 92.00, None),
# --- inclined low orbits: deliberately not sun-synchronous --------------
("GPM Core", "earthdata", "leo", 407, 65.00, None),
("Sentinel-6", "eumetsat", "leo", 1336, 66.00, None),
]
#: One colour per backend, so a viewer can read which provider reaches what.
BACKEND_COLOUR = {
"goes": "#ffb54d",
"eumetsat": "#7ec8ff",
"jaxa": "#ff8f8f",
"asf": "#b6f2b6",
"stac": "#d9b3ff",
"earthdata": "#ffe680",
"firms": "#ff9de0",
}
print(f"{len(SATELLITES)} satellites across {len(BACKEND_COLOUR)} backends")
for kind in ("geo", "sso", "leo"):
rows = [s for s in SATELLITES if s[2] == kind]
print(
f" {kind}: {len(rows):>2} ({', '.join(r[0] for r in rows[:4])}{', ...' if len(rows) > 4 else ''})"
)
#: One row per spacecraft earthlens can reach. `alt_km` and `incl_deg` are the
#: published mission values (sources at the foot of the notebook); `lon_deg` is
#: the parked longitude for geostationary satellites and `None` otherwise.
SATELLITES = [
# --- geostationary: parked over a fixed longitude at 35,786 km -----------
("GOES-19", "goes", "geo", 35786, 0.0, -75.2),
("GOES-18", "goes", "geo", 35786, 0.0, -137.0),
("GOES-16", "goes", "geo", 35786, 0.0, -104.7),
("Meteosat-12", "eumetsat", "geo", 35786, 0.0, 0.0),
("Meteosat-11", "eumetsat", "geo", 35786, 0.0, 9.5),
("Meteosat-9", "eumetsat", "geo", 35786, 0.0, 45.5),
("Himawari-9", "jaxa", "geo", 35786, 0.0, 140.7),
# --- sun-synchronous: near-polar, retrograde, fixed local crossing time --
("Sentinel-1A", "asf", "sso", 693, 98.18, None),
("Sentinel-1C", "asf", "sso", 693, 98.18, None),
("Sentinel-2A", "stac", "sso", 786, 98.62, None),
("Sentinel-2B", "stac", "sso", 786, 98.62, None),
("Sentinel-2C", "stac", "sso", 786, 98.62, None),
("Sentinel-3A", "eumetsat", "sso", 815, 98.60, None),
("Sentinel-3B", "eumetsat", "sso", 815, 98.60, None),
("Sentinel-5P", "eumetsat", "sso", 824, 98.74, None),
("Landsat 8", "stac", "sso", 705, 98.20, None),
("Landsat 9", "stac", "sso", 705, 98.20, None),
("Terra", "earthdata", "sso", 694, 98.20, None),
("Aqua", "earthdata", "sso", 705, 98.20, None),
("Suomi-NPP", "earthdata", "sso", 824, 98.74, None),
("NOAA-20", "firms", "sso", 824, 98.79, None),
("NOAA-21", "firms", "sso", 833, 98.80, None),
("Metop-B", "eumetsat", "sso", 817, 98.70, None),
("Metop-C", "eumetsat", "sso", 817, 98.70, None),
("GCOM-C", "jaxa", "sso", 798, 98.60, None),
("GCOM-W", "jaxa", "sso", 700, 98.20, None),
("ALOS-2", "asf", "sso", 628, 97.90, None),
("GOSAT-2", "jaxa", "sso", 613, 97.84, None),
("SMAP", "asf", "sso", 685, 98.00, None),
("PACE", "earthdata", "sso", 677, 98.00, None),
("OCO-2", "earthdata", "sso", 705, 98.20, None),
# ICESat-2 sits in this block because the rows are ordered by mission,
# not by kind, and the index sets each orbital plane. Its 92 deg
# inclination is deliberately NOT sun-synchronous: the orbit precesses
# so ATLAS samples the poles across every local time rather than one.
("ICESat-2", "earthdata", "leo", 496, 92.00, None),
# --- inclined low orbits: deliberately not sun-synchronous --------------
("GPM Core", "earthdata", "leo", 407, 65.00, None),
("Sentinel-6", "eumetsat", "leo", 1336, 66.00, None),
]
#: One colour per backend, so a viewer can read which provider reaches what.
BACKEND_COLOUR = {
"goes": "#ffb54d",
"eumetsat": "#7ec8ff",
"jaxa": "#ff8f8f",
"asf": "#b6f2b6",
"stac": "#d9b3ff",
"earthdata": "#ffe680",
"firms": "#ff9de0",
}
print(f"{len(SATELLITES)} satellites across {len(BACKEND_COLOUR)} backends")
for kind in ("geo", "sso", "leo"):
rows = [s for s in SATELLITES if s[2] == kind]
print(
f" {kind}: {len(rows):>2} ({', '.join(r[0] for r in rows[:4])}{', ...' if len(rows) > 4 else ''})"
)
From elements to positions¶
Two pieces of physics do all the work.
Kepler's third law gives the period from the orbital radius, so no period is invented:
$$T = 2\pi\sqrt{a^3/\mu}$$
with $\mu = 398{,}600.4418\ \mathrm{km^3/s^2}$. At geostationary radius this returns 23 h 56 min — one sidereal day — which is the check that the maths is right: the satellite's period matches Earth's rotation, so it hangs over one longitude.
The orbital plane is built from the inclination and the right ascension of the ascending node, then the satellite is placed on it by its phase angle. A geostationary satellite is the degenerate case: zero inclination, and the node fixed to the longitude it is parked over.
In [ ]:
Copied!
#: Radial compression. True scale puts the whole low-orbit fleet in a band 1%
#: of the frame wide while geostationary sits 6.6 radii out, so nothing reads.
#: Altitudes above the surface are raised to this power, which keeps the
#: ordering and the shape of the gap but pulls the two populations together.
RADIAL_GAMMA = 0.34
#: Where geostationary orbit is drawn, in Earth radii (its true value is 6.617).
GEO_DRAWN = 3.20
def draw_radius(true_r):
"""Compress a true orbital radius into the radius actually drawn.
Only the drawing is compressed: periods, inclinations and the relative
speeds all come from the true radius, so the physics on screen is unchanged.
Earth stays at 1.0 and the ordering of the orbits is preserved.
Args:
true_r: The real orbital radius, in Earth radii.
Returns:
float: The radius to draw at, in Earth radii.
Examples:
- Earth's surface is the fixed point of the mapping:
```python
>>> round(draw_radius(1.0), 3)
1.0
```
- Geostationary lands where GEO_DRAWN says:
```python
>>> round(draw_radius(6.617), 2)
3.2
```
"""
k = (GEO_DRAWN - 1.0) / (orbit_radius(35786) - 1.0) ** RADIAL_GAMMA
return 1.0 + k * max(true_r - 1.0, 0.0) ** RADIAL_GAMMA
def orbit_radius(alt_km):
"""Orbital radius in Earth radii, from an altitude above the surface.
Args:
alt_km: Altitude above mean sea level, in kilometres.
Returns:
float: The orbital radius expressed in Earth radii.
"""
return (R_EARTH_KM + alt_km) / R_EARTH_KM
def orbit_period_s(alt_km):
"""Orbital period from Kepler's third law.
Nothing here is tuned by eye: the relative speeds in the animation follow
from this, which is why a geostationary satellite ends up matching Earth's
rotation rather than being placed to look as though it does.
Args:
alt_km: Altitude above mean sea level, in kilometres.
Returns:
float: The period in seconds.
"""
a = R_EARTH_KM + alt_km
return 2.0 * np.pi * np.sqrt(a**3 / MU)
def plane_basis(incl_deg, raan_deg):
"""Two orthonormal vectors spanning an orbital plane, in tilted world axes.
The plane is built from its inclination and the right ascension of its
ascending node, then tilted with Earth so the poles line up with the drawn
globe.
Args:
incl_deg: Inclination to the equator, in degrees. Values above 90 are
retrograde, which is what makes an orbit sun-synchronous.
raan_deg: Right ascension of the ascending node, in degrees.
Returns:
tuple: `(u, v)` -- unit vectors; the satellite sits at
`r * (u cos(phase) + v sin(phase))`.
"""
i, o = np.deg2rad(incl_deg), np.deg2rad(raan_deg)
u = np.array([np.cos(o), np.sin(o), 0.0])
v = np.array(
[
-np.sin(o) * np.cos(i),
np.cos(o) * np.cos(i),
np.sin(i),
]
)
ux, uz = tilt_y(u[0], u[2])
vx, vz = tilt_y(v[0], v[2])
return np.array([ux, u[1], uz]), np.array([vx, v[1], vz])
def satellite_track(sat, index, n_planes):
"""Everything needed to place and draw one satellite over time.
Args:
sat: A row from `SATELLITES`.
index: The satellite's position in the list, used to spread the orbital
planes so the low-Earth shell reads as a shell rather than a line.
n_planes: How many satellites share the spreading.
Returns:
dict: Radius, period, plane basis and starting phase.
"""
name, backend, kind, alt_km, incl, lon = sat
r = orbit_radius(alt_km) # true, drives the period
r_draw = draw_radius(r) # compressed, drives the drawing
period = orbit_period_s(alt_km)
if kind == "geo":
# A parked longitude fixes the node; zero inclination makes the plane
# the equator itself, so the satellite hangs over one meridian.
raan, phase0 = lon, 0.0
else:
# Real inclination, presentational phasing: spread the nodes and the
# starting positions so the shell is legible.
raan = (index * 360.0 / n_planes) % 360.0
phase0 = (index * 2.399963) % (2 * np.pi) # golden angle, in radians
u, v = plane_basis(incl, raan)
return {
"name": name,
"backend": backend,
"kind": kind,
"r": r,
"r_draw": r_draw,
"period": period,
"u": u,
"v": v,
"phase0": phase0,
"lon": lon,
"colour": BACKEND_COLOUR.get(backend, "#cccccc"),
}
#: Plane index counted over the non-geostationary rows only. A geostationary
#: satellite takes its node from its parked longitude and throws the spread
#: value away, so counting all 34 rows spent the first seven slots on planes
#: nothing was drawn in -- leaving an 85-degree wedge of the shell empty, in a
#: spread whose whole purpose is to fill it evenly.
_PLANE_INDEX = {
sat[0]: j for j, sat in enumerate(s for s in SATELLITES if s[2] != "geo")
}
TRACKS = [
satellite_track(sat, _PLANE_INDEX.get(sat[0], 0), len(_PLANE_INDEX))
for sat in SATELLITES
]
#: The catalogue split by the one thing about a satellite that never changes.
#: Every frame draws these two populations separately -- beams, swaths and
#: near-side names belong to the low orbiters, the parked-longitude labels to the
#: geostationary ring -- so the split is taken once here rather than re-derived
#: inside the draw loop on each of the 480 frames.
GEO_TRACKS = [t for t in TRACKS if t["kind"] == "geo"]
LOW_ORBIT_TRACKS = [t for t in TRACKS if t["kind"] != "geo"]
geo = next(t for t in TRACKS if t["kind"] == "geo")
geo_check_h = 2.0 * np.pi * np.sqrt((R_EQUATORIAL_KM + 35786.0) ** 3 / MU) / 3600.0
print(f"geostationary period : {geo_check_h:.5f} h (one sidereal day = 23.93447 h)")
print(
f" {geo['period'] / 3600:.3f} h as drawn, on the mean radius"
)
print(f"geostationary radius : {geo['r']:.3f} true, drawn at {geo['r_draw']:.2f}")
low = min(TRACKS, key=lambda t: t["r"])
print(
f"lowest orbit : {low['name']} at {low['r']:.3f} true, "
f"drawn at {low['r_draw']:.2f}, {low['period'] / 60:.1f} min"
)
#: Radial compression. True scale puts the whole low-orbit fleet in a band 1%
#: of the frame wide while geostationary sits 6.6 radii out, so nothing reads.
#: Altitudes above the surface are raised to this power, which keeps the
#: ordering and the shape of the gap but pulls the two populations together.
RADIAL_GAMMA = 0.34
#: Where geostationary orbit is drawn, in Earth radii (its true value is 6.617).
GEO_DRAWN = 3.20
def draw_radius(true_r):
"""Compress a true orbital radius into the radius actually drawn.
Only the drawing is compressed: periods, inclinations and the relative
speeds all come from the true radius, so the physics on screen is unchanged.
Earth stays at 1.0 and the ordering of the orbits is preserved.
Args:
true_r: The real orbital radius, in Earth radii.
Returns:
float: The radius to draw at, in Earth radii.
Examples:
- Earth's surface is the fixed point of the mapping:
```python
>>> round(draw_radius(1.0), 3)
1.0
```
- Geostationary lands where GEO_DRAWN says:
```python
>>> round(draw_radius(6.617), 2)
3.2
```
"""
k = (GEO_DRAWN - 1.0) / (orbit_radius(35786) - 1.0) ** RADIAL_GAMMA
return 1.0 + k * max(true_r - 1.0, 0.0) ** RADIAL_GAMMA
def orbit_radius(alt_km):
"""Orbital radius in Earth radii, from an altitude above the surface.
Args:
alt_km: Altitude above mean sea level, in kilometres.
Returns:
float: The orbital radius expressed in Earth radii.
"""
return (R_EARTH_KM + alt_km) / R_EARTH_KM
def orbit_period_s(alt_km):
"""Orbital period from Kepler's third law.
Nothing here is tuned by eye: the relative speeds in the animation follow
from this, which is why a geostationary satellite ends up matching Earth's
rotation rather than being placed to look as though it does.
Args:
alt_km: Altitude above mean sea level, in kilometres.
Returns:
float: The period in seconds.
"""
a = R_EARTH_KM + alt_km
return 2.0 * np.pi * np.sqrt(a**3 / MU)
def plane_basis(incl_deg, raan_deg):
"""Two orthonormal vectors spanning an orbital plane, in tilted world axes.
The plane is built from its inclination and the right ascension of its
ascending node, then tilted with Earth so the poles line up with the drawn
globe.
Args:
incl_deg: Inclination to the equator, in degrees. Values above 90 are
retrograde, which is what makes an orbit sun-synchronous.
raan_deg: Right ascension of the ascending node, in degrees.
Returns:
tuple: `(u, v)` -- unit vectors; the satellite sits at
`r * (u cos(phase) + v sin(phase))`.
"""
i, o = np.deg2rad(incl_deg), np.deg2rad(raan_deg)
u = np.array([np.cos(o), np.sin(o), 0.0])
v = np.array(
[
-np.sin(o) * np.cos(i),
np.cos(o) * np.cos(i),
np.sin(i),
]
)
ux, uz = tilt_y(u[0], u[2])
vx, vz = tilt_y(v[0], v[2])
return np.array([ux, u[1], uz]), np.array([vx, v[1], vz])
def satellite_track(sat, index, n_planes):
"""Everything needed to place and draw one satellite over time.
Args:
sat: A row from `SATELLITES`.
index: The satellite's position in the list, used to spread the orbital
planes so the low-Earth shell reads as a shell rather than a line.
n_planes: How many satellites share the spreading.
Returns:
dict: Radius, period, plane basis and starting phase.
"""
name, backend, kind, alt_km, incl, lon = sat
r = orbit_radius(alt_km) # true, drives the period
r_draw = draw_radius(r) # compressed, drives the drawing
period = orbit_period_s(alt_km)
if kind == "geo":
# A parked longitude fixes the node; zero inclination makes the plane
# the equator itself, so the satellite hangs over one meridian.
raan, phase0 = lon, 0.0
else:
# Real inclination, presentational phasing: spread the nodes and the
# starting positions so the shell is legible.
raan = (index * 360.0 / n_planes) % 360.0
phase0 = (index * 2.399963) % (2 * np.pi) # golden angle, in radians
u, v = plane_basis(incl, raan)
return {
"name": name,
"backend": backend,
"kind": kind,
"r": r,
"r_draw": r_draw,
"period": period,
"u": u,
"v": v,
"phase0": phase0,
"lon": lon,
"colour": BACKEND_COLOUR.get(backend, "#cccccc"),
}
#: Plane index counted over the non-geostationary rows only. A geostationary
#: satellite takes its node from its parked longitude and throws the spread
#: value away, so counting all 34 rows spent the first seven slots on planes
#: nothing was drawn in -- leaving an 85-degree wedge of the shell empty, in a
#: spread whose whole purpose is to fill it evenly.
_PLANE_INDEX = {
sat[0]: j for j, sat in enumerate(s for s in SATELLITES if s[2] != "geo")
}
TRACKS = [
satellite_track(sat, _PLANE_INDEX.get(sat[0], 0), len(_PLANE_INDEX))
for sat in SATELLITES
]
#: The catalogue split by the one thing about a satellite that never changes.
#: Every frame draws these two populations separately -- beams, swaths and
#: near-side names belong to the low orbiters, the parked-longitude labels to the
#: geostationary ring -- so the split is taken once here rather than re-derived
#: inside the draw loop on each of the 480 frames.
GEO_TRACKS = [t for t in TRACKS if t["kind"] == "geo"]
LOW_ORBIT_TRACKS = [t for t in TRACKS if t["kind"] != "geo"]
geo = next(t for t in TRACKS if t["kind"] == "geo")
geo_check_h = 2.0 * np.pi * np.sqrt((R_EQUATORIAL_KM + 35786.0) ** 3 / MU) / 3600.0
print(f"geostationary period : {geo_check_h:.5f} h (one sidereal day = 23.93447 h)")
print(
f" {geo['period'] / 3600:.3f} h as drawn, on the mean radius"
)
print(f"geostationary radius : {geo['r']:.3f} true, drawn at {geo['r_draw']:.2f}")
low = min(TRACKS, key=lambda t: t["r"])
print(
f"lowest orbit : {low['name']} at {low['r']:.3f} true, "
f"drawn at {low['r_draw']:.2f}, {low['period'] / 60:.1f} min"
)
In [ ]:
Copied!
def position(track, t_hours, spin):
"""Where one satellite is at time `t_hours`.
Args:
track: A row from `TRACKS`.
t_hours: Hours since the start of the clip.
spin: Earth's rotation angle at that moment, in radians. Geostationary
satellites are carried by it; everything else ignores it.
Returns:
numpy.ndarray: The `(x, y, z)` position in Earth radii.
"""
if track["kind"] == "geo":
# Turning with Earth is the whole point, so drive the angle from the
# surface rotation rather than integrating the orbit separately. The
# equatorial plane plus the parked longitude is all it takes.
u, v = plane_basis(0.0, 0.0)
ang = spin + np.deg2rad(track["lon"])
return track["r_draw"] * (u * np.cos(ang) + v * np.sin(ang))
ang = track["phase0"] + 2 * np.pi * (t_hours * 3600.0) / track["period"]
return track["r_draw"] * (track["u"] * np.cos(ang) + track["v"] * np.sin(ang))
def orbit_ring(track, n=240):
"""The full closed orbit as a polyline, for drawing the path itself.
Args:
track: A row from `TRACKS`.
n: How many points to sample around the circle.
Returns:
tuple: `(x, y, z)` arrays tracing the orbit once.
"""
ang = np.linspace(0.0, 2 * np.pi, n)
pts = track["r_draw"] * (
np.outer(np.cos(ang), track["u"]) + np.outer(np.sin(ang), track["v"])
)
return pts[:, 0], pts[:, 1], pts[:, 2]
def position(track, t_hours, spin):
"""Where one satellite is at time `t_hours`.
Args:
track: A row from `TRACKS`.
t_hours: Hours since the start of the clip.
spin: Earth's rotation angle at that moment, in radians. Geostationary
satellites are carried by it; everything else ignores it.
Returns:
numpy.ndarray: The `(x, y, z)` position in Earth radii.
"""
if track["kind"] == "geo":
# Turning with Earth is the whole point, so drive the angle from the
# surface rotation rather than integrating the orbit separately. The
# equatorial plane plus the parked longitude is all it takes.
u, v = plane_basis(0.0, 0.0)
ang = spin + np.deg2rad(track["lon"])
return track["r_draw"] * (u * np.cos(ang) + v * np.sin(ang))
ang = track["phase0"] + 2 * np.pi * (t_hours * 3600.0) / track["period"]
return track["r_draw"] * (track["u"] * np.cos(ang) + track["v"] * np.sin(ang))
def orbit_ring(track, n=240):
"""The full closed orbit as a polyline, for drawing the path itself.
Args:
track: A row from `TRACKS`.
n: How many points to sample around the circle.
Returns:
tuple: `(x, y, z)` arrays tracing the orbit once.
"""
ang = np.linspace(0.0, 2 * np.pi, n)
pts = track["r_draw"] * (
np.outer(np.cos(ang), track["u"]) + np.outer(np.sin(ang), track["v"])
)
return pts[:, 0], pts[:, 1], pts[:, 2]
In [ ]:
Copied!
BRAND_DIR = Path("../../_images/branding/earthlens-brand-kit/logo")
BRAND_MARK = BRAND_DIR / "earthlens-lockup-stacked-overlay-full-transparent.png"
# the transparent variant: this clip is black under the mark (mean luminance 0.8), so the
# shadow alone separates it and no plate is needed -- a plate would read as a box on black
LOGO_FRAC = 0.11 # of figure width. The geostationary ring sweeps into both
# bottom corners as it turns and grazes the mark at this size -- about 2% of the
# mark's box -- but what crosses is one thin orbit line drawn at alpha 0.05-0.25,
# passing behind an opaque logo. Sizing down to clear it entirely put the
# wordmark at roughly 50 px in the 720-wide gif, which is the copy most likely to
# be reposted without any surrounding attribution, and too small to read there.
LOGO_MARGIN_Y = 0.0 # of figure height -- tucked low, clear of the orbit rings
LOGO_MARGIN = 0.025
def _mark_with_shadow(path, blur=0.065, strength=0.9):
"""The mark pre-composited over a blurred dark copy of itself.
A shadow rather than a panel: the mark keeps its transparent background, but the dark halo
gives the near-white wordmark something to sit against. Over sunlit cloud the bare mark
measures ~1 contrast, which is invisible.
Args:
path: Path to the RGBA lockup.
blur: Shadow blur radius, as a fraction of mark width.
strength: Shadow opacity, 0..1.
Returns:
tuple: `(rgba, grow)` -- the composited image and how much wider it is than the mark,
so the caller can size the axes by the mark rather than by the padded canvas.
"""
mark = plt.imread(path)
if mark.shape[2] == 3: # no alpha: treat black as transparent
mark = np.dstack([mark, mark.max(axis=2)])
# Both the padding and the blur are fractions of the *mark*, measured before
# the padding is added. Taking sigma from mark.shape[1] afterwards made the
# effective blur 1.3x the number passed in, and left the pad at 2.3 sigma,
# clipping the shadow's tail. The default is raised to match the sigma the
# mark was tuned against, so it looks the same and now means what it says.
width = mark.shape[1]
pad = int(round(width * blur * 3))
mark = np.pad(mark, ((pad, pad), (pad, pad), (0, 0)))
shade = gaussian_filter(mark[..., 3], sigma=width * blur) * strength
out = np.zeros_like(mark)
out[..., 3] = np.clip(mark[..., 3] + shade * (1.0 - mark[..., 3]), 0.0, 1.0)
share = np.divide(mark[..., 3], np.maximum(out[..., 3], 1e-6))
out[..., :3] = mark[..., :3] * share[..., None] # shadow contributes black
return out, mark.shape[1] / width
def stamp_logo(fig, frac=LOGO_FRAC, margin=LOGO_MARGIN):
"""Draw the earthlens mark in the figure's bottom-right corner.
Sized as a fraction of the figure, not in pixels, because the MP4 and the GIF are saved
from the same figure at different dpi -- a pixel size would come out inconsistent.
Args:
fig: The figure to stamp. The mark is added as its own axes, so it sits
above everything already drawn.
frac: Mark width as a fraction of figure width.
margin: Gap from the figure edge, as a fraction of figure width.
Returns:
matplotlib.axes.Axes: The inset axes holding the mark.
"""
rgba, grow = _mark_with_shadow(BRAND_MARK)
fig_w, fig_h = fig.get_size_inches()
width = frac * grow
height = width * (rgba.shape[0] / rgba.shape[1]) * (fig_w / fig_h)
inset = fig.add_axes([1.0 - margin - width, LOGO_MARGIN_Y, width, height], zorder=6)
inset.set_axis_off()
inset.patch.set_alpha(0.0)
inset.imshow(rgba, interpolation="antialiased")
return inset
BRAND_DIR = Path("../../_images/branding/earthlens-brand-kit/logo")
BRAND_MARK = BRAND_DIR / "earthlens-lockup-stacked-overlay-full-transparent.png"
# the transparent variant: this clip is black under the mark (mean luminance 0.8), so the
# shadow alone separates it and no plate is needed -- a plate would read as a box on black
LOGO_FRAC = 0.11 # of figure width. The geostationary ring sweeps into both
# bottom corners as it turns and grazes the mark at this size -- about 2% of the
# mark's box -- but what crosses is one thin orbit line drawn at alpha 0.05-0.25,
# passing behind an opaque logo. Sizing down to clear it entirely put the
# wordmark at roughly 50 px in the 720-wide gif, which is the copy most likely to
# be reposted without any surrounding attribution, and too small to read there.
LOGO_MARGIN_Y = 0.0 # of figure height -- tucked low, clear of the orbit rings
LOGO_MARGIN = 0.025
def _mark_with_shadow(path, blur=0.065, strength=0.9):
"""The mark pre-composited over a blurred dark copy of itself.
A shadow rather than a panel: the mark keeps its transparent background, but the dark halo
gives the near-white wordmark something to sit against. Over sunlit cloud the bare mark
measures ~1 contrast, which is invisible.
Args:
path: Path to the RGBA lockup.
blur: Shadow blur radius, as a fraction of mark width.
strength: Shadow opacity, 0..1.
Returns:
tuple: `(rgba, grow)` -- the composited image and how much wider it is than the mark,
so the caller can size the axes by the mark rather than by the padded canvas.
"""
mark = plt.imread(path)
if mark.shape[2] == 3: # no alpha: treat black as transparent
mark = np.dstack([mark, mark.max(axis=2)])
# Both the padding and the blur are fractions of the *mark*, measured before
# the padding is added. Taking sigma from mark.shape[1] afterwards made the
# effective blur 1.3x the number passed in, and left the pad at 2.3 sigma,
# clipping the shadow's tail. The default is raised to match the sigma the
# mark was tuned against, so it looks the same and now means what it says.
width = mark.shape[1]
pad = int(round(width * blur * 3))
mark = np.pad(mark, ((pad, pad), (pad, pad), (0, 0)))
shade = gaussian_filter(mark[..., 3], sigma=width * blur) * strength
out = np.zeros_like(mark)
out[..., 3] = np.clip(mark[..., 3] + shade * (1.0 - mark[..., 3]), 0.0, 1.0)
share = np.divide(mark[..., 3], np.maximum(out[..., 3], 1e-6))
out[..., :3] = mark[..., :3] * share[..., None] # shadow contributes black
return out, mark.shape[1] / width
def stamp_logo(fig, frac=LOGO_FRAC, margin=LOGO_MARGIN):
"""Draw the earthlens mark in the figure's bottom-right corner.
Sized as a fraction of the figure, not in pixels, because the MP4 and the GIF are saved
from the same figure at different dpi -- a pixel size would come out inconsistent.
Args:
fig: The figure to stamp. The mark is added as its own axes, so it sits
above everything already drawn.
frac: Mark width as a fraction of figure width.
margin: Gap from the figure edge, as a fraction of figure width.
Returns:
matplotlib.axes.Axes: The inset axes holding the mark.
"""
rgba, grow = _mark_with_shadow(BRAND_MARK)
fig_w, fig_h = fig.get_size_inches()
width = frac * grow
height = width * (rgba.shape[0] / rgba.shape[1]) * (fig_w / fig_h)
inset = fig.add_axes([1.0 - margin - width, LOGO_MARGIN_Y, width, height], zorder=6)
inset.set_axis_off()
inset.patch.set_alpha(0.0)
inset.imshow(rgba, interpolation="antialiased")
return inset
In [ ]:
Copied!
#: How far out each geostationary label sits, as a multiple of the orbit
#: radius. Every label sits at the same offset so none looks detached from
#: its satellite. Meteosat-12 (0 deg) and Meteosat-11 (9.5E) are close
#: enough on the ring that their labels can touch when both face the
#: camera; the map is kept as the hook for staggering a pair if that ever
#: needs solving again.
LABEL_OFFSET = {
"GOES-19": 1.10,
"GOES-18": 1.10,
"GOES-16": 1.10,
"Himawari-9": 1.10,
"Meteosat-12": 1.10,
"Meteosat-9": 1.10,
"Meteosat-11": 1.10,
}
#: Earth's rotation axis in scene coordinates -- the globe is drawn spun about
#: z and then tilted, so the pole ends up here. A patch of ground painted at one
#: moment has to be carried around this axis to stay stuck to the surface.
SPIN_AXIS = np.array([np.sin(TILT), 0.0, np.cos(TILT)])
#: Across-track swath of the instrument each mission is flown for, in km. These
#: differ by more than an order of magnitude, which is the point: VIIRS sweeps a
#: sixth of the planet's circumference in one pass while OCO-2 and ICESat-2 draw
#: a line thinner than the coastlines. Missions absent from this map get
#: `DEFAULT_SWATH_KM`.
SWATH_KM = {
"Sentinel-1A": 250,
"Sentinel-1C": 250,
"Sentinel-2A": 290,
"Sentinel-2B": 290,
"Sentinel-2C": 290,
"Sentinel-3A": 1270,
"Sentinel-3B": 1270,
"Sentinel-5P": 2600,
"Landsat 8": 185,
"Landsat 9": 185,
"Terra": 2330,
"Aqua": 2330,
"Suomi-NPP": 3060,
"NOAA-20": 3060,
"NOAA-21": 3060,
"Metop-B": 2900,
"Metop-C": 2900,
"GCOM-C": 1150,
"GCOM-W": 1450,
"ALOS-2": 350,
"GOSAT-2": 920,
"SMAP": 1000,
"PACE": 2663,
"OCO-2": 10,
"ICESat-2": 6,
"GPM Core": 885,
"Sentinel-6": 30,
}
DEFAULT_SWATH_KM = 300
#: How much of the recent past each trail keeps, in minutes, and how finely it
#: is sampled. A whole orbit is about 100 minutes, so this is roughly a third of
#: a circuit -- enough to read the track's curve without wrapping the globe.
TRAIL_MINUTES = 14.0
TRAIL_SAMPLES = 72
TRAIL_PEAK_ALPHA = 0.52
def rotate_about(vecs, axis, ang):
"""Rotate row vectors about `axis` by `ang` radians (Rodrigues' formula).
Args:
vecs: An `(n, 3)` array of vectors.
axis: A unit vector to rotate around.
ang: An `(n,)` array of angles, one per row of `vecs`.
Returns:
numpy.ndarray: The rotated `(n, 3)` array.
"""
c, s = np.cos(ang)[:, None], np.sin(ang)[:, None]
dot = (vecs @ axis)[:, None]
return vecs * c + np.cross(axis, vecs) * s + axis * dot * (1.0 - c)
#: Largest across-track step drawn as a straight chord, in radians. A chord
#: across an angle `a` sags below the surface by 1 - cos(a/2), so at 0.06 rad
#: that is 25 km on a 6371 km sphere -- under a screen pixel here. A 3060 km
#: swath spans 0.48 rad in one piece and sags 170 km, which draws the band as a
#: plate cutting into the planet rather than a strip lying on it.
MAX_CHORD_RAD = 0.06
def _across_steps(half):
"""How many spans an across-track half-angle should be cut into.
Args:
half: Half the swath's angular width, in radians.
Returns:
int: Span count, at least 1. Narrow swaths stay a single quad, so the
subdivision costs nothing where it buys nothing.
"""
return max(1, int(np.ceil(2.0 * half / MAX_CHORD_RAD)))
def ground_swath(track, t_now):
"""The strip of ground `track` has swept over the last `TRAIL_MINUTES`.
The sub-satellite point is the satellite's direction from Earth's centre. Each
sample is then carried around `SPIN_AXIS` by however far Earth has turned since
it was laid down, so the strip stays fixed to the surface rather than sliding
with the satellite. Samples on the far side of the globe are dropped, which is
what keeps the strip from painting through Earth.
Args:
track: A row from `TRACKS`.
t_now: Hours since the start of the clip.
Returns:
tuple: `(quads, alphas)` for a `Poly3DCollection`, or `None` before the
trail has any length.
"""
t0 = t_now - TRAIL_MINUTES / 60.0
if t_now - max(t0, 0.0) < 1e-4:
return None
ts = np.linspace(max(t0, 0.0), t_now, TRAIL_SAMPLES)
spins = 2 * np.pi * ts / 23.9344696
pts = np.array([position(track, t, s) for t, s in zip(ts, spins)])
d = pts / np.linalg.norm(pts, axis=1)[:, None]
d = rotate_about(d, SPIN_AXIS, spins[-1] - spins)
# Cross-track direction: perpendicular to both the ground point and the
# direction of travel along the track.
step = np.gradient(d, axis=0)
side = np.cross(d, step)
norm = np.linalg.norm(side, axis=1)[:, None]
side = np.divide(side, norm, out=np.zeros_like(side), where=norm > 1e-12)
half = SWATH_KM.get(track["name"], DEFAULT_SWATH_KM) / 2.0 / R_EARTH_KM
n_across = _across_steps(half)
offsets = np.linspace(-half, half, n_across + 1)
edges = [
(d * np.cos(a) + side * np.sin(a)) * 1.002 for a in offsets
] # each follows the sphere
front = d @ EYE > 0.02 # drop the far side, and the limb where it degenerates
age = np.linspace(0.0, 1.0, TRAIL_SAMPLES) ** 1.1 * TRAIL_PEAK_ALPHA
quads, alphas = [], []
for k in range(TRAIL_SAMPLES - 1):
if not (front[k] and front[k + 1]):
continue
for m in range(n_across):
quads.append(
[edges[m][k], edges[m + 1][k], edges[m + 1][k + 1], edges[m][k + 1]]
)
alphas.append(age[k])
if not quads:
return None
return np.array(quads), np.array(alphas)
#: How wide the top of a beam is, as a fraction of its footprint width, and how
#: opaque the sheet is. A push-broom instrument images a line across its track, so
#: its field of view is a flat fan rather than a cone: nearly a line at the
#: spacecraft, opening out to the full swath at the ground. Giving the top a small
#: but non-zero width makes it read as a trapezoid instead of a spike.
BEAM_TOP_FRAC = 0.07
BEAM_ALPHA = 0.15
def sensor_beam(track, t_hours, spin):
"""The trapezoidal field of view between a satellite and the ground it is imaging.
Args:
track: A row from `TRACKS`.
t_hours: Hours since the start of the clip.
spin: Earth's rotation angle at that moment, in radians.
Returns:
numpy.ndarray: The trapezoid's four corners, or `None` when the footprint
is on the far side of the globe and the sheet would draw through Earth.
"""
dt = 1.0 / 3600.0 # one second on
spin_ahead = spin + 2 * np.pi * dt / 23.9344696
here = position(track, t_hours, spin)
ahead = position(track, t_hours + dt, spin_ahead)
down = here / np.linalg.norm(here)
if float(down @ EYE) <= 0.02:
return None
# Both samples carried into the same Earth-fixed frame before differencing.
# Taking `ahead` at the current spin measured the inertial velocity, while
# ground_swath measures the ground track -- so the beam's edge and the
# leading edge of the trail it paints were skewed by a few degrees against
# each other, which is exactly the angle a beam should not be wrong by.
d_ahead = ahead / np.linalg.norm(ahead)
d_ahead = rotate_about(d_ahead[None, :], SPIN_AXIS, np.array([spin - spin_ahead]))[
0
]
across = np.cross(down, d_ahead - down)
norm = np.linalg.norm(across)
if norm < 1e-12:
return None
across /= norm
half = SWATH_KM.get(track["name"], DEFAULT_SWATH_KM) / 2.0 / R_EARTH_KM
# The ground edge is an arc, cut finely enough that it lies on the surface
# instead of chording under it.
n_across = _across_steps(half)
arc = np.array(
[
(down * np.cos(a) + across * np.sin(a)) * 1.004
for a in np.linspace(half, -half, n_across + 1)
]
)
top = across * (BEAM_TOP_FRAC * half)
return np.vstack([[here - top], [here + top], arc])
#: Depth band over which a low-orbit label fades in, as a fraction of the orbit
#: radius toward the camera. Below the first value the satellite is on the far
#: side and unlabelled; above the second the label is at full strength.
NEAR_LABEL_FADE = (0.10, 0.45)
HOURS = 1.5 # Everything runs at its true rate; the window is short instead.
# A low orbiter covers a little under one circuit here, which is slow enough to
# follow, and the ratio between orbit rate and Earth's spin stays correct -- so
# the ground tracks are real. Slowing the satellites instead would have shown
# the same on-screen speed (0.25 x 6.0 h = 1.0 x 1.5 h) but drawn swaths that
# drift west four times too fast.
N_FRAMES = 480
CLIP_FPS = 24
#: How much faster than real time the clip runs -- the simulated window divided by
#: the screen time. It is stated on the clock because a bare elapsed time does not
#: tell a viewer what to do with it: the point is that a low orbiter completing one
#: circuit as the clock reaches 1.5 h is showing them its orbital period.
SPEEDUP = round(HOURS * 3600 / (N_FRAMES / CLIP_FPS))
#: RENDER_DPI 160 on the 12x8 inch figure gives a 1920x1280 master. The bitrate
#: has to rise with the resolution: an mp4's size is bitrate x duration no matter
#: how many pixels are in it, so at the 3600 default a 1920-wide clip simply
#: compresses away the mesh detail the render spent hours producing.
RENDER_DPI = 160
RENDER_BITRATE = 12000
ELEV, AZIM = 22.0, -58.0
#: Unit vector from the origin toward the camera. Projecting a point onto it
#: gives how far that point lies toward the viewer: +1 straight in front of
#: Earth, 0 on the limb, -1 straight behind. Both the ring fade and the
#: near-side label test are this projection, so they agree with what is
#: actually in front of what.
EYE = np.array(
[
np.cos(np.deg2rad(ELEV)) * np.cos(np.deg2rad(AZIM)),
np.cos(np.deg2rad(ELEV)) * np.sin(np.deg2rad(AZIM)),
np.sin(np.deg2rad(ELEV)),
]
)
#: The painting order. Matplotlib's own 3D depth sorting is off
#: (computed_zorder = False), so every layer's place is stated here rather than
#: left to the order the calls happen to be made in. The globe sits in the
#: middle: anything on the far side goes below it, anything nearer goes above.
Z_BEHIND_GLOBE = 2
Z_GLOBE = 4
Z_NEAR_RING = 4.5
Z_SURFACE = 5 # swaths and beams, lying on the near face
Z_NEAR_SAT = 6
Z_LABEL = 7
#: The geostationary ring is drawn once, in a tone of its own: it is shared by
#: seven satellites, so no single provider's colour is the honest one for it.
GEO_RING_COLOUR = "#93a4b8"
LIM = 3.35 # tight around the drawn geostationary ring
fig = plt.figure(figsize=(12, 8), facecolor="black")
# Fitted so the scene fills 93% x 83% of the frame, centred, with the
# geostationary labels just inside the edges rather than clipped. The
# projected scene is a wide, shallow ellipse, so the axes overflows the
# figure on every side to leave no dead margin.
ax = fig.add_axes([-0.345, -0.437, 1.665, 1.873], projection="3d")
# mplot3d sorts artists by mean depth, which lets the globe paint over orbits
# in front of it; honour the zorder the drawing assigns instead.
ax.computed_zorder = False
def draw_frame(i):
"""Draw one frame: Earth, every orbit, and every satellite on it.
Clears and redraws the shared axis in place rather than building a new one,
so the figure the writer is grabbing from stays the same object throughout
the render.
Args:
i: Frame index, from 0 to `N_FRAMES - 1`. Time and Earth's rotation are
both derived from it, so the frame is fully determined by this.
Returns:
tuple: Empty. Matplotlib's animation callbacks return the artists to
re-blit; this redraws the whole axis instead, so there are none.
"""
ax.clear()
ax.set_axis_off()
ax.set_facecolor("black")
ax.set_xlim(-LIM, LIM)
ax.set_ylim(-LIM, LIM)
ax.set_zlim(-LIM * 0.62, LIM * 0.62)
ax.set_box_aspect((1, 1, 0.62))
ax.view_init(elev=ELEV, azim=AZIM)
t_hours = HOURS * i / (N_FRAMES - 1)
spin = 2 * np.pi * t_hours / 23.9344696 # one sidereal day
textured_earth(ax, spin)
# The orbit paths, faded continuously by depth so the far side of each ring
# recedes instead of competing with the near side. 34 rings at similar radii
# otherwise read as one flat tangle. A per-segment alpha ramp is used rather
# than a near/far split, which leaves a hard seam where the two halves meet.
#
# Each ring is split by depth against the globe's own zorder, though: a ring
# is a closed loop with half of it behind Earth, so drawing the whole thing
# at one zorder is wrong either way round. Held entirely behind, the near
# half vanished into the planet it should pass in front of.
geo_ring_drawn = False
for track in TRACKS:
# plane_basis(0, lon) rotates one equatorial circle, so every
# geostationary row traces identical geometry. Drawing all seven stacked
# them into a line about seven times its own alpha -- the brightest thing
# in the frame, and colour-neutral from the overlap rather than any
# provider's colour. One ring, drawn once, in a neutral tone.
if track["kind"] == "geo":
if geo_ring_drawn:
continue
geo_ring_drawn = True
colour = GEO_RING_COLOUR if track["kind"] == "geo" else track["colour"]
ox, oy, oz = orbit_ring(track)
pts = np.stack([ox, oy, oz], axis=1)
segs = np.stack([pts[:-1], pts[1:]], axis=1)
depth = ox * EYE[0] + oy * EYE[1] + oz * EYE[2]
depth = (depth[:-1] + depth[1:]) / 2 / track["r_draw"] # -1 far .. +1 near
rgba = np.tile(to_rgba(colour), (len(segs), 1))
rgba[:, 3] = 0.045 + 0.20 * (depth + 1) / 2
near = depth > 0
for mask, z in ((~near, Z_BEHIND_GLOBE), (near, Z_NEAR_RING)):
if mask.any():
ax.add_collection3d(
Line3DCollection(
segs[mask], colors=rgba[mask], linewidths=0.20, zorder=z
)
)
# The instrument's field of view: a flat sheet from the spacecraft down to the
# strip it is imaging right now. Drawn before the trail so the trail's older,
# fainter ground reads as behind it.
for track in LOW_ORBIT_TRACKS:
quad = sensor_beam(track, t_hours, spin)
if quad is None:
continue
ax.add_collection3d(
Poly3DCollection(
[quad],
facecolors=[to_rgba(track["colour"], BEAM_ALPHA)],
edgecolors="none",
zorder=Z_SURFACE,
)
)
# The ground each low orbiter has swept recently, fading out behind it. This
# is the point of a sun-synchronous orbit: a narrow strip at a time, a new one
# every pass. Geostationary satellites are skipped -- they stare at one face.
for track in LOW_ORBIT_TRACKS:
swept = ground_swath(track, t_hours)
if swept is None:
continue
quads, alphas = swept
rgba = np.tile(to_rgba(track["colour"]), (len(quads), 1))
rgba[:, 3] = alphas
ax.add_collection3d(
Poly3DCollection(
quads, facecolors=rgba, edgecolors="none", zorder=Z_SURFACE
)
)
# A satellite on the far side belongs behind the globe. Painted above it,
# roughly a third of the fleet floated over the disc at any moment, which
# reads as a flat sticker rather than an orbit. Held below the globe's
# zorder it is hidden only where Earth actually covers it -- one that is
# behind but outside the limb still shows, which is the truth of it.
for track in TRACKS:
p = position(track, t_hours, spin)
ax.scatter(
[p[0]],
[p[1]],
[p[2]],
marker=SATELLITE,
s=70,
c=track["colour"],
edgecolors="none",
depthshade=False,
zorder=Z_NEAR_SAT if float(p @ EYE) > 0 else Z_BEHIND_GLOBE,
)
# Name the low-orbit satellites on the near side only. Labelling all 27
# doubles the text on screen and the far half is already dimmed, so the
# names would sit on satellites the eye is not following. A satellite fades
# its label in as it comes round, rather than popping, so the set of names
# changes smoothly instead of flickering.
for track in LOW_ORBIT_TRACKS:
p = position(track, t_hours, spin)
depth = (p[0] * EYE[0] + p[1] * EYE[1] + p[2] * EYE[2]) / track["r_draw"]
if depth <= NEAR_LABEL_FADE[0]:
continue
a = np.clip(
(depth - NEAR_LABEL_FADE[0]) / (NEAR_LABEL_FADE[1] - NEAR_LABEL_FADE[0]),
0.0,
1.0,
)
ax.text(
p[0],
p[1],
p[2] + 0.085,
track["name"],
color=track["colour"],
fontsize=7.0,
ha="center",
va="bottom",
alpha=float(a),
zorder=Z_LABEL,
path_effects=[
pe.withStroke(linewidth=2.0, foreground="black", alpha=float(a))
],
)
# Name the geostationary satellites. Meteosat-12/-9/-11 sit within 9.5 deg
# of each other; LABEL_OFFSET is the hook for pushing their labels to
# of overlapping on the ring.
for track in GEO_TRACKS:
p = position(track, t_hours, spin)
out = LABEL_OFFSET.get(track["name"], 1.10)
ax.text(
p[0] * out,
p[1] * out,
p[2] * out + 0.06,
track["name"],
color=track["colour"],
fontsize=8.5,
ha="center",
va="center",
zorder=Z_LABEL,
path_effects=[pe.withStroke(linewidth=2.2, foreground="black")],
)
ax.text2D(
0.022,
0.035,
f"{t_hours:4.1f} h simulated · {SPEEDUP}× real time",
transform=fig.transFigure,
color="#9fb4c7",
fontsize=13,
ha="left",
family="monospace",
)
return ()
stamp_logo(fig)
draw_frame(0)
fig
#: How far out each geostationary label sits, as a multiple of the orbit
#: radius. Every label sits at the same offset so none looks detached from
#: its satellite. Meteosat-12 (0 deg) and Meteosat-11 (9.5E) are close
#: enough on the ring that their labels can touch when both face the
#: camera; the map is kept as the hook for staggering a pair if that ever
#: needs solving again.
LABEL_OFFSET = {
"GOES-19": 1.10,
"GOES-18": 1.10,
"GOES-16": 1.10,
"Himawari-9": 1.10,
"Meteosat-12": 1.10,
"Meteosat-9": 1.10,
"Meteosat-11": 1.10,
}
#: Earth's rotation axis in scene coordinates -- the globe is drawn spun about
#: z and then tilted, so the pole ends up here. A patch of ground painted at one
#: moment has to be carried around this axis to stay stuck to the surface.
SPIN_AXIS = np.array([np.sin(TILT), 0.0, np.cos(TILT)])
#: Across-track swath of the instrument each mission is flown for, in km. These
#: differ by more than an order of magnitude, which is the point: VIIRS sweeps a
#: sixth of the planet's circumference in one pass while OCO-2 and ICESat-2 draw
#: a line thinner than the coastlines. Missions absent from this map get
#: `DEFAULT_SWATH_KM`.
SWATH_KM = {
"Sentinel-1A": 250,
"Sentinel-1C": 250,
"Sentinel-2A": 290,
"Sentinel-2B": 290,
"Sentinel-2C": 290,
"Sentinel-3A": 1270,
"Sentinel-3B": 1270,
"Sentinel-5P": 2600,
"Landsat 8": 185,
"Landsat 9": 185,
"Terra": 2330,
"Aqua": 2330,
"Suomi-NPP": 3060,
"NOAA-20": 3060,
"NOAA-21": 3060,
"Metop-B": 2900,
"Metop-C": 2900,
"GCOM-C": 1150,
"GCOM-W": 1450,
"ALOS-2": 350,
"GOSAT-2": 920,
"SMAP": 1000,
"PACE": 2663,
"OCO-2": 10,
"ICESat-2": 6,
"GPM Core": 885,
"Sentinel-6": 30,
}
DEFAULT_SWATH_KM = 300
#: How much of the recent past each trail keeps, in minutes, and how finely it
#: is sampled. A whole orbit is about 100 minutes, so this is roughly a third of
#: a circuit -- enough to read the track's curve without wrapping the globe.
TRAIL_MINUTES = 14.0
TRAIL_SAMPLES = 72
TRAIL_PEAK_ALPHA = 0.52
def rotate_about(vecs, axis, ang):
"""Rotate row vectors about `axis` by `ang` radians (Rodrigues' formula).
Args:
vecs: An `(n, 3)` array of vectors.
axis: A unit vector to rotate around.
ang: An `(n,)` array of angles, one per row of `vecs`.
Returns:
numpy.ndarray: The rotated `(n, 3)` array.
"""
c, s = np.cos(ang)[:, None], np.sin(ang)[:, None]
dot = (vecs @ axis)[:, None]
return vecs * c + np.cross(axis, vecs) * s + axis * dot * (1.0 - c)
#: Largest across-track step drawn as a straight chord, in radians. A chord
#: across an angle `a` sags below the surface by 1 - cos(a/2), so at 0.06 rad
#: that is 25 km on a 6371 km sphere -- under a screen pixel here. A 3060 km
#: swath spans 0.48 rad in one piece and sags 170 km, which draws the band as a
#: plate cutting into the planet rather than a strip lying on it.
MAX_CHORD_RAD = 0.06
def _across_steps(half):
"""How many spans an across-track half-angle should be cut into.
Args:
half: Half the swath's angular width, in radians.
Returns:
int: Span count, at least 1. Narrow swaths stay a single quad, so the
subdivision costs nothing where it buys nothing.
"""
return max(1, int(np.ceil(2.0 * half / MAX_CHORD_RAD)))
def ground_swath(track, t_now):
"""The strip of ground `track` has swept over the last `TRAIL_MINUTES`.
The sub-satellite point is the satellite's direction from Earth's centre. Each
sample is then carried around `SPIN_AXIS` by however far Earth has turned since
it was laid down, so the strip stays fixed to the surface rather than sliding
with the satellite. Samples on the far side of the globe are dropped, which is
what keeps the strip from painting through Earth.
Args:
track: A row from `TRACKS`.
t_now: Hours since the start of the clip.
Returns:
tuple: `(quads, alphas)` for a `Poly3DCollection`, or `None` before the
trail has any length.
"""
t0 = t_now - TRAIL_MINUTES / 60.0
if t_now - max(t0, 0.0) < 1e-4:
return None
ts = np.linspace(max(t0, 0.0), t_now, TRAIL_SAMPLES)
spins = 2 * np.pi * ts / 23.9344696
pts = np.array([position(track, t, s) for t, s in zip(ts, spins)])
d = pts / np.linalg.norm(pts, axis=1)[:, None]
d = rotate_about(d, SPIN_AXIS, spins[-1] - spins)
# Cross-track direction: perpendicular to both the ground point and the
# direction of travel along the track.
step = np.gradient(d, axis=0)
side = np.cross(d, step)
norm = np.linalg.norm(side, axis=1)[:, None]
side = np.divide(side, norm, out=np.zeros_like(side), where=norm > 1e-12)
half = SWATH_KM.get(track["name"], DEFAULT_SWATH_KM) / 2.0 / R_EARTH_KM
n_across = _across_steps(half)
offsets = np.linspace(-half, half, n_across + 1)
edges = [
(d * np.cos(a) + side * np.sin(a)) * 1.002 for a in offsets
] # each follows the sphere
front = d @ EYE > 0.02 # drop the far side, and the limb where it degenerates
age = np.linspace(0.0, 1.0, TRAIL_SAMPLES) ** 1.1 * TRAIL_PEAK_ALPHA
quads, alphas = [], []
for k in range(TRAIL_SAMPLES - 1):
if not (front[k] and front[k + 1]):
continue
for m in range(n_across):
quads.append(
[edges[m][k], edges[m + 1][k], edges[m + 1][k + 1], edges[m][k + 1]]
)
alphas.append(age[k])
if not quads:
return None
return np.array(quads), np.array(alphas)
#: How wide the top of a beam is, as a fraction of its footprint width, and how
#: opaque the sheet is. A push-broom instrument images a line across its track, so
#: its field of view is a flat fan rather than a cone: nearly a line at the
#: spacecraft, opening out to the full swath at the ground. Giving the top a small
#: but non-zero width makes it read as a trapezoid instead of a spike.
BEAM_TOP_FRAC = 0.07
BEAM_ALPHA = 0.15
def sensor_beam(track, t_hours, spin):
"""The trapezoidal field of view between a satellite and the ground it is imaging.
Args:
track: A row from `TRACKS`.
t_hours: Hours since the start of the clip.
spin: Earth's rotation angle at that moment, in radians.
Returns:
numpy.ndarray: The trapezoid's four corners, or `None` when the footprint
is on the far side of the globe and the sheet would draw through Earth.
"""
dt = 1.0 / 3600.0 # one second on
spin_ahead = spin + 2 * np.pi * dt / 23.9344696
here = position(track, t_hours, spin)
ahead = position(track, t_hours + dt, spin_ahead)
down = here / np.linalg.norm(here)
if float(down @ EYE) <= 0.02:
return None
# Both samples carried into the same Earth-fixed frame before differencing.
# Taking `ahead` at the current spin measured the inertial velocity, while
# ground_swath measures the ground track -- so the beam's edge and the
# leading edge of the trail it paints were skewed by a few degrees against
# each other, which is exactly the angle a beam should not be wrong by.
d_ahead = ahead / np.linalg.norm(ahead)
d_ahead = rotate_about(d_ahead[None, :], SPIN_AXIS, np.array([spin - spin_ahead]))[
0
]
across = np.cross(down, d_ahead - down)
norm = np.linalg.norm(across)
if norm < 1e-12:
return None
across /= norm
half = SWATH_KM.get(track["name"], DEFAULT_SWATH_KM) / 2.0 / R_EARTH_KM
# The ground edge is an arc, cut finely enough that it lies on the surface
# instead of chording under it.
n_across = _across_steps(half)
arc = np.array(
[
(down * np.cos(a) + across * np.sin(a)) * 1.004
for a in np.linspace(half, -half, n_across + 1)
]
)
top = across * (BEAM_TOP_FRAC * half)
return np.vstack([[here - top], [here + top], arc])
#: Depth band over which a low-orbit label fades in, as a fraction of the orbit
#: radius toward the camera. Below the first value the satellite is on the far
#: side and unlabelled; above the second the label is at full strength.
NEAR_LABEL_FADE = (0.10, 0.45)
HOURS = 1.5 # Everything runs at its true rate; the window is short instead.
# A low orbiter covers a little under one circuit here, which is slow enough to
# follow, and the ratio between orbit rate and Earth's spin stays correct -- so
# the ground tracks are real. Slowing the satellites instead would have shown
# the same on-screen speed (0.25 x 6.0 h = 1.0 x 1.5 h) but drawn swaths that
# drift west four times too fast.
N_FRAMES = 480
CLIP_FPS = 24
#: How much faster than real time the clip runs -- the simulated window divided by
#: the screen time. It is stated on the clock because a bare elapsed time does not
#: tell a viewer what to do with it: the point is that a low orbiter completing one
#: circuit as the clock reaches 1.5 h is showing them its orbital period.
SPEEDUP = round(HOURS * 3600 / (N_FRAMES / CLIP_FPS))
#: RENDER_DPI 160 on the 12x8 inch figure gives a 1920x1280 master. The bitrate
#: has to rise with the resolution: an mp4's size is bitrate x duration no matter
#: how many pixels are in it, so at the 3600 default a 1920-wide clip simply
#: compresses away the mesh detail the render spent hours producing.
RENDER_DPI = 160
RENDER_BITRATE = 12000
ELEV, AZIM = 22.0, -58.0
#: Unit vector from the origin toward the camera. Projecting a point onto it
#: gives how far that point lies toward the viewer: +1 straight in front of
#: Earth, 0 on the limb, -1 straight behind. Both the ring fade and the
#: near-side label test are this projection, so they agree with what is
#: actually in front of what.
EYE = np.array(
[
np.cos(np.deg2rad(ELEV)) * np.cos(np.deg2rad(AZIM)),
np.cos(np.deg2rad(ELEV)) * np.sin(np.deg2rad(AZIM)),
np.sin(np.deg2rad(ELEV)),
]
)
#: The painting order. Matplotlib's own 3D depth sorting is off
#: (computed_zorder = False), so every layer's place is stated here rather than
#: left to the order the calls happen to be made in. The globe sits in the
#: middle: anything on the far side goes below it, anything nearer goes above.
Z_BEHIND_GLOBE = 2
Z_GLOBE = 4
Z_NEAR_RING = 4.5
Z_SURFACE = 5 # swaths and beams, lying on the near face
Z_NEAR_SAT = 6
Z_LABEL = 7
#: The geostationary ring is drawn once, in a tone of its own: it is shared by
#: seven satellites, so no single provider's colour is the honest one for it.
GEO_RING_COLOUR = "#93a4b8"
LIM = 3.35 # tight around the drawn geostationary ring
fig = plt.figure(figsize=(12, 8), facecolor="black")
# Fitted so the scene fills 93% x 83% of the frame, centred, with the
# geostationary labels just inside the edges rather than clipped. The
# projected scene is a wide, shallow ellipse, so the axes overflows the
# figure on every side to leave no dead margin.
ax = fig.add_axes([-0.345, -0.437, 1.665, 1.873], projection="3d")
# mplot3d sorts artists by mean depth, which lets the globe paint over orbits
# in front of it; honour the zorder the drawing assigns instead.
ax.computed_zorder = False
def draw_frame(i):
"""Draw one frame: Earth, every orbit, and every satellite on it.
Clears and redraws the shared axis in place rather than building a new one,
so the figure the writer is grabbing from stays the same object throughout
the render.
Args:
i: Frame index, from 0 to `N_FRAMES - 1`. Time and Earth's rotation are
both derived from it, so the frame is fully determined by this.
Returns:
tuple: Empty. Matplotlib's animation callbacks return the artists to
re-blit; this redraws the whole axis instead, so there are none.
"""
ax.clear()
ax.set_axis_off()
ax.set_facecolor("black")
ax.set_xlim(-LIM, LIM)
ax.set_ylim(-LIM, LIM)
ax.set_zlim(-LIM * 0.62, LIM * 0.62)
ax.set_box_aspect((1, 1, 0.62))
ax.view_init(elev=ELEV, azim=AZIM)
t_hours = HOURS * i / (N_FRAMES - 1)
spin = 2 * np.pi * t_hours / 23.9344696 # one sidereal day
textured_earth(ax, spin)
# The orbit paths, faded continuously by depth so the far side of each ring
# recedes instead of competing with the near side. 34 rings at similar radii
# otherwise read as one flat tangle. A per-segment alpha ramp is used rather
# than a near/far split, which leaves a hard seam where the two halves meet.
#
# Each ring is split by depth against the globe's own zorder, though: a ring
# is a closed loop with half of it behind Earth, so drawing the whole thing
# at one zorder is wrong either way round. Held entirely behind, the near
# half vanished into the planet it should pass in front of.
geo_ring_drawn = False
for track in TRACKS:
# plane_basis(0, lon) rotates one equatorial circle, so every
# geostationary row traces identical geometry. Drawing all seven stacked
# them into a line about seven times its own alpha -- the brightest thing
# in the frame, and colour-neutral from the overlap rather than any
# provider's colour. One ring, drawn once, in a neutral tone.
if track["kind"] == "geo":
if geo_ring_drawn:
continue
geo_ring_drawn = True
colour = GEO_RING_COLOUR if track["kind"] == "geo" else track["colour"]
ox, oy, oz = orbit_ring(track)
pts = np.stack([ox, oy, oz], axis=1)
segs = np.stack([pts[:-1], pts[1:]], axis=1)
depth = ox * EYE[0] + oy * EYE[1] + oz * EYE[2]
depth = (depth[:-1] + depth[1:]) / 2 / track["r_draw"] # -1 far .. +1 near
rgba = np.tile(to_rgba(colour), (len(segs), 1))
rgba[:, 3] = 0.045 + 0.20 * (depth + 1) / 2
near = depth > 0
for mask, z in ((~near, Z_BEHIND_GLOBE), (near, Z_NEAR_RING)):
if mask.any():
ax.add_collection3d(
Line3DCollection(
segs[mask], colors=rgba[mask], linewidths=0.20, zorder=z
)
)
# The instrument's field of view: a flat sheet from the spacecraft down to the
# strip it is imaging right now. Drawn before the trail so the trail's older,
# fainter ground reads as behind it.
for track in LOW_ORBIT_TRACKS:
quad = sensor_beam(track, t_hours, spin)
if quad is None:
continue
ax.add_collection3d(
Poly3DCollection(
[quad],
facecolors=[to_rgba(track["colour"], BEAM_ALPHA)],
edgecolors="none",
zorder=Z_SURFACE,
)
)
# The ground each low orbiter has swept recently, fading out behind it. This
# is the point of a sun-synchronous orbit: a narrow strip at a time, a new one
# every pass. Geostationary satellites are skipped -- they stare at one face.
for track in LOW_ORBIT_TRACKS:
swept = ground_swath(track, t_hours)
if swept is None:
continue
quads, alphas = swept
rgba = np.tile(to_rgba(track["colour"]), (len(quads), 1))
rgba[:, 3] = alphas
ax.add_collection3d(
Poly3DCollection(
quads, facecolors=rgba, edgecolors="none", zorder=Z_SURFACE
)
)
# A satellite on the far side belongs behind the globe. Painted above it,
# roughly a third of the fleet floated over the disc at any moment, which
# reads as a flat sticker rather than an orbit. Held below the globe's
# zorder it is hidden only where Earth actually covers it -- one that is
# behind but outside the limb still shows, which is the truth of it.
for track in TRACKS:
p = position(track, t_hours, spin)
ax.scatter(
[p[0]],
[p[1]],
[p[2]],
marker=SATELLITE,
s=70,
c=track["colour"],
edgecolors="none",
depthshade=False,
zorder=Z_NEAR_SAT if float(p @ EYE) > 0 else Z_BEHIND_GLOBE,
)
# Name the low-orbit satellites on the near side only. Labelling all 27
# doubles the text on screen and the far half is already dimmed, so the
# names would sit on satellites the eye is not following. A satellite fades
# its label in as it comes round, rather than popping, so the set of names
# changes smoothly instead of flickering.
for track in LOW_ORBIT_TRACKS:
p = position(track, t_hours, spin)
depth = (p[0] * EYE[0] + p[1] * EYE[1] + p[2] * EYE[2]) / track["r_draw"]
if depth <= NEAR_LABEL_FADE[0]:
continue
a = np.clip(
(depth - NEAR_LABEL_FADE[0]) / (NEAR_LABEL_FADE[1] - NEAR_LABEL_FADE[0]),
0.0,
1.0,
)
ax.text(
p[0],
p[1],
p[2] + 0.085,
track["name"],
color=track["colour"],
fontsize=7.0,
ha="center",
va="bottom",
alpha=float(a),
zorder=Z_LABEL,
path_effects=[
pe.withStroke(linewidth=2.0, foreground="black", alpha=float(a))
],
)
# Name the geostationary satellites. Meteosat-12/-9/-11 sit within 9.5 deg
# of each other; LABEL_OFFSET is the hook for pushing their labels to
# of overlapping on the ring.
for track in GEO_TRACKS:
p = position(track, t_hours, spin)
out = LABEL_OFFSET.get(track["name"], 1.10)
ax.text(
p[0] * out,
p[1] * out,
p[2] * out + 0.06,
track["name"],
color=track["colour"],
fontsize=8.5,
ha="center",
va="center",
zorder=Z_LABEL,
path_effects=[pe.withStroke(linewidth=2.2, foreground="black")],
)
ax.text2D(
0.022,
0.035,
f"{t_hours:4.1f} h simulated · {SPEEDUP}× real time",
transform=fig.transFigure,
color="#9fb4c7",
fontsize=13,
ha="left",
family="monospace",
)
return ()
stamp_logo(fig)
draw_frame(0)
fig
What to notice¶
- The geostationary ring does not move. Three GOES, three Meteosat and Himawari-9 sit motionless above their
longitudes while Earth turns underneath at the same rate.
positiontakes the shortcut here and carries them with the surface rotation instead of integrating their orbit — it is entitled to, because Kepler's period at that radius comes out within half a minute of the sidereal day, which the check above prints. The ring is parked by construction; what Kepler establishes is that construction is honest. - The low-Earth shell is retrograde. Sun-synchronous orbits are inclined just past 90°, so they cross the poles going backwards relative to Earth's spin. That is what keeps their ground track crossing the equator at the same local solar time every day.
- ICESat-2 is the odd one out among the near-polar satellites at 92°, and GPM is odder still at 65° — it was deliberately given a low inclination to sample the tropics at every hour of the day, which a sun-synchronous orbit can never do.
- The gap is the story. Almost everything sits in a thin shell close to the surface; a handful sit six times further out. There is very little in between, and that is a real feature of how Earth observation is done — the radial compression above squeezes that gap to fit the frame, but it cannot invent one that is not there.
In [ ]:
Copied!
# The clips are published from docs/_images/, which the site serves and git
# tracks. The notebook's own out/ is ignored and excluded from the build, so a
# clip left there could not be referenced by the README or the documentation.
#
# The -1920 master is the exception: it is gitignored. Every published size is
# transcoded from it, so nothing ever references the master itself, and a
# re-render writes a byte-different 27 MB copy that video codecs cannot
# delta-compress. Committing it would cost its full size again every time.
OUT_DIR = Path("../../_images/animation")
OUT_DIR.mkdir(parents=True, exist_ok=True)
# Every published basemap in one pass. render_master draws every frame once per
# clip; publish_clip derives the rest by transcoding, so only the master is
# expensive. Running the whole tuple is what makes the two published sets
# comparable -- they come out of one state of this notebook, rather than one
# being whatever the last run with that constant set happened to produce.
written = []
for name in PUBLISHED_BASEMAPS:
select_basemap(name)
out = OUT_DIR / f"earthlens-satellites-{name}-1920"
master = render_master(
fig,
draw_frame,
N_FRAMES,
str(out),
fps=CLIP_FPS,
bitrate=RENDER_BITRATE,
dpi=RENDER_DPI,
)
web, gif, webp = publish_clip(master)
written.append((name, gif, [Path(master), Path(web), Path(gif), Path(webp)]))
plt.close("all")
for name, _, files in written:
print(name)
for f in files:
print(f" {f.name:<44} {f.stat().st_size / 1e6:6.1f} MB")
display(*[Image(filename=gif) for _, gif, _ in written])
# The clips are published from docs/_images/, which the site serves and git
# tracks. The notebook's own out/ is ignored and excluded from the build, so a
# clip left there could not be referenced by the README or the documentation.
#
# The -1920 master is the exception: it is gitignored. Every published size is
# transcoded from it, so nothing ever references the master itself, and a
# re-render writes a byte-different 27 MB copy that video codecs cannot
# delta-compress. Committing it would cost its full size again every time.
OUT_DIR = Path("../../_images/animation")
OUT_DIR.mkdir(parents=True, exist_ok=True)
# Every published basemap in one pass. render_master draws every frame once per
# clip; publish_clip derives the rest by transcoding, so only the master is
# expensive. Running the whole tuple is what makes the two published sets
# comparable -- they come out of one state of this notebook, rather than one
# being whatever the last run with that constant set happened to produce.
written = []
for name in PUBLISHED_BASEMAPS:
select_basemap(name)
out = OUT_DIR / f"earthlens-satellites-{name}-1920"
master = render_master(
fig,
draw_frame,
N_FRAMES,
str(out),
fps=CLIP_FPS,
bitrate=RENDER_BITRATE,
dpi=RENDER_DPI,
)
web, gif, webp = publish_clip(master)
written.append((name, gif, [Path(master), Path(web), Path(gif), Path(webp)]))
plt.close("all")
for name, _, files in written:
print(name)
for f in files:
print(f" {f.name:<44} {f.stat().st_size / 1e6:6.1f} MB")
display(*[Image(filename=gif) for _, gif, _ in written])
Sources for the orbital elements¶
- GOES positions (GOES-19 at 75.2°W, GOES-18 at 137.0°W, GOES-16 at 104.7°W): NOAA NESDIS
- Meteosat positions (Meteosat-12 prime at 0°, Meteosat-11 rapid-scan at 9.5°E): EUMETSAT Meteosat series, MSG services
- Meteosat-9 (45.5°E, Indian Ocean Data Coverage since 1 June 2022): EUMETSAT IODC
- Himawari-9 (140.7°E): eoPortal — Himawari-8/9
- Sentinel-1 (693 km, 98.18°): Sentinel-1 orbit
- Sentinel-2 (786 km, 98.62°): Sentinel-2 orbit
- Sentinel-3 (815 km, 98.6°): eoPortal — Sentinel-3
- Sentinel-5P (824 km, 98.74°): eoPortal — Sentinel-5P
- Sentinel-6 (1336 km, 66°): eoPortal — Sentinel-6
- Landsat 8/9 (705 km, 98.2°): USGS Landsat 8, NASA Landsat 9
- Terra (694 km, lowered by two retrograde manoeuvres in October 2022) and Aqua (705 km): NSIDC — ongoing changes in Terra and Aqua orbits
- NOAA-20 (824 km, 98.79°): NOAA-20
- NOAA-21 (833 km, 98.80°): NOAA-21
- Suomi-NPP (824 km, 98.74°): eoPortal — Suomi NPP
- OCO-2 (705 km, 98.2°): eoPortal — OCO-2
- ICESat-2 (496 km, 92.0°, deliberately not sun-synchronous): ICESat-2 technical specs (NASA states ~500 km)
- GPM Core (407 km, 65°): NASA GPM Core Observatory
- PACE (676.5 km, 98°): NASA SVS — PACE orbit
- Metop (817 km): EUMETSAT Metop series
- GCOM-C (798 km) / GCOM-W (700 km): eoPortal — GCOM-C
- ALOS-2 (628 km, 97.9°): eoPortal — ALOS-2
- GOSAT-2 (613 km, 97.84°): eoPortal — GOSAT-2
- SMAP (685 km, 98°): eoPortal — SMAP