Opening NISAR GCOV HDF5 Files in Python

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

#1: Search for and download NISAR GCOV HDF5 products from ASF

#authenticate with Earthdata Login
#create an ASF session object (stores your login + settings for requests)
session = asf.ASFSession()
#prompt for username/password
session.auth_with_creds(input('EDL Username'), getpass('EDL Password'))

#2: Define the search window and area of interest 

#insert start and end dates. Note that dates and times are inclusive. 
start_date = datetime(2025, 12, 4)
end_date = datetime(2025, 12, 5)
# POINT or POLYGON as WKT (well-known-text) using lat/lon pairs 
area_of_interest = "POLYGON((40.9131 12.3904,41.8891 12.3904,41.8891 13.2454,40.9131 13.2454,40.9131 12.3904))" 

#3: Build search options 

#these filters narrow results to NISAR GCOV products over the AOI and time window
opts = asf.ASFSearchOptions(
    **{
        "maxResults": 250,                    # cap the number of returned results
        "intersectsWith": area_of_interest,   # spatial filter (WKT geometry)
        "start": start_date,                  # temporal filter (start)
        "end": end_date,                      # temporal filter (end)
        "processingLevel": ["GCOV"],          # product type / processing level
        "dataset": ["NISAR"],                 # mission/dataset name
        "productionConfiguration": ["PR"],    # production config (e.g., provisional)
        "session": session,                   # use authenticated session from above
    }
)

#4: Run the search for HDF5 files
#execute search
response = asf.search(opts=opts)

#exclude QA_STATE files 
pattern = r'^(?!.*QA_STATS).*'

#extract files that end with .h5 
hdf5_files = response.find_urls(extension='.h5', pattern=pattern, directAccess=False)
hdf5_files

import re

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

latest_version_dict = {}

for url in hdf5_files:
    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_files = [i[1] for i in latest_version_dict.values()]
print(f"Retained {len(hdf5_files)} GCOV products:")
hdf5_files

# Create a folder and download the files
data_dir = Path.home() / "GCOV_data"
data_dir.mkdir(exist_ok=True)
print(f'data_dir: {data_dir}')

#download files
asf.download_urls(hdf5_files, data_dir, session=session, processes=4)

import h5py

#path to the NISAR GCOV HDF5 file
h5_path = list(data_dir.glob("NISAR_L2_*_GCOV*.h5"))[0]

#open the GCOV HDF5 file in read mode ("r")
#HDF5 files behave like Python dictionaries, where keys correspond to groups
f = h5py.File(h5_path, "r")

#inspect the top-level group names in the file
list(f)

list(f["science"])

#list the contents of the /science/LSAR group
#this shows the available product types and supporting subgroups
list(f["science/LSAR"])

#list the contents of the GCOV group
#this reveals how GCOV data are organized internally
list(f["science/LSAR/GCOV"])

#list the contents of the grids group
#this shows how GCOV raster data are organized by frequency
list(f["science/LSAR/GCOV/grids"])

#list the datasets available within the frequencyA group
#these correspond to specific polarization combinations and supporting layers
list(f["science/LSAR/GCOV/grids/frequencyA"])

#access a covariance (measurement) dataset by its full HDF5 path
#HVHV is a real-valued covariance term for the HV–HV polarization combination
hvhv = f["science/LSAR/GCOV/grids/frequencyA/HVHV"]

#inspect basic properties of the dataset without loading the full array
hvhv.shape, hvhv.dtype

#list all attributes attached to the HVHV dataset
list(hvhv.attrs)

#select key dataset attributes
#before reading the dataset into memory, we inspect a small number of attributes that directly affect how the data should be handled in analysis and visualization

print("Fill value:", hvhv.attrs.get("_FillValue"))
print("Minimum valid value:", hvhv.attrs.get("valid_min"))

#read the full HVHV dataset into memory as a NumPy array
hvhv_data = hvhv[()]

#confirm the array shape
hvhv_data.shape

#read the x and y coordinate vectors for the frequencyA grid
x = f["science/LSAR/GCOV/grids/frequencyA/xCoordinates"][:]
y = f["science/LSAR/GCOV/grids/frequencyA/yCoordinates"][:]

#inspect coordinate ranges and lengths
print(f"x range: {x.min():.1f} to {x.max():.1f}")
print(f"y range: {y.min():.1f} to {y.max():.1f}")
print(f"number of x coordinates (columns): {len(x)}")
print(f"number of y coordinates (rows): {len(y)}")

#read the grid spacing (pixel size) for the frequencyA grid
dx = f["science/LSAR/GCOV/grids/frequencyA/xCoordinateSpacing"][()]
dy = f["science/LSAR/GCOV/grids/frequencyA/yCoordinateSpacing"][()]

#print
print(f"x spacing (dx): {float(dx)}")
print(f"y spacing (dy): {float(dy)}")

#read the projection information for the frequencyA grid
projection = f["science/LSAR/GCOV/grids/frequencyA/projection"][()]

projection

import numpy as np
import rasterio
from rasterio.transform import from_origin

#convert the projection value into a CRS string for rasterio 
crs = f"EPSG:{int(projection)}"

#build pixel footprint
#x.min() gives the left edge of the raster
#y.max() gives the top edge of the raster
#dx is the pixel width
#abs(dy) is the pixel height (use absolute value to ensure positive values)
transform = from_origin(
    x.min(),          #x-coordinate of the upper-left corner
    y.max(),          #y-coordinate of the upper-left corner
    float(dx),        #pixel width (meters)
    float(abs(dy))    #pixel height (meters)
)

#name of the output GeoTIFF file
out_tif = data_dir / "GCOV_frequencyA_HVHV.tif"

#open a new GeoTIFF file 
with rasterio.open(
    out_tif,
    "w",                      #write mode
    driver="GTiff",           #output file format
    height=hvhv_data.shape[0],#number of rows
    width=hvhv_data.shape[1], #number of columns
    count=1,                  #single-band raster
    dtype=hvhv_data.dtype,    #data type matches the NumPy array
    crs=crs,                  #coordinate reference system
    transform=transform,      #georeferencing information
    nodata=np.nan,            #missing values stored as NaNs
) as dst:
    #write the HVHV data array to band 1 of the GeoTIFF
    dst.write(hvhv_data, 1)

#print the output filename for confirmation of file creation
out_tif

#quality check
with rasterio.open(out_tif) as src:
    print("CRS:", src.crs)
    print("Resolution:", src.res)
    print("Bounds:", src.bounds)