[Data]: IMERG Half Hourly in Cloud-native format

[mccabete]

Dataset Name

GPM IMERG Final Precipitation L3 Half Hourly 0.1 degree x 0.1 degree V07

Dataset Description

This is IMERG's multi-sensor precipitation product at the half-hourly resolution. It is stored in HDF5 file format.

Requestor Name/Affiliation

Tempest McCabe, NASA GSFC, University of Maryland

Platform/Method/Sensor

GPM

URL or DOI to Dataset Description

10.5067/GPM/IMERG/3B-HH/07

URL to Download or Access the Data

https://disc.gsfc.nasa.gov/datasets/GPM_3IMERGHH_07/summary?keywords=%22IMERG%20final%22

Data License

Unknown, but NASA product so freely available and usable.

Intended Science Use Case

I will be extracting a precipitation time series for ~4000 fires in CONUS. This involves a spatial and temporal subset, and then an aggregation to 12-hourly. This is different from, but related to, the conversion this team did for #1281. My workflow was previously built around the already-cloud-optimized raw daily product (10.5067/GPM/IMERGDF/DAY/07) but we realized that the daily aggregation happened in UTC, so if we wanted out time series to align with local-time satellite overpasses, we needed to go to half-hourly.

Unlike the daily product, the half-hourly product is in HDF5 and I think isn't cloud optimized, so when I try to do my spatial and temporal subsetting with xarray, all the files are read into memory. For example, I tried to do a point-extraction for two week-long fire events and it took ~ 2GB. (7 days * 48 * 2.8MB per file) * 2 fires = 1.9GB. I am not sure how to build a workflow that can handle that for 4000 fires.

It the MAAP team is the correct team to contact about this? Would earthaccess be more appropriate? A comment on a different earthaccess discussion suggests to me that some folks may be already working on some element of this.

Format of the Data

HDF5

Approximate Size of the Data

2.8MB per 30 min file, 1998-01-01 to 2025-09-30 in total.

Date Needed By

05/25/26

Additional Information

I would be working in a spatial range of CONUS 2012-2025.

[wildintellect]

@mccabete I think this may be available via Virtual Zarr https://source.coop/bkr/imerg
@hrodmn @abarciauskas-bgse can you confirm this is the correct data, and provide a link to an access example?

[hrodmn]

I don't know if that dataset on source coop is reliable. I tried accessing the precipitation values and it's coming back as all zeros which seems... wrong 🤔.

Here is how I opened it on the MAAP Hub, I think it's the right dataset.

import xarray as xr
import fsspec

zarr_url = "s3://bkr/imerg/imerg_final.zarr/"

mapper = fsspec.get_mapper(
    zarr_url,
    client_kwargs={"endpoint_url": "https://data.source.coop"},
    anon=True,
)

ds = xr.open_zarr(mapper, consolidated=True)

I also tried using earthaccess + fsspec + xarray to open a bunch of granules and can confirm that it is slow:

import earthaccess
import fsspec
import s3fs
import xarray as xr

# authenticate with earthaccess
auth = earthaccess.login()

# run a search with earthaccess (this just gets a month's worth of granules)
results = earthaccess.search_data(
    short_name="GPM_3IMERGHH",
    temporal=("2025-01-01T00:00:00Z", "2025-01-31T00:00:00Z"),
    # count=1,
    cloud_hosted=True,
)
print(len(results))

s3_links = [g.data_links(access="direct")[0] for g in results]
print(s3_links[0])

# get S3 credentials for the GESDISC bucket
creds = auth.get_s3_credentials(endpoint=results[0].get_s3_credentials_endpoint())

fs_s3 = s3fs.S3FileSystem(
    anon=False,
    key=creds["accessKeyId"],
    secret=creds["secretAccessKey"],
    token=creds["sessionToken"],
)

readers = [
    fs_s3.open(
        s3_link,
        mode="rb",
        cache_type="blockcache",
        block_size=8 * 1024 * 1024,
    )
    for s3_link in s3_links
]

# this has not finished running yet :/
ds = xr.open_mfdataset(
    readers,
    group="Grid",
    engine="h5netcdf",
    decode_timedelta=False,
    phony_dims="sort",
    parallel=True,
)

Virtualization or conversion to zarr (since the dataset is not being updated) would be really helpful here - since the source coop zarr copy does not seem to be correct.

[wildintellect]

@mccabete one question I have - do you plan to create output for each file individually, or were you planning to build a sparse datacube of the time series for all of CONUS? Just trying to think about the data access patterns, and how the work might be parallelized.

[abarciauskas-bgse]

Ryan Abernathey virtualized GPM IMERG half-hourly but it almost certainly needs to be updated (ref: https://github.com/earth-mover/icechunk-nasa/blob/main/design-docs/icechunk-stores.md). This store was generated as a prototype and (1) has no infrastructure to keep it up-to-date and (2) is using an old version of icechunk, I think even pre-1.0 release, and we should generate a new virtual icechunk store using icechunk 2.0.

Given the use case, this is probably a good use case for virtual dataset --> secondary dataset in native zarr format rechunked for timeseries optimization.

[mccabete]

@wildintellect Don't know if I know enough to answer, but probably sparse. I would be extracting a short (days-months) time series for each fire, tagging it with the fire-id, and writing it out to a geoparquet. The fires are usually smaller than the resolution of the product, so it's ~2/3 point extractions, and 1/3 two grid cells that I take the max of.

So an individual file for say, May 6th 2017 may be called 400 times if there were 400 fires burning that day, but I wouldn't want all of CONUS for all for 2012-2025. Is that helpful?

[wildintellect]

@mccabete yes that helps, sounds to me like you want 1 geoparquet per fire (or at least 1 row, though 1 row per time is more likely). This is good, because it means you can run 400 jobs in parallel. We may end up reading the same chunk a few times, but that's not a big deal.

[mccabete]

Here is some example code for a single fire accessed through earthaccess:

minx  = -101.639160
miny  =  35.738127
maxx  = -101.566161
maxy  =  35.792090
#Search for the granule by DOI
results = earthaccess.search_data(
    doi='10.5067/GPM/IMERG/3B-HH/07', # Half hourly
    temporal=("2019-03-13", "2019-03-14"),
    bounding_box = (minx, miny, maxx, maxy)
)

fn = earthaccess.open(results)

ds = xr.open_mfdataset(fn, engine='h5netcdf',
                       group='Grid',
                       backend_kwargs={'phony_dims': 'sort'},
                       combine='nested',
                       concat_dim='time',
                       preprocess = preprocess_with_bounds)

### extract the times to local 1:30 pm and 1:30 am
max_lon = ds['lon'].values.max()
min_lon = ds['lon'].values.min()

if ((max_lon - min_lon) >= 0.3):
    print(f"ERROR: attempting to find local solar time across a distance between {max_lon} degrees lon and {min_lon} degrees lon.")
else:
    longitude = np.mean(ds['lon'].values)  # Taking mean so that multiple close values of longitude work the same as a point
    desired_lst_start_hour = 1.5  # 1:30 AM
    lst_to_utc_offset_hours = longitude / 15.0  
    utc_anchor_hour = (desired_lst_start_hour - lst_to_utc_offset_hours) % 12
    offset_timedelta = pd.Timedelta(hours=utc_anchor_hour)
    
    # Resample the data natively in UTC
    resampled_da = ds.resample(time="12h", offset=offset_timedelta).sum()

    if (resampled_da["lon"].values.size > 1):
        spatial_dims = ["lon", "lat"]
        agg = getattr(resampled_da['precipitation'], "max")(dim=spatial_dims) 

df = agg.to_dataframe.reset_index()