Basic Polarimetric Techniques with PolSARtools

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",
)

# Option 1: Use Earthdata Login credentials stored in a .netrc file
session = asf.ASFSession()

# Option 2: Uncomment if needed: Enter your NASA Earthdata Login credentials
#session.auth_with_creds(input('EDL Username'), getpass('EDL Password'))

start_date = datetime(2026, 7, 3)
end_date = datetime(2026, 7, 5)

area_of_interest = "POLYGON((-147.762 64.8592,-147.6954 64.8592,-147.6954 64.8789,-147.762 64.8789,-147.762 64.8592))" # POINT or POLYGON as WKT (well-known-text)

pattern = r'^(?!.*QA_STATS).*'

opts=asf.ASFSearchOptions(**{
    "maxResults": 250,
    "intersectsWith": area_of_interest,
    "flightDirection": "DESCENDING",
    "start": start_date,
    "end": end_date,
    "processingLevel": ["GSLC"],
    "dataset": ["NISAR"],
    "productionConfiguration": ["PR"],
    "dataMaturity" : "PROVISIONAL",
    '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

if len(hdf5_urls) == 1:
    print("Single GSLC product found. Skipping CRID filter.")
    
else:
    pattern = re.compile(
        r"(NISAR_L2_PR_GSLC.*?)_P([0-9]{5})_N_F_J_[0-9]{3}\.h5$")

    latest_version_dict = {}

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

        match = pattern.search(fname)

        if not match:
            print(f"No match: {fname}")
            continue

        product, crid = match.groups()

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

    hdf5_urls = [item[1] for item in latest_version_dict.values()]

print(f"Retained {len(hdf5_urls)} GSLC products:")
hdf5_urls

%%time 

from pathlib import Path

# Create a folder and download the files
data_dir = Path.home() / "polSAR_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

# Select the first downloaded GSLC file
in_file = str(gslc_files[0])

print(in_file)
print(type(in_file))

import polsartools as pst

pst.import_nisar_gslc(
    in_file,
    mat="C3",
    azlks=20,
    rglks=10,
    fmt="tif",
    cog=False,
    ovr=[2, 4, 8, 16],
    comp=False,
    recip=False
)

C3Folder = os.path.join(in_file.split(".h5")[0], "C3")

window_size = 3

pst.filter_refined_lee(
    C3Folder,
    win=window_size,
    fmt="tif",
    sub_dir=True,
    cog=False,
    ovr=[2, 4, 8, 16],
    comp=False
)

filtered_C3 = os.path.join(
    in_file.split(".h5")[0],
    f"rlee_{window_size}x{window_size}",
    "C3"
)

pst.pauli_rgb(filtered_C3, save_tif=True)

# Assign the HH, HV, and VV channels to the RGB channels
hh_path = next(Path(filtered_C3).glob("C11*.tif"))
hv_path = next(Path(filtered_C3).glob("C22*.tif"))
vv_path = next(Path(filtered_C3).glob("C33*.tif"))

# Create an RGB composite
pst.rgb(
    str(hh_path),  # Red: HH
    str(hv_path),  # Green: HV
    str(vv_path),  # Blue: VV
    norm_each=True,
    lp=5,
    gp=95,
    save_tif=True
)

%%time 
import sys
project_root = Path.cwd().parent
sys.path.insert(0, str(project_root))

from util.subset_h5 import subset_h5

area_of_interest = "POLYGON((-147.9245 64.7827,-147.779 64.7827,-147.779 64.8382,-147.9245 64.8382,-147.9245 64.7827))" # desired AOI for subset

subset_file = subset_h5(
    in_file,
    area_of_interest,
)

print(subset_file)

subset_in_file = str(subset_file)

pst.import_nisar_gslc(
    subset_in_file,
    mat="C3",
    azlks=2,
    rglks=2,
    fmt="tif",
    cog=False,
    ovr=[2, 4, 8, 16],
    comp=False,
    recip=False
)

subset_C3Folder = os.path.join(
    subset_in_file.split(".h5")[0],
    "C3"
)
subset_window_size = 3

pst.filter_refined_lee(
    subset_C3Folder,
    win=subset_window_size,
    fmt="tif",
    sub_dir=True,
    cog=False,
    ovr=[2, 4, 8, 16],
    comp=False
)

filtered_subset_C3 = os.path.join(
    subset_in_file.split(".h5")[0],
    f"rlee_{subset_window_size}x{subset_window_size}",
    "C3"
)

pst.pauli_rgb(filtered_subset_C3, save_tif=True)

# Assign the HH, HV, and VV channels to the RGB channels
hh_path = next(Path(filtered_subset_C3).glob("C11*.tif"))
hv_path = next(Path(filtered_subset_C3).glob("C22*.tif"))
vv_path = next(Path(filtered_subset_C3).glob("C33*.tif"))

# Create an HH/HV/VV RGB composite
pst.rgb(
    str(hh_path),  # Red: HH
    str(hv_path),  # Green: HV
    str(vv_path),  # Blue: VV
    norm_each=True,
    lp=5,
    gp=95,
    save_tif=True
)

# Scene-wide GeoTIFFs
scene_pauli_tif = Path(filtered_C3) / "PauliRGB.tif"
scene_rgb_tif = Path(filtered_C3) / "RGB.tif"

# High res subsetted GeoTIFFs
subset_pauli_tif = Path(filtered_subset_C3) / "PauliRGB.tif"
subset_rgb_tif = Path(filtered_subset_C3) / "RGB.tif"

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

def prepare_rgb_overlay(tif_path):

    with rasterio.open(tif_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs,
            "EPSG:4326",
            src.width,
            src.height,
            *src.bounds,
        )

        rgb = np.zeros((3, height, width), dtype=np.uint8)

        for band in range(3):
            reproject(
                source=rasterio.band(src, band + 1),
                destination=rgb[band],
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs="EPSG:4326",
                resampling=Resampling.nearest,
            )

    rgb = np.moveaxis(rgb, 0, -1)

    # Make pixels with no RGB values transparent
    alpha = np.where(np.any(rgb > 0, axis=2), 255, 0).astype(np.uint8)
    rgba = np.dstack((rgb, alpha))

    west = transform.c
    north = transform.f
    east = west + transform.a * width
    south = north + transform.e * height

    bounds = [[south, west], [north, east]]

    return rgba, bounds

scene_pauli_img, scene_pauli_bounds = prepare_rgb_overlay(scene_pauli_tif)
scene_rgb_img, scene_rgb_bounds = prepare_rgb_overlay(scene_rgb_tif)

subset_pauli_img, subset_pauli_bounds = prepare_rgb_overlay(subset_pauli_tif)
subset_rgb_img, subset_rgb_bounds = prepare_rgb_overlay(subset_rgb_tif)

import folium

m = folium.Map(
    tiles="Esri.WorldImagery",
)

folium.raster_layers.ImageOverlay(
    image=scene_pauli_img,
    bounds=scene_pauli_bounds,
    name="Scene-wide Pauli RGB",
    opacity=0.8,
    show=True,
).add_to(m)

folium.raster_layers.ImageOverlay(
    image=scene_rgb_img,
    bounds=scene_rgb_bounds,
    name="Scene-wide HH/HV/VV RGB",
    opacity=0.8,
    show=False,
).add_to(m)

folium.raster_layers.ImageOverlay(
    image=subset_pauli_img,
    bounds=subset_pauli_bounds,
    name="Subset Pauli RGB",
    opacity=0.9,
    show=False,
).add_to(m)

folium.raster_layers.ImageOverlay(
    image=subset_rgb_img,
    bounds=subset_rgb_bounds,
    name="Subset HH/HV/VV RGB",
    opacity=0.9,
    show=False,
).add_to(m)

folium.LayerControl(collapsed=False).add_to(m)

m.fit_bounds(scene_pauli_bounds)

m