Plot GLM data on a map

[blaylockbk]

Had someone ask about one of my old GLM scripts. Apparently NWS was using some of my old code to make these plots. I had better luck having an LLM create a new script from scratch that gets GLM data from AWS for a time range and plots the flash locations on a map.

import os
from contextlib import suppress
from datetime import datetime, timedelta

import cartopy.crs as ccrs
import fsspec
import matplotlib.pyplot as plt
import xarray as xr


def get_glm_files(
    satellite: str,
    start_time: datetime,
    end_time: datetime,
) -> list[str]:
    """
    List GOES GLM files from AWS Open Data within the given time range.

    Parameters
    ----------
    satellite : {"noaa-goes16", "noaa-goes17", "noaa-goes18", "noaa-goes19"}
    """
    fs = fsspec.filesystem("s3", anon=True)
    base_url = f"s3://{satellite}/GLM-L2-LCFA"

    files = []
    current_time = start_time

    while current_time <= end_time:
        year = current_time.strftime("%Y")
        day_of_year = current_time.strftime("%j")
        hour = current_time.strftime("%H")
        path = f"{base_url}/{year}/{day_of_year}/{hour}/"

        with suppress(FileNotFoundError):
            for file in fs.ls(path):
                filename = os.path.basename(file)
                file_time_str = filename.split("_")[3]
                file_time = datetime.strptime(file_time_str, "s%Y%j%H%M%S%f")
                if start_time <= file_time <= end_time:
                    files.append(f"s3://{file}")

        current_time += timedelta(hours=1)

    return files


def load_glm(files: list[str]) -> list[xr.Dataset]:
    """
    Load GLM NetCDF files from S3 into memory.
    """
    fs = fsspec.filesystem("s3", anon=True)
    datasets = []

    for f in files:
        with fs.open(f) as file:
            ds = xr.open_dataset(file, engine="h5netcdf").load()
            datasets.append(ds)

    return datasets


def plot_glm_events(datasets: list[xr.Dataset]):
    """
    Plot GLM flash event locations on a map using Cartopy.
    """
    fig = plt.figure(figsize=(10, 6))
    ax = plt.axes(projection=ccrs.PlateCarree())
    ax.coastlines()

    for ds in datasets:
        ax.scatter(
            ds.flash_lon.values,
            ds.flash_lat.values,
            s=1,
            color="red",
            transform=ccrs.PlateCarree(),
            label="GLM Flashes",
        )

    ax.set_title("GOES-19 GLM Lightning Flashes")
    plt.legend()
    plt.tight_layout()
    plt.show()


# Example usage
if __name__ == "__main__":
    start = datetime(2025, 7, 1, 22, 0)
    end = datetime(2025, 7, 1, 22, 1)
    satellite = "noaa-goes19"

    print(f"Fetching GLM files from {start} to {end}...")
    files = get_glm_files(satellite, start, end)

    print(f"Found {len(files)} files.")

    if files:
        datasets = load_glm(files)
        print("Data loaded.")
        plot_glm_events(datasets)
    else:
        print("No files found for the specified time range.")