Sub-daily hurricanes in true colour: Ida (2021) & Milton (2024)¶

The polar-orbiting satellites in the other showcases fly over a spot only once or twice a day — fine for a yearly time-lapse, useless for watching a hurricane spin. For that you need a geostationary satellite, parked over one longitude and imaging the same disc every few minutes. Over the Americas that is NOAA's GOES-R series: its Advanced Baseline Imager (ABI) scans the continental US every 5 minutes. This notebook builds two full-life-cycle true-colour time-lapses from the earthlens goes backend — following each storm from before it even has a whirl until it falls apart: Hurricane Ida (Aug 2021) from a ragged Caribbean disturbance, through its Louisiana landfall, to its decaying remnants over the eastern US; and Hurricane Milton (Oct 2024) from a broad Gulf low, through its explosive run at Florida, to its exit into the Atlantic — embedded as animations.

Two things make GOES different from the Landsat/MODIS showcases:

• True colour needs a synthesised green. ABI has red (C02), blue (C01) and a "veggie" near-infrared (C03) band, but no green — so green is reconstructed with the standard CIMSS weights 0.45·red + 0.10·veggie + 0.45·blue. The granules also arrive in the satellite's geostationary projection, so each frame is reprojected to lat/lon. Both steps are pure GIS plumbing that really belongs in pyramids; here they live in one small render helper (a good pyramids port candidate).
• Visible light only works in daylight. To follow each storm over several days — from a disorganised cloud cluster to a major hurricane at the coast — we sample a handful of daytime hours (UTC) on each day and let the nights fall out; the storm's growth and track carry across the day-to-day jumps.

Setup¶

The goes backend reads raw ABI NetCDF granules from NOAA's public S3 buckets — anonymous, no credentials. pyramids builds and animates the frames; the render helper below reads each band and reprojects it from the satellite's geostationary grid to lat/lon with pyramids.netcdf.NetCDF directly (no raw driver calls). The MP4 export uses ffmpeg, which cleopatra (≥0.22) finds on the PATH or falls back to the copy bundled with imageio-ffmpeg — so no separate install is needed.

In [ ]:

import base64
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from cleopatra.glyphs.gridded.array_glyph import FrameLabel
from IPython.display import HTML
from pyramids.dataset import Dataset
from pyramids.dataset.collection import DatasetCollection
from pyramids.netcdf import NetCDF

from earthlens.core import EarthLens

OUT = Path("out") / "hurricanes"
OUT.mkdir(parents=True, exist_ok=True)

import base64 from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd from cleopatra.glyphs.gridded.array_glyph import FrameLabel from IPython.display import HTML from pyramids.dataset import Dataset from pyramids.dataset.collection import DatasetCollection from pyramids.netcdf import NetCDF from earthlens.core import EarthLens OUT = Path("out") / "hurricanes" OUT.mkdir(parents=True, exist_ok=True)

From a raw GOES granule to a true-colour frame¶

A GOES MCMIP granule is a stack of 16 ABI bands in the satellite's geostationary grid. To make one true-colour frame we pull the three visible-ish bands, reproject each to lat/lon over the storm box, synthesise the missing green, gamma-stretch, and write an ordinary RGB GeoTIFF — which then plugs straight into the same DatasetCollection animation the other showcases use.

pyramids.netcdf.NetCDF reads a band, applies its CF scale_factor/add_offset (read_array(unpack=True)), and reprojects the geostationary grid to lat/lon (to_crs(4326)) with no manual CRS wrangling. The one rough edge left: a windowed read_array(bbox=...) on a to_crs()'d geostationary variable currently swaps width/height (pyramids#719, filed after the crash it used to hit was fixed upstream) — so we read the full reprojected band and crop it with plain numpy slicing instead.

In [ ]:

