Producing an Interferogram from Two NISAR GSLCs

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, 21)
end_date = datetime(2025, 12, 7)
area_of_interest = "POLYGON((40.1106 12.4104,41.8474 12.4104,41.8474 14.2538,40.1106 14.2538,40.1106 12.4104))" # 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,
    "relativeOrbit": 172,
    "frame": 8,
    "processingLevel": [
        "GSLC"
    ],
    "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)} GSLC products:")
hdf5_urls

import re

pattern = re.compile(
    r"(NISAR_L2_PR_GSLC.*?)_X([0-9]{5})_N_F_J_[0-9]{3}\.h5$"
)

latest_version_dict = {}

for url in hdf5_urls:
    fname = url.split("/")[-1]

    m = pattern.search(fname)
    if not m:
        print(f"No match: {fname}")
        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)} GSLC products:")
hdf5_urls

from pathlib import Path

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

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

gslc_files = sorted(data_dir.glob("*.h5"))
gslc_files

from getpass import getpass

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

%%time
import fsspec
import xarray as xr

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

group_path = "/science/LSAR/GSLC/grids/frequencyA"

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

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

import xarray as xr

gslc_dt = [
    xr.open_datatree(
        f,
        engine="h5netcdf",
        decode_timedelta=False,
        phony_dims="access",
        chunks="auto",
        group="/science/LSAR/GSLC/grids/frequencyA",
    )
    for f in gslc_files
]

gslc_dt[0]

gslc_dt[1]

pol = "HH"

ref_file = gslc_dt[0][pol]
sec_file = gslc_dt[1][pol]

#define epsg for subsetting and writing exports 

epsg = int(gslc_dt[0].projection.attrs["epsg_code"])

epsg

#define area of interest 

wkt = "POLYGON((40.5327 13.4179,40.8702 13.4179,40.8702 13.7066,40.5327 13.7066,40.5327 13.4179))"

import rioxarray
from shapely import wkt as shapely_wkt
from rasterio.warp import transform_bounds

epsg = int(gslc_dt[0].projection.attrs["epsg_code"])
crs = f"EPSG:{epsg}"

ref_file = ref_file.rio.write_crs(crs)
sec_file = sec_file.rio.write_crs(crs)

geom = shapely_wkt.loads(wkt)

bounds = transform_bounds(
    "EPSG:4326",
    crs,
    *geom.bounds
)

ref_subset = ref_file.rio.clip_box(*bounds)
sec_subset = sec_file.rio.clip_box(*bounds)

import matplotlib.pyplot as plt
from pyproj import Transformer
import contextily as cx

# Bounds from rioxarray
xmin, ymin, xmax, ymax = ref_file.rio.bounds()
axmin, aymin, axmax, aymax = ref_subset.rio.bounds()

transformer = Transformer.from_crs(
    ref_file.rio.crs,
    "EPSG:3857",
    always_xy=True
)

# Transform bounds
xmin_wm, ymin_wm = transformer.transform(xmin, ymin)
xmax_wm, ymax_wm = transformer.transform(xmax, ymax)

axmin_wm, aymin_wm = transformer.transform(axmin, aymin)
axmax_wm, aymax_wm = transformer.transform(axmax, aymax)

fig, ax = plt.subplots(figsize=(10, 10))

ax.plot(
    [xmin_wm, xmax_wm, xmax_wm, xmin_wm, xmin_wm],
    [ymin_wm, ymin_wm, ymax_wm, ymax_wm, ymin_wm],
    linewidth=2,
    color="yellow",
    label="Full raster"
)

ax.plot(
    [axmin_wm, axmax_wm, axmax_wm, axmin_wm, axmin_wm],
    [aymin_wm, aymin_wm, aymax_wm, aymax_wm, aymin_wm],
    linewidth=3,
    color="red",
    label="Subset AOI"
)

cx.add_basemap(
    ax,
    source=cx.providers.Esri.WorldImagery,
)

ax.legend(loc="upper right")
ax.set_axis_off()
plt.show()

import numpy as np

ifg_subset = ref_subset * np.conj(sec_subset)

ifg_subset

phase_subset = xr.apply_ufunc(
    np.angle,
    ifg_subset,
    dask="allowed",
)

phase_subset.plot.imshow(
    cmap="jet",
    vmin=-np.pi,
    vmax=np.pi,
    size=10,
    cbar_kwargs={"label": "Phase (radians)"}
)

plt.title("Wrapped Subsetted Interferogram")

ifg_scene = ref_file * np.conj(sec_file)

ifg_scene

phase_scene = xr.apply_ufunc(
    np.angle,
    ifg_scene,
    dask="allowed",
)

%%time
# multilook
phase_multilook = ifg_scene.coarsen(
    yCoordinates=10,
    xCoordinates=10,
    boundary="trim"
).mean()

# Find phase angles
phase_multilook_radians = xr.apply_ufunc(
    np.angle,
    phase_multilook,
    dask="allowed",
)

# plot
phase_multilook_radians.plot.imshow(
    cmap="jet",
    vmin=-np.pi,
    vmax=np.pi,
    size=10
)

phase_subset_export = phase_subset.rio.write_crs(f"EPSG:{epsg}")

phase_scene_export = phase_scene.rio.write_crs(f"EPSG:{epsg}")

phase_subset_export.rio.to_raster(
    "phase_subset.tif",
    compress="deflate"
)

phase_subset_export.to_netcdf(
    "phase_subset.h5",
    engine="h5netcdf"
)

%%time
phase_scene_export.rio.to_raster(
    "phase_scene.tif",
    driver="GTiff",
    tiled=True,
    windowed=True,
    BIGTIFF="YES",
)

%%time
# ~8MB chunks
phase_scene_h5 = phase_scene.chunk({
    "yCoordinates": 1440,
    "xCoordinates": 1440
})

phase_scene_h5.to_netcdf(
    "phase_scene.h5",
    engine="h5netcdf"
)