Search and Stream NISAR Data with asf_search

from pathlib import Path
import sys
project_root = Path.cwd().parent
sys.path.insert(0, str(project_root))

from util.aoi_coverage import plot_aoi_coverage
from util.wkt import rectangular_aoi_select, rectangular_wkt_bounds

aoi_map = rectangular_aoi_select()
print("Click the ■ button to select an area of interest\n") 
display(aoi_map["map"])

from datetime import datetime

start_date = datetime(2025, 11, 22)
end_date = datetime(2028, 3, 5)
area_of_interest = "POLYGON((-90.31002 29.902604, -90.250797 29.902604, -90.250797 29.924476, -90.31002 29.924476, -90.31002 29.902604))" # POINT or POLYGON as WKT (well-known-text)
pattern = r'^(?!.*QA_STATS).*'

product_type = "GCOV" # "RSLC", "RIFG", "RUNW", "ROFF", "GSLC", "GCOV", "GUNW", "SME2", "GOFF"

import os
import asf_search as asf
from pprint import pprint

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

session = asf.ASFSession()

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

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

import re
from util.aoi_coverage import plot_aoi_coverage

pattern = re.compile(r"(NISAR_L[123]_PR_[GRS][COSUIM][OFLNE][CGFVW2](?:_[^_]+){9,11})_([PX]\d{5})")

latest_version_dict = {}

for url in 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)

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

gpolygons = [r.umm['SpatialExtent']['HorizontalSpatialDomain']['Geometry']['GPolygons'][0] for r in response]
display(plot_aoi_coverage(gpolygons, area_of_interest, product_type))

from getpass import getpass

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

%%time
import fsspec
import xarray as xr
import rioxarray
from urllib.parse import urlparse

product_type = urlparse(urls[0]).path.split('_')[2]
# load only the frequency A data
if product_type[0] == "G":
    group_path = f"/science/LSAR/{product_type}/grids/frequencyA"
elif product_type[0] == "R":
    group_path = f"/science/LSAR/{product_type}/swaths/frequencyA"
elif product_type[0] == "S":
    group_path = f"/science/LSAR/{product_type}/grids"
print(f"group: {group_path}\n")

fs = fsspec.filesystem(
    "http",
    headers = {"Authorization": f"Bearer {token}"},
    cache_type = "background",
    block_size = 8 * 1024 * 1024, # 8 MB
)

dt = xr.open_datatree(
        fs.open(urls[0], "rb"),
        engine="h5netcdf",
        decode_timedelta=False,
        phony_dims="access",
        chunks="auto",
        group=group_path, # comment group_path to load the full H5 dataset
)
dt

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

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

import json
import s3fs
import urllib

prefix = f"NISAR_L2_{product_type}_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 rioxarray

product_type = urlparse(urls[0]).path.split('_')[2]
# load only the frequency A data
if product_type[0] == "G":
    group_path = f"/science/LSAR/{product_type}/grids/frequencyA"
elif product_type[0] == "R":
    group_path = f"/science/LSAR/{product_type}/swaths/frequencyA"
elif product_type[0] == "S":
    group_path = f"/science/LSAR/{product_type}/grids/radarData/frequencyA"
print(f"group: {group_path}\n")

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

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

file = fs.open(s3_urls[0], "rb", **kwargs) 

dt = xr.open_datatree(
    file,
    engine="h5netcdf",
    decode_timedelta=False,
    phony_dims="access",
    chunks="auto",
    group=group_path, # comment group_path to load the full H5 dataset
)
dt

%%time
import fsspec
import xarray as xr
import rioxarray
from urllib.parse import urlparse

product_type = urlparse(urls[0]).path.split('_')[2]
# load only the frequency A data
if product_type[0] == "G": # geocoded level-2 data
    group_path = f"/science/LSAR/{product_type}/grids/frequencyA"
elif product_type[0] == "R": # range-doppler level-1 data
    group_path = f"/science/LSAR/{product_type}/swaths/frequencyA"
elif product_type == "SME2":
    group_path = f"/science/LSAR/{product_type}/grids"

fs = fsspec.filesystem(
    "http",
    headers = {"Authorization": f"Bearer {token}"},
    cache_type = "background",
    block_size = 8 * 1024 * 1024, # 8 MB
)

files = [fs.open(url, "rb") for url in urls]

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

from shapely import wkt

subset_datatrees = []
for tree in datatrees:
    if product_type == "GCOV":
        projection = tree.projection.attrs['epsg_code'].item()
        var = tree.HHHH
    elif product_type == "GSLC":
        projection = tree.projection.attrs['epsg_code'].item()
        var = tree.HH
    elif product_type == "SME2":
        projection = tree.projection.attrs['epsg_code'].item()
        var = tree.soilMoisture
    else:
        raise Exception("The spatial subsetting and time series example requires GSLC, GCOV, or SME2 data")
    
    var = var.rio.write_crs(projection)

    # subset to area_of_interest
    geom = wkt.loads(area_of_interest)
    min_lon, min_lat, max_lon, max_lat = geom.bounds
    subset_datatrees.append(var.rio.clip_box(
        minx=min_lon, miny=min_lat,
        maxx=max_lon, maxy=max_lat, 
        crs="EPSG:4326"
    ))

subset_datatrees[0]

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 urls]

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

dataarrays[0].time

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

import numpy as np

if product_type != "SME2":
    ts = 10 * np.log10(ts.where(ts > 0)) # Convert power to dB scale for plotting
    
if product_type == "GSLC":
    ts = ts.real

vmin = ts.isel(time=0).quantile(0.05, skipna=True).values.item()
vmax = ts.isel(time=0).quantile(0.95, skipna=True).values.item()

if product_type == "SME2":
    cmap = "YlGnBu"
else:
    cmap = "gray"   

ts.isel(time=0).plot(vmin=vmin, vmax=vmax, cmap=cmap)

for f in files:
    f.close()