def goes_band_to_latlon(nc_path, band, bbox):
    """Read one ABI CMI band from a granule and reproject it to lat/lon over bbox."""
    variable = NetCDF.read_file(nc_path).get_variable(band)
    reprojected = variable.to_crs(4326)
    full = np.squeeze(np.asarray(reprojected.read_array(unpack=True), dtype="float32"))

    # Crop to bbox by hand: read_array(bbox=...) on a to_crs()'d geostationary
    # variable currently swaps width/height (pyramids#719).
    west, south, east, north = bbox
    origin_x, pixel_w, _, origin_y, _, pixel_h = reprojected.geotransform
    col0, col1 = round((west - origin_x) / pixel_w), round((east - origin_x) / pixel_w)
    row0, row1 = (
        round((north - origin_y) / pixel_h),
        round((south - origin_y) / pixel_h),
    )
    cropped = full[row0:row1, col0:col1]
    geo = (west, pixel_w, 0.0, north, 0.0, pixel_h)
    return cropped, geo


def render_true_colour(nc_path, dst, bbox):
    """Write a gamma-stretched true-colour RGB GeoTIFF from one MCMIP granule."""
    red, geo = goes_band_to_latlon(nc_path, "CMI_C02", bbox)  # 0.64 um red
    veggie, _ = goes_band_to_latlon(nc_path, "CMI_C03", bbox)  # 0.86 um near-IR
    blue, _ = goes_band_to_latlon(nc_path, "CMI_C01", bbox)  # 0.47 um blue
    green = 0.45 * red + 0.10 * veggie + 0.45 * blue  # CIMSS synthetic green
    rgb = np.power(np.clip([red, green, blue], 0.0, 1.0), 1 / 2.2)
    rgb = np.nan_to_num(rgb * 255.0).astype("uint8")
    scene = Dataset.create_from_array(arr=rgb, geo=geo, epsg=4326, no_data_value=0)
    scene.to_file(str(dst))

def goes_band_to_latlon(nc_path, band, bbox): """Read one ABI CMI band from a granule and reproject it to lat/lon over bbox.""" variable = NetCDF.read_file(nc_path).get_variable(band) reprojected = variable.to_crs(4326) full = np.squeeze(np.asarray(reprojected.read_array(unpack=True), dtype="float32")) # Crop to bbox by hand: read_array(bbox=...) on a to_crs()'d geostationary # variable currently swaps width/height (pyramids#719). west, south, east, north = bbox origin_x, pixel_w, _, origin_y, _, pixel_h = reprojected.geotransform col0, col1 = round((west - origin_x) / pixel_w), round((east - origin_x) / pixel_w) row0, row1 = ( round((north - origin_y) / pixel_h), round((south - origin_y) / pixel_h), ) cropped = full[row0:row1, col0:col1] geo = (west, pixel_w, 0.0, north, 0.0, pixel_h) return cropped, geo def render_true_colour(nc_path, dst, bbox): """Write a gamma-stretched true-colour RGB GeoTIFF from one MCMIP granule.""" red, geo = goes_band_to_latlon(nc_path, "CMI_C02", bbox) # 0.64 um red veggie, _ = goes_band_to_latlon(nc_path, "CMI_C03", bbox) # 0.86 um near-IR blue, _ = goes_band_to_latlon(nc_path, "CMI_C01", bbox) # 0.47 um blue green = 0.45 * red + 0.10 * veggie + 0.45 * blue # CIMSS synthetic green rgb = np.power(np.clip([red, green, blue], 0.0, 1.0), 1 / 2.2) rgb = np.nan_to_num(rgb * 255.0).astype("uint8") scene = Dataset.create_from_array(arr=rgb, geo=geo, epsg=4326, no_data_value=0) scene.to_file(str(dst))

Fetch the daytime frames and animate each storm¶

