Showcase — Hurricane Harvey's record rainfall (NASA GPM, Earthdata)¶
When Hurricane Harvey stalled over southeast Texas in late August 2017 it produced the largest rainfall total of any tropical cyclone in US history — over 1.5 m (60 in) of rain, triggering catastrophic flooding in Houston. NASA's GPM IMERG merges a constellation of microwave and infrared satellites into a global half-hourly/daily precipitation grid. This notebook uses the earthlens Earthdata backend to pull the daily IMERG rainfall over Texas across the event and accumulate the storm total.
Uses NASA Earthdata; set
EARTHDATA_USERNAME/EARTHDATA_PASSWORD(a free Earthdata Login); the download skips gracefully without them. IMERG granules are global, so we crop to Texas after download.
What this notebook does¶
- Download daily GPM IMERG precipitation for 25–31 Aug 2017.
- Accumulate the daily rainfall into the storm total over Texas.
- Map the total and report the peak — among the highest tropical-cyclone rainfall totals ever recorded.
Setup¶
First the imports and a bit of noise-suppression. pyramids does the cropping and
accumulation of the downloaded granules; earthlens provides the unified
EarthLens entry point.
import tempfile
import matplotlib.pyplot as plt
import numpy as np
from cleopatra.styling.colorbar import ColorBar
from cleopatra.styling.params import DataStyle
from loguru import logger
from pyramids.dataset import Dataset, GeoReference
from pyramids.netcdf import NetCDF
from earthlens.base import AuthenticationError
from earthlens.core import EarthLens
logger.disable("earthlens")
logger.disable("pyramids")
1 · Download the daily IMERG granules¶
We build the request first — source, dataset, variable, the event date window,
and a bounding box over Texas. This is the daily IMERG Final Run v07, which
serves one global granule per day (we crop to Texas after download).
el = EarthLens(
data_source="earthdata",
dataset="GPM_3IMERGDF_07",
variables=["precipitation"],
start="2017-08-25",
end="2017-08-31",
aoi=[-99.0, 27.0, -93.0, 32.0],
path=tempfile.mkdtemp(),
)
Authentication and the download are kept as separate steps so each is easy to
read and re-run: authenticate() resolves the Earthdata login, then download()
fetches the daily granules. If credentials are missing the notebook skips the
download gracefully instead of erroring.
paths = None
try:
el.authenticate()
except AuthenticationError as exc: # noqa: BLE001
print("Earthdata download skipped (needs credentials):", exc)
else:
paths = el.download(progress_bar=False)
print(f"{len(paths)} daily IMERG granules")
2 · Accumulate the storm-total rainfall¶
Each granule carries precipitation in mm/day on a global 0.1° grid. We crop each day to a box over southeast Texas and sum the days into the storm total.
BOX = [-99.0, 27.0, -93.0, 32.0] # SE Texas: min_lon, min_lat, max_lon, max_lat
total = None
geo_ref = None
if paths:
for granule in sorted(paths):
nc = NetCDF.read_file(granule, read_only=True)
# The daily V07 granules are flat NetCDF-4: `precipitation` sits at the
# root, and `nc.group_names` is empty. (The *monthly* IMERG product is
# HDF5 with everything under a `Grid` group -- a different file layout,
# not a different variable name.)
field = nc.get_variable("precipitation")
clipped = field.crop(bbox=BOX, epsg=4326)
values = clipped.read_array(masked=True).astype("float64")
if values.ndim == 3:
values = values[0]
total = values if total is None else total + values
nc.close()
# Every granule shares the same crop box/grid, so the geo-reference only
# needs computing once -- from the last granule's clipped field, after
# the loop rather than on every iteration.
geo_ref = GeoReference(geo=clipped.geotransform, epsg=clipped.epsg)
print(
f"storm-total rainfall: peak {float(total.max()):.0f} mm "
f"({float(total.max()) / 25.4:.0f} in) over the cropped box"
)
3 · Map the storm total¶
The heaviest accumulation sits over the Houston–Beaumont corridor, where Harvey's rainbands stalled for days.
Note the scale. IMERG's 0.1° cells peak in the hundreds of mm over this box (re-verify against float(total.max()), printed above, on the next live run), while the figure usually quoted for Harvey — over 1.5 m — is a point total from rain gauges. A ~11 km grid cell averages that extreme away, so the two numbers are measuring different things rather than disagreeing.
if total is not None:
storm_total = Dataset.from_array(
total.filled(np.nan), no_data_value=np.nan, geo_ref=geo_ref
)
glyph = storm_total.plot(
data_style=DataStyle(style="total_precipitation"),
colorbar=ColorBar(label="rainfall (mm)"),
title="Hurricane Harvey storm-total rainfall, 25–31 Aug 2017 (GPM IMERG)",
)
glyph.ax.plot(-95.37, 29.76, "k*", ms=14)
glyph.ax.annotate("Houston", (-95.37, 29.76), fontsize=10, fontweight="bold")
glyph.ax.set_xlabel("longitude")
glyph.ax.set_ylabel("latitude")
plt.show()
else:
print("no granules — set EARTHDATA_USERNAME / EARTHDATA_PASSWORD and rerun")
Recap¶
The earthlens Earthdata backend fetches NASA's GPM IMERG granules; a few lines of pyramids crop and accumulate them into the storm-total rainfall map that defined Harvey as the wettest US tropical cyclone on record.
Try it yourself¶
- Switch to the half-hourly product (
GPM_3IMERGHHL_07) for the rainfall-rate time series at a point. - Re-point at another storm (Florence 2018, Ida 2021) or a monsoon.
- See the Earthdata backend reference.