Search and Stream NISAR GCOV Data with asf_search

import os
import asf_search as asf
from datetime import datetime
from getpass import getpass

import warnings
warnings.filterwarnings(
    "ignore",
    message="Parsing dates involving a day of month without a year specified",
)

session = asf.ASFSession()

start_date = datetime(2025, 11, 22)
end_date = datetime(2026, 3, 5)
area_of_interest = "POLYGON((40.9131 12.3904,41.8891 12.3904,41.8891 13.2454,40.9131 13.2454,40.9131 12.3904))" # POINT or POLYGON as WKT (well-known-text)
pattern = r'^(?!.*QA_STATS).*'

opts=asf.ASFSearchOptions(**{
    "maxResults": 250,
    "intersectsWith": area_of_interest,
    "flightDirection": "ASCENDING",
    "start": start_date,
    "end": end_date,
    "processingLevel": [
        "GCOV"
    ],
    "dataset": [
        "NISAR"
    ],
    "productionConfiguration": [
        "PR"
    ],
    'session': session,
})

response = asf.search(opts=opts)
hdf5_urls = response.find_urls(extension='.h5', pattern=pattern, directAccess=False)
print(f"Found {len(hdf5_urls)} GCOV products:")
hdf5_urls

import re

pattern = re.compile(r"(NISAR_L2_PR_GCOV(?:_[^_]+){9})_(X\d{5})")

latest_version_dict = {}

for url in hdf5_urls:
    m = pattern.search(url)
    if not m:
        continue

    product, crid = m.groups()

    if product not in latest_version_dict or crid > latest_version_dict[product][0]:
        latest_version_dict[product] = (crid, url)

hdf5_urls = [i[1] for i in latest_version_dict.values()]
print(f"Retained {len(hdf5_urls)} GCOV products:")
hdf5_urls

from getpass import getpass

token = getpass("Enter your EDL Bearer Token")

%%time
import asf_search as asf
import fsspec
import rioxarray
import xarray as xr

fs = fsspec.filesystem(
    "http",
    headers = {"Authorization": f"Bearer {token}"},
    block_size = 16 * 512 * 512,
)

with fs.open(hdf5_urls[0], "rb") as f:
    dt = xr.open_datatree(
        f,
        engine="h5netcdf",
        decode_timedelta=False,
        phony_dims="access",
        chunks="auto",
    )

    ### Perform any calculations and save any computed results for future access here, inside the `with` block ###
    frequencyA = dt["/science/LSAR/GCOV/grids/frequencyA"] # access the frequency A data
    projection = frequencyA.projection.attrs['epsg_code'].item() # access the GCOV product's projection
    hhhh = frequencyA.HHHH # access frequency A's HHHH band
    hhhh = hhhh.rio.write_crs(projection) # write the project to the HHHH data for easy lat/lon subsetting

    # subset the data
    subset_hhhh = hhhh.rio.clip_box(
        minx=40.8463, miny=13.2553,
        maxx=40.8574, maxy=13.2684, 
        crs="EPSG:4326"
    )

    # save the mean of HHHH subset for use outside of the `with` block
    subset_hhhh_mean = subset_hhhh.mean().to_numpy().item()
    print(f'subset_hhhh_mean: {subset_hhhh_mean}\n')

from urllib.parse import urlparse

bucket = "sds-n-cumulus-prod-nisar-products"

s3_urls = [f"s3://{bucket}/{'/'.join(urlparse(url).path.split('/')[2:])}" for url in hdf5_urls]
s3_urls

import json
import s3fs
import urllib

prefix = "NISAR_L2_GCOV_BETA_V1"

event = {
    "CredentialsEndpoint": "https://nisar.asf.earthdatacloud.nasa.gov/s3credentials",
    "BearerToken": token,
    "Bucket": bucket,
    "Prefix": prefix,
}

# Get temporary download credentials
tea_url = event["CredentialsEndpoint"]
bearer_token = event["BearerToken"]
req = urllib.request.Request(
    url=tea_url,
    headers={"Authorization": f"Bearer {bearer_token}"}
)
with urllib.request.urlopen(req) as f:
    creds = json.loads(f.read().decode())

%%time
import xarray as xr
import s3fs
import rioxarray

s3_url = s3_urls[0]

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

kwargs = {
    "cache_type": "background",
    "block_size": 16 * 512 * 512,  # 16 MB
}

with fs.open(s3_url, "rb", **kwargs) as f:
    dt = xr.open_datatree(
        f, 
        engine="h5netcdf", 
        decode_timedelta=False, 
        chunks="auto",
        phony_dims="access"
    )

    ### Perform any calculations and save any computed results for future access here, inside the `with` block ###
    frequencyA = dt["/science/LSAR/GCOV/grids/frequencyA"] # access the frequency A data
    projection = frequencyA.projection.attrs['epsg_code'].item() # access the GCOV product's projection
    hhhh = frequencyA.HHHH # access frequency A's HHHH band
    hhhh = hhhh.rio.write_crs(projection) # write the project to the HHHH data for easy lat/lon subsetting

    # subset the data
    subset_hhhh = hhhh.rio.clip_box(
        minx=40.8463, miny=13.2553,
        maxx=40.8574, maxy=13.2684, 
        crs="EPSG:4326"
    )

    # save the mean of HHHH subset for use outside of the `with` block
    subset_hhhh_mean = subset_hhhh.mean().to_numpy().item()
    print(f'subset_hhhh_mean: {subset_hhhh_mean}\n')

dt

dt.compute()

subset_hhhh_mean

%%time
import xarray as xr
import rioxarray

# Explore the DataTree rendering above in Step 4 for a complete list of available groups 
group_path = "/science/LSAR/GCOV/grids/frequencyA" # change this to any GCOV HDF5 group you wish

kwargs = {
    "cache_type": "background",
    "block_size": 16 * 1024 * 1024,  # 16 MB
}

files = [fs.open(url, "rb", **kwargs) for url in s3_urls]

datatrees = [
    xr.open_datatree(
        f,
        engine="h5netcdf",
        decode_timedelta=False,
        phony_dims="access",
        chunks="auto",
        group=group_path,
    )
    for f in files
]

for tree in datatrees:
    projection = tree.projection.attrs['epsg_code'].item()
    hhhh = tree.HHHH
    hhhh = hhhh.rio.write_crs(projection)
    
    subset_hhhh = hhhh.rio.clip_box(
        minx=40.8463, miny=13.2553,
        maxx=40.8574, maxy=13.2684, 
        crs="EPSG:4326"
    )
    subset_hhhh_mean = subset_hhhh.mean().to_numpy().item()
    print(subset_hhhh_mean)

import re
from urllib.parse import urlparse
from datetime import datetime
from pathlib import PurePosixPath

NISAR_TS_RE = re.compile(r"_(\d{8}T\d{6})_")

def nisar_start_time_from_url(s3_url: str) -> datetime:
    path = urlparse(s3_url).path
    name = PurePosixPath(path).name
    
    m = NISAR_TS_RE.search(name)
    if not m:
        raise ValueError(f"No NISAR timestamp found in: {s3_url}")
    
    return datetime.strptime(m.group(1), "%Y%m%dT%H%M%S")

dts = [nisar_start_time_from_url(url) for url in s3_urls]
dts

dataarrays = [
    tree.ds.assign_coords(time=dt).expand_dims(time=1)
    for dt, tree in zip(dts, datatrees)
]

dataarrays[0].time

for da in dataarrays:
    print(da.dims)

ts = xr.concat(dataarrays, dim="time")
ts

for f in files:
    f.close()

ts.compute()