For each storm we walk the days of its life and, on each day, a few daytime hours (satellite="16" — GOES-16 was GOES-East for both dates); for every slot we pull the one ABI CONUS granule, render it to a true-colour GeoTIFF, and stamp the frame with its date. Then we animate the stack with pyramids: collection.plot(rgb=…) → glyph.animate(...) with the time label lifted into the top margin, and write the result as a crisp, light MP4 with glyph.save_animation(..., crf=26) — the right format for a LinkedIn post, where it uploads as a native video (the crf keeps the file small; cleopatra pads odd dimensions and forces yuv420p for you). Each panel is sized to its own bounding-box aspect so the map fills the frame with barely any white margin. To give it the clean look of an operational weather-centre chart, one call — glyph.add_reference_map(style="ecmwf-dark", extent=…) (cleopatra ≥0.22) — draws grey coastlines + country borders, a dashed lon/lat graticule and °W/°N labels, georeferenced to the storm's box, so it is unmistakable which coastline the storm is bearing down on. It is called after animate() — calling it first triggers a cleopatra blitting bug where the saved video's image never advances past frame 0 (see the code comment below).

In [ ]:

# Each storm is followed from a bit before it organises through to landfall.
# Geostationary visible bands only work in daylight, so we sample a few daytime
# hours (UTC) on each day and let the night gaps fall out — the storm's growth and
# track carry across the day-to-day jumps. `bbox` is wide enough to hold the whole
# track (Ida spins up in the Caribbean, Milton in the south-western Gulf). The
# saved MP4 is trimmed to a tight border with an ffmpeg crop that is computed
# from the axes' actual rendered position each run (below), so it adapts to
# whatever margins the current cleopatra layout produces.
STORMS = [
    {
        "name": "Ida",
        "satellite": "16",
        "days": ("2021-08-24", "2021-09-02"),
        "hours": [14, 16, 18, 20],
        "bbox": [-99.0, 15.0, -66.0, 46.0],
    },
    {
        "name": "Milton",
        "satellite": "16",
        "days": ("2024-10-04", "2024-10-11"),
        "hours": [14, 16, 18, 20],
        "bbox": [-98.0, 16.0, -64.0, 36.0],
    },
]
DPI = 200

