Lightning over the Americas, hour by hour: 373 million GOES-19 flashes

[jvdd]

One year of lightning over the Americas, observed by the GOES-19 Geostationary Lightning Mapper (GLM): 373M flashes from September 2025 to August 2026, grouped by local solar time.
-> Brush the hour histogram and the map shows where lightning strikes at those hours. Zoom the map and it re-bins to the new view.

• Data: NOAA GOES-19 GLM Level 2 flashes, public on AWS: https://registry.opendata.aws/noaa-goes/ (bucket noaa-goes19, prefix GLM-L2-LCFA/)
• Install: pip install "flexviz>=0.1.0b5" h5py (color_norm="log" is new in 0.1.0b5)
• Size: one year is about 770 GB of NetCDF files to download. The scripts read them in memory and keep only the flashes: 4.8 GB of Parquet, 3.6 GB after cleaning.

Steps

1. python etl_glm.py 2025-09-01 2026-08-31 downloads the files and writes one Parquet file of flashes per day. It runs every 7th day first, so a stopped run still covers all seasons.
2. python clean.py removes GLM artifacts and adds the local solar time (UTC + longitude / 15).
3. python dashboard.py opens the live dashboard. It keeps the 373 million rows in memory (about 4.5 GB), so a brush answers in about a second.

Cleaning

clean.py removes 5.8 % of the flashes:

• Stuck pixels. A few GLM pixels, mostly over open ocean, report up to tens of thousands of one-pixel "flashes" a day at one spot, always with the same flash sizes. A real storm's flashes vary in size, so the script drops the (0.1° cell, flash area) pairs of those pixels on every day.
• False detections in August 2026. On 17, 19, 23, 24 and 25 August, GLM reported hundreds of thousands of weak false flashes over the south-east Pacific. Those 5 days are left out.

NOAA's archive has almost no data for 15-17 July 2026.

GLM misses more flashes in daylight than at night, because sunlit cloud tops hide weak flashes. So the afternoon peak in this data is, if anything, smaller than the real one.

The video

The video is a headless-browser capture of this dashboard in its dark theme. The brush moves in 15-minute steps over one-hour windows, with a fixed log color range (1 to 2,000 flashes per 0.1° cell), so the hours compare with each other. The title and the clock are overlays.

Code

etl_glm.py: download and extract the flashes

"""Extract GOES-19 GLM L2 flashes to one Parquet file per UTC day.

NetCDF files are read from memory and never written to disk. Only flashes
with quality flag 0 are kept. Re-running skips days that already exist.
Days run every-7th-day first, so a stopped run still covers all seasons.

Usage: python etl_glm.py 2025-09-01 2026-08-31
"""

import datetime as dt
import io
import re
import sys
import time
import urllib.request
from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor
from pathlib import Path

import h5py
import numpy as np
import polars as pl

BUCKET = "https://noaa-goes19.s3.amazonaws.com"
OUT = Path(__file__).parent / "data" / "flashes"
LOG = Path(__file__).parent / "data" / "etl.log"
HOUR_WORKERS = 8
DOWNLOAD_THREADS = 16


def _get(url: str, tries: int = 4) -> bytes | None:
    for attempt in range(tries):
        try:
            with urllib.request.urlopen(url, timeout=60) as r:
                return r.read()
        except Exception:
            time.sleep(2**attempt)
    return None


def _list_keys(prefix: str) -> list[str]:
    body = _get(f"{BUCKET}/?list-type=2&prefix={prefix}&max-keys=1000")
    if body is None:
        raise RuntimeError(f"listing failed: {prefix}")
    return re.findall(r"<Key>([^<]+)</Key>", body.decode())


def _unpack(ds) -> np.ndarray:
    # GLM stores these as scaled, unsigned 16-bit integers.
    raw = ds[:].view(np.uint16).astype(np.float64)
    return raw * float(ds.attrs["scale_factor"][0]) + float(ds.attrs["add_offset"][0])


