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¶

1. Download daily GPM IMERG precipitation for 25–31 Aug 2017.
2. Accumulate the daily rainfall into the storm total over Texas.
3. 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. xarray does the cropping and accumulation of the downloaded granules; earthlens provides the unified EarthLens entry point.

In [1]:

import os
import tempfile

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from loguru import logger

from earthlens import EarthLens
from earthlens.base import AuthenticationError

logger.disable("earthlens")
logger.disable("pyramids")

import os import tempfile import matplotlib.pyplot as plt import numpy as np import xarray as xr from loguru import logger from earthlens import EarthLens from earthlens.base import AuthenticationError 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).

In [2]:

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(),
)

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.

In [3]:

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")

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")

C:\python-environments\pixi\envs\earthlens-7782220113782294872\envs\dev\Lib\site-packages\earthaccess\results.py:348: FutureWarning: As of version 1.0, `DataGranule.size` will be accessed as an attribute; e.g. use `DataCollection.size` **not** `DataCollection.size()`
  self["size"] = self.size()
C:\python-environments\pixi\envs\earthlens-7782220113782294872\envs\dev\Lib\site-packages\earthaccess\store.py:838: FutureWarning: As of version 1.0, `DataGranule.size` will be accessed as an attribute; e.g. use `DataCollection.size` **not** `DataCollection.size()`
  total_size = round(sum(granule.size() for granule in granules) / 1024, 2)

7 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.

In [4]:

total = None
if paths:
    for p in sorted(paths):
        da = xr.open_dataset(p)["precipitation"].squeeze()  # dims (lon, lat), mm/day
        da = da.sel(lon=slice(-99, -93), lat=slice(27, 32))  # crop to SE Texas
        total = da.copy() if total is None else total + da.values
    total = total.transpose("lat", "lon")
    print(
        f"storm-total rainfall: peak {float(total.max()):.0f} mm "
        f"({float(total.max())/25.4:.0f} in) over the cropped box"
    )

total = None if paths: for p in sorted(paths): da = xr.open_dataset(p)["precipitation"].squeeze() # dims (lon, lat), mm/day da = da.sel(lon=slice(-99, -93), lat=slice(27, 32)) # crop to SE Texas total = da.copy() if total is None else total + da.values total = total.transpose("lat", "lon") print( f"storm-total rainfall: peak {float(total.max()):.0f} mm " f"({float(total.max())/25.4:.0f} in) over the cropped box" )

C:\python-environments\pixi\envs\earthlens-7782220113782294872\envs\dev\Lib\site-packages\xarray\backends\plugins.py:110: RuntimeWarning: Engine 'cfgrib' loading failed:
cannot load library 'C:\python-environments\pixi\envs\earthlens-7782220113782294872\envs\dev\lib\libeccodes.dll': error 0x7e
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)

storm-total rainfall: peak 861 mm (34 in) over the cropped box

3 · Map the storm total¶

The bullseye of >1 m sits over the Houston–Beaumont corridor, where Harvey's rainbands stalled for days.

In [5]:

if total is not None:
    fig, ax = plt.subplots(figsize=(9, 7))
    im = ax.pcolormesh(total.lon, total.lat, total.values, cmap="turbo", shading="auto")
    ax.plot(-95.37, 29.76, "k*", ms=14)
    ax.annotate("Houston", (-95.37, 29.76), fontsize=10, fontweight="bold")
    ax.set(
        xlabel="longitude",
        ylabel="latitude",
        title="Hurricane Harvey storm-total rainfall, 25–31 Aug 2017 (GPM IMERG)",
    )
    fig.colorbar(im, ax=ax, label="rainfall (mm)")
    plt.show()
else:
    print("no granules — set EARTHDATA_USERNAME / EARTHDATA_PASSWORD and rerun")

if total is not None: fig, ax = plt.subplots(figsize=(9, 7)) im = ax.pcolormesh(total.lon, total.lat, total.values, cmap="turbo", shading="auto") ax.plot(-95.37, 29.76, "k*", ms=14) ax.annotate("Houston", (-95.37, 29.76), fontsize=10, fontweight="bold") ax.set( xlabel="longitude", ylabel="latitude", title="Hurricane Harvey storm-total rainfall, 25–31 Aug 2017 (GPM IMERG)", ) fig.colorbar(im, ax=ax, label="rainfall (mm)") 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 xarray 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.