videos = {}
for storm in STORMS:
    days = pd.date_range(storm["days"][0], storm["days"][1], freq="1D")
    slots = [day + pd.Timedelta(hours=hour) for day in days for hour in storm["hours"]]
    out_dir = OUT / storm["name"].lower()
    out_dir.mkdir(parents=True, exist_ok=True)
    labels, paths = [], []
    for slot in slots:
        job = EarthLens(
            data_source="goes",
            dataset="abi-l2-mcmip",
            satellite=storm["satellite"],
            domain="C",
            start=slot.strftime("%Y-%m-%d %H:%M"),
            end=(slot + pd.Timedelta(minutes=6)).strftime("%Y-%m-%d %H:%M"),
            fmt="%Y-%m-%d %H:%M",
            lat_lim=[storm["bbox"][1], storm["bbox"][3]],
            lon_lim=[storm["bbox"][0], storm["bbox"][2]],
            path=out_dir / "granules",
        )
        got = job.download(progress_bar=False)
        if not got:
            continue  # occasional GOES scan gap at this exact slot — skip the frame
        frame = out_dir / f"{slot.strftime('%Y%m%d_%H%M')}.tif"
        render_true_colour(got[0], frame, storm["bbox"])
        labels.append(slot.strftime("%b%d %HZ"))
        paths.append(frame)

    # Size each frame to its own bounding-box aspect (Milton is much wider than
    # tall) at a LinkedIn-friendly ~1100 px width, so the map fills the panel with
    # barely any white margin.
    west, south, east, north = storm["bbox"]
    fig_w = 11.0
    figsize = (fig_w, fig_w * (north - south) / (east - west) / 0.905 + 0.32)

    cube = DatasetCollection.from_files(paths)
    glyph = cube.plot(
        rgb_options={"rgb": [0, 1, 2], "surface_reflectance": 255}, figsize=figsize
    )

    # Float the per-frame scan-time label above the frame instead of riding on the
    # map. An explicit FrameLabel location switches to data coordinates and is
    # already non-clipping; ax.transData is live, so setting it here resolves
    # correctly once drawn — no transform/clip_on/subplots_adjust juggling needed.
    #
    # Called BEFORE add_reference_map() (cleopatra bug: with blit=True,
    # add_reference_map()'s basemap artists drawn before animate()'s own
    # ax.imshow() poison FuncAnimation's cached blit background, so the image
    # never updates past frame 0 in the saved video while the frame label still
    # advances -- reproduced against both cleopatra 0.25.0 and 0.26.0, filed
    # upstream. Swapping the call order (animate() first) avoids it.
    label_location = [west, north + 0.02 * (north - south)]
    glyph.animate(
        labels,
        interval=250,
        frame_label=FrameLabel(location=label_location),
    )

    # Dress the frames as an ECMWF-style reference map in one call (cleopatra >=0.22):
    # grey Natural Earth coastlines + borders, a dashed lon/lat graticule, °W/°N
    # labels and a subtle frame. `"ecmwf-dark"` keeps the greys light so they read on
    # the dark satellite background; `extent=` georeferences the RGB/animate axes
    # (which are otherwise in pixel coordinates). Called after animate() -- see the
    # blit note above.
    glyph.add_reference_map(style="dark", extent=[west, south, east, north])

    # cleopatra's internal tight_layout() (inside animate()) only tightens, it
    # doesn't centre -- the y-axis °N tick labels eat space on the left with
    # nothing matching on the right, so the map ends up hard against the right
    # edge. Force equal left/right margins (measured empirically) so the panel
    # sits centred instead. `storm["crop"]` then trims that (now-equal) margin
    # down to a thin, LinkedIn-tight border on all four sides.
    glyph.fig.subplots_adjust(left=0.121, right=0.879, bottom=0.05, top=0.92)

    # Compute the ffmpeg crop rectangle from the axes' ACTUAL rendered position
    # each run, so it tracks whatever margins the current cleopatra layout and
    # reference-map style produce (the previous hand-measured crop broke whenever
    # a cleopatra release shifted the panel). Pad a few px, clamp to the frame,
    # and keep both dimensions even for the yuv420p encoder.
    frame_w, frame_h = (glyph.fig.get_size_inches() * DPI).round().astype(int)
    ax_box = glyph.ax.get_position()  # figure fractions, y=0 at the bottom
    ax_left, ax_right = ax_box.x0 * frame_w, ax_box.x1 * frame_w
    ax_top, ax_bottom = (1 - ax_box.y1) * frame_h, (1 - ax_box.y0) * frame_h
    pad = 8
    crop_x = max(int(ax_left) - pad, 0)
    crop_y = max(int(ax_top) - pad, 0)
    crop_w = min(int(round(ax_right - ax_left)) + 2 * pad, frame_w - crop_x)
    crop_h = min(int(round(ax_bottom - ax_top)) + 2 * pad, frame_h - crop_y)
    crop_w -= crop_w % 2
    crop_h -= crop_h % 2
    crop = f"crop={crop_w}:{crop_h}:{crop_x}:{crop_y}"

    # Write a crisp, LinkedIn-ready MP4: 200 dpi (roughly 2200 px wide) for a
    # sharp look on large feeds, crf=20 (visually near-lossless, still a small
    # file), preset="slow" trades a bit of encode time for better compression at
    # that quality, and the crop filter trims the figure margins down to a tight
    # ~30 px border so the map fills the frame edge-to-edge instead of sitting in
    # a wide grey mat.
    mp4_path = out_dir / f"{storm['name'].lower()}.mp4"
    glyph.save_animation(
        str(mp4_path),
        fps=2,
        dpi=DPI,
        crf=20,
        preset="slow",
        extra_args=["-vf", crop],
    )
    plt.close("all")
    videos[storm["name"]] = mp4_path
list(videos)