def _parse(blob: bytes) -> dict[str, np.ndarray] | None:
    with h5py.File(io.BytesIO(blob), "r") as f:
        if f["flash_lat"].shape[0] == 0:
            return None
        keep = f["flash_quality_flag"][:] == 0
        offset_ds = f["flash_time_offset_of_first_event"]
        units = offset_ds.attrs["units"].decode()
        base = np.datetime64(
            units.removeprefix("seconds since ").replace(" ", "T"), "ms"
        )
        offset_ms = np.round(_unpack(offset_ds) * 1000).astype("timedelta64[ms]")
        return {
            "time": (base + offset_ms)[keep],
            "lat": f["flash_lat"][:][keep],
            "lon": f["flash_lon"][:][keep],
            "energy_fj": (_unpack(f["flash_energy"]) * 1e15).astype(np.float32)[keep],
            "area_km2": (_unpack(f["flash_area"]) / 1e6).astype(np.float32)[keep],
        }


def _hour(prefix: str) -> tuple[list[dict], int, int, int]:
    keys = _list_keys(prefix)
    with ThreadPoolExecutor(DOWNLOAD_THREADS) as pool:
        blobs = list(pool.map(lambda k: _get(f"{BUCKET}/{k}"), keys))
    parts = [p for b in blobs if b is not None and (p := _parse(b)) is not None]
    missing = sum(b is None for b in blobs)
    return parts, len(keys), missing, sum(len(b) for b in blobs if b)


def _day(day: dt.date, pool: ProcessPoolExecutor) -> str:
    out = OUT / f"{day.isoformat()}.parquet"
    if out.exists():
        return f"{day} skip"
    t0 = time.time()
    doy = day.timetuple().tm_yday
    prefixes = [f"GLM-L2-LCFA/{day.year}/{doy:03d}/{h:02d}/" for h in range(24)]
    parts, n_files, n_missing, n_bytes = [], 0, 0, 0
    for p, n, m, b in pool.map(_hour, prefixes):
        parts += p
        n_files += n
        n_missing += m
        n_bytes += b
    if not parts:
        return f"{day} no flashes in {n_files} files"
    df = pl.DataFrame(
        {k: np.concatenate([p[k] for p in parts]) for k in parts[0]}
    ).sort("time")
    tmp = out.with_suffix(".tmp")
    df.write_parquet(tmp, compression="zstd")
    tmp.rename(out)
    return (
        f"{day} flashes={df.height} files={n_files} missing={n_missing} "
        f"MB={n_bytes / 1e6:.0f} s={time.time() - t0:.0f}"
    )


def main() -> None:
    start, end = (dt.date.fromisoformat(a) for a in sys.argv[1:3])
    OUT.mkdir(parents=True, exist_ok=True)
    days = [start + dt.timedelta(d) for d in range((end - start).days + 1)]
    # Every 7th day first: a partial run is still a sample of every season.
    days.sort(key=lambda d: ((d - start).days % 7, d))
    with ProcessPoolExecutor(HOUR_WORKERS) as pool, LOG.open("a") as log:
        for day in days:
            line = _day(day, pool)
            print(line, flush=True)
            log.write(line + "\n")
            log.flush()


if __name__ == "__main__":
    main()

clean.py: remove artifacts, add local solar time

"""Remove stuck-pixel artifacts and add local solar time.

GLM has a few stuck pixels, mostly over open ocean, that report hundreds to
tens of thousands of one-pixel "flashes" a day at one spot, always with the
same one or two flash areas. A real storm's flashes vary in size. So a
cell-day with at least 200 flashes but fewer distinct areas than 2 % of the
flash count is a stuck pixel (artifacts stay below 1.2 %, real storms above
3.3 %). Its (0.1 degree cell, flash area) pairs are then dropped on every
day: that also catches the quiet days on which the pixel fires fewer than
200 times, and it keeps the real flashes in the same cell, which have other
areas.

In August 2026 GLM also reported hundreds of thousands of weak false flashes
over the south-east Pacific (136-72 W, 30-5 S), mostly between 8 and 12 local
time. On every other day that box gets at most 20,550 such morning flashes, so
a day with more than 30,000 is left out (5 days: 2026-08-17, 19, 23, 24, 25).

Usage: python clean.py   (rewrites data/clean/, prints what it removed)
"""

from pathlib import Path

import polars as pl

DATA = Path(__file__).parent / "data"
RAW, CLEAN = DATA / "flashes", DATA / "clean"
CELL = ["lon10", "lat10"]
KEY = [*CELL, "area_km2"]
FALSE_MORNING = (
    pl.col("lon").is_between(-136, -72)
    & pl.col("lat").is_between(-30, -5)
    & pl.col("local_hour").is_between(8, 12)
)