# Each storm is followed from a bit before it organises through to landfall. # Geostationary visible bands only work in daylight, so we sample a few daytime # hours (UTC) on each day and let the night gaps fall out — the storm's growth and # track carry across the day-to-day jumps. `bbox` is wide enough to hold the whole # track (Ida spins up in the Caribbean, Milton in the south-western Gulf). The # saved MP4 is trimmed to a tight border with an ffmpeg crop that is computed # from the axes' actual rendered position each run (below), so it adapts to # whatever margins the current cleopatra layout produces. STORMS = [ { "name": "Ida", "satellite": "16", "days": ("2021-08-24", "2021-09-02"), "hours": [14, 16, 18, 20], "bbox": [-99.0, 15.0, -66.0, 46.0], }, { "name": "Milton", "satellite": "16", "days": ("2024-10-04", "2024-10-11"), "hours": [14, 16, 18, 20], "bbox": [-98.0, 16.0, -64.0, 36.0], }, ] DPI = 200 videos = {} for storm in STORMS: days = pd.date_range(storm["days"][0], storm["days"][1], freq="1D") slots = [day + pd.Timedelta(hours=hour) for day in days for hour in storm["hours"]] out_dir = OUT / storm["name"].lower() out_dir.mkdir(parents=True, exist_ok=True) labels, paths = [], [] for slot in slots: job = EarthLens( data_source="goes", dataset="abi-l2-mcmip", satellite=storm["satellite"], domain="C", start=slot.strftime("%Y-%m-%d %H:%M"), end=(slot + pd.Timedelta(minutes=6)).strftime("%Y-%m-%d %H:%M"), fmt="%Y-%m-%d %H:%M", lat_lim=[storm["bbox"][1], storm["bbox"][3]], lon_lim=[storm["bbox"][0], storm["bbox"][2]], path=out_dir / "granules", ) got = job.download(progress_bar=False) if not got: continue # occasional GOES scan gap at this exact slot — skip the frame frame = out_dir / f"{slot.strftime('%Y%m%d_%H%M')}.tif" render_true_colour(got[0], frame, storm["bbox"]) labels.append(slot.strftime("%b%d %HZ")) paths.append(frame) # Size each frame to its own bounding-box aspect (Milton is much wider than # tall) at a LinkedIn-friendly ~1100 px width, so the map fills the panel with # barely any white margin. west, south, east, north = storm["bbox"] fig_w = 11.0 figsize = (fig_w, fig_w * (north - south) / (east - west) / 0.905 + 0.32) cube = DatasetCollection.from_files(paths) glyph = cube.plot( rgb_options={"rgb": [0, 1, 2], "surface_reflectance": 255}, figsize=figsize ) # Float the per-frame scan-time label above the frame instead of riding on the # map. An explicit FrameLabel location switches to data coordinates and is # already non-clipping; ax.transData is live, so setting it here resolves # correctly once drawn — no transform/clip_on/subplots_adjust juggling needed. # # Called BEFORE add_reference_map() (cleopatra bug: with blit=True, # add_reference_map()'s basemap artists drawn before animate()'s own # ax.imshow() poison FuncAnimation's cached blit background, so the image # never updates past frame 0 in the saved video while the frame label still # advances -- reproduced against both cleopatra 0.25.0 and 0.26.0, filed # upstream. Swapping the call order (animate() first) avoids it. label_location = [west, north + 0.02 * (north - south)] glyph.animate( labels, interval=250, frame_label=FrameLabel(location=label_location), ) # Dress the frames as an ECMWF-style reference map in one call (cleopatra >=0.22): # grey Natural Earth coastlines + borders, a dashed lon/lat graticule, °W/°N # labels and a subtle frame. `"ecmwf-dark"` keeps the greys light so they read on # the dark satellite background; `extent=` georeferences the RGB/animate axes # (which are otherwise in pixel coordinates). Called after animate() -- see the # blit note above. glyph.add_reference_map(style="dark", extent=[west, south, east, north]) # cleopatra's internal tight_layout() (inside animate()) only tightens, it # doesn't centre -- the y-axis °N tick labels eat space on the left with # nothing matching on the right, so the map ends up hard against the right # edge. Force equal left/right margins (measured empirically) so the panel # sits centred instead. `storm["crop"]` then trims that (now-equal) margin # down to a thin, LinkedIn-tight border on all four sides. glyph.fig.subplots_adjust(left=0.121, right=0.879, bottom=0.05, top=0.92) # Compute the ffmpeg crop rectangle from the axes' ACTUAL rendered position # each run, so it tracks whatever margins the current cleopatra layout and # reference-map style produce (the previous hand-measured crop broke whenever # a cleopatra release shifted the panel). Pad a few px, clamp to the frame, # and keep both dimensions even for the yuv420p encoder. frame_w, frame_h = (glyph.fig.get_size_inches() * DPI).round().astype(int) ax_box = glyph.ax.get_position() # figure fractions, y=0 at the bottom ax_left, ax_right = ax_box.x0 * frame_w, ax_box.x1 * frame_w ax_top, ax_bottom = (1 - ax_box.y1) * frame_h, (1 - ax_box.y0) * frame_h pad = 8 crop_x = max(int(ax_left) - pad, 0) crop_y = max(int(ax_top) - pad, 0) crop_w = min(int(round(ax_right - ax_left)) + 2 * pad, frame_w - crop_x) crop_h = min(int(round(ax_bottom - ax_top)) + 2 * pad, frame_h - crop_y) crop_w -= crop_w % 2 crop_h -= crop_h % 2 crop = f"crop={crop_w}:{crop_h}:{crop_x}:{crop_y}" # Write a crisp, LinkedIn-ready MP4: 200 dpi (roughly 2200 px wide) for a # sharp look on large feeds, crf=20 (visually near-lossless, still a small # file), preset="slow" trades a bit of encode time for better compression at # that quality, and the crop filter trims the figure margins down to a tight # ~30 px border so the map fills the frame edge-to-edge instead of sitting in # a wide grey mat. mp4_path = out_dir / f"{storm['name'].lower()}.mp4" glyph.save_animation( str(mp4_path), fps=2, dpi=DPI, crf=20, preset="slow", extra_args=["-vf", crop], ) plt.close("all") videos[storm["name"]] = mp4_path list(videos)