def with_cell(df: pl.DataFrame) -> pl.DataFrame:
    return df.with_columns(
        lon10=(pl.col("lon") * 10).floor(), lat10=(pl.col("lat") * 10).floor()
    )


def stuck_pairs(df: pl.DataFrame) -> pl.DataFrame:
    """(cell, area) pairs of the stuck pixels in one day."""
    cells = (
        df.group_by(CELL)
        .agg(n=pl.len(), areas=pl.col("area_km2").n_unique())
        .filter((pl.col("n") >= 200) & (pl.col("areas") < 0.02 * pl.col("n")))
    )
    return df.join(cells.select(CELL), on=CELL).select(KEY).unique()


def main() -> None:
    days = sorted(RAW.glob("*.parquet"))
    stuck = pl.concat(
        [stuck_pairs(with_cell(pl.read_parquet(src))) for src in days]
    ).unique()
    print(
        f"stuck (cell, area) pairs: {stuck.height} in "
        f"{stuck.select(CELL).unique().height} cells"
    )

    CLEAN.mkdir(exist_ok=True)
    total = removed = 0
    t = pl.col("time")
    utc_hour = t.dt.hour() + t.dt.minute() / 60 + t.dt.second() / 3600
    skipped = []
    for src in days:
        df = with_cell(pl.read_parquet(src))
        out = df.join(stuck, on=KEY, how="anti").select(
            "lat",
            "lon",
            local_hour=((utc_hour + pl.col("lon") / 15) % 24).cast(pl.Float32),
        )
        total += df.height
        if out.filter(FALSE_MORNING).height > 30_000:
            skipped.append(src.stem)
            (CLEAN / src.name).unlink(missing_ok=True)
            removed += df.height
            continue
        out.write_parquet(CLEAN / src.name, compression="zstd")
        removed += df.height - out.height
    print(f"days left out: {skipped}")
    print(
        f"flashes {total}, removed {removed} ({100 * removed / total:.2f} %), "
        f"kept {total - removed}"
    )


if __name__ == "__main__":
    main()

dashboard.py: the live dashboard

"""Live FlexViz dashboard: one year of GOES-19 lightning flashes.

Brush the hour histogram and the map shows where lightning strikes at those
hours. Zoom the map and it re-bins to the new view.

Run: python dashboard.py   (data/clean/ comes from etl_glm.py and clean.py)
"""

import polars as pl
from flexviz import Dashboard, GridItem, LayoutSpec

# 373 million flashes: lat, lon and local solar hour. collect() keeps them in
# RAM (about 4.5 GB), so a brush answers in about a second. Pass the LazyFrame
# instead to stream them from disk.
flashes = pl.scan_parquet("data/clean/*.parquet").collect()

dash = Dashboard(flashes, cache=True)

lightning_map = dash.add_figure(title="Lightning flashes, Sep 2025 – Aug 2026")
lightning_map.add_histogram2d(
    x="lon",
    y="lat",
    x_bins=200,
    y_bins=100,
    # x_bins=1320,  # about 0.1 degree per bin
    # y_bins=1150,
    color_scale="Electric",
    color_norm="log",
)
lightning_map.update_layout(
    plot_bgcolor="black",
    xaxis={"visible": False},
    yaxis={"visible": False, "scaleanchor": "x"},
)

hours = dash.add_figure(title="Local solar time")
hours.add_histogram(x="local_hour", bins=96).xlabel("hour of day")

map_uid, hours_uid = (f.uid for f in dash.to_spec().figures)
dash.show(
    port=8050,
    source_name="lightning",
    # A live-brush cube for the 1320x1150 map would hold ~170M cells: too big
    # for the browser, so brushes use the normal request on mouse-up.
    live_brush="off",
    layout=LayoutSpec(
        grid_items=[
            GridItem(fig_uid=map_uid, x=0, y=0, w=12, h=8),
            GridItem(fig_uid=hours_uid, x=0, y=8, w=12, h=2),
        ]
    ),
)

[jvdd]

Reddit post: https://www.reddit.com/r/dataisbeautiful/comments/1wsd5lb/oc_lightning_throughout_the_day_over_the_americas/