Embed both animations (stacked, large) as base64 MP4 data URIs so they travel inside the notebook.

In [ ]:

figures = ""
for name, video_path in videos.items():
    encoded = base64.b64encode(video_path.read_bytes()).decode()
    figures += (
        '<figure style="margin:14px auto;max-width:900px;text-align:center">'
        f'<video src="data:video/mp4;base64,{encoded}" style="width:100%"'
        ' autoplay loop muted playsinline controls></video>'
        f'<figcaption><b>Hurricane {name}</b></figcaption></figure>'
    )
HTML(figures)

figures = "" for name, video_path in videos.items(): encoded = base64.b64encode(video_path.read_bytes()).decode() figures += ( '<figure style="margin:14px auto;max-width:900px;text-align:center">' f'<video src="data:video/mp4;base64,{encoded}" style="width:100%"' ' autoplay loop muted playsinline controls></video>' f'<figcaption><b>Hurricane {name}</b></figcaption></figure>' ) HTML(figures)

What you are watching¶

• Ida (24 Aug – 2 Sep 2021) — opens as scattered, whirl-less convection over the central Caribbean, organises into a tropical storm, crosses western Cuba, rapidly intensifies over the warm Gulf into a compact Category 4 with a crisp eye, slams into Louisiana near Port Fourchon, then — inland — loses its eye and unwinds into a big sheared comma that sweeps up over the eastern US and finally smears out into a frontal band off the Northeast, gone.
• Milton (4–11 Oct 2024) — starts as a broad, formless low in the south-western Gulf, coils up during an explosive rapid-intensification burst to a pinhole-eyed Category 5, sweeps east across the Gulf, crosses the Florida peninsula (the actual landfall near Sarasota came after dark, between two daytime frames), and races off into the Atlantic where it comes apart and slides off the eastern edge of GOES-East's view.
• Each day contributes a few daytime frames; the jumps between days are the nights we skipped.

Make it your own: change STORMS to any storm — its days range, the daytime hours to sample, and a bbox wide enough for the track (satellite="18" is GOES-West for the eastern Pacific; domain="F" is the full disc, domain="M1"/"M2" the 1-minute mesoscale sectors); add more hours for a smoother loop.