Sentinel-1 GRD Geolocation with EOPFZARR Driver

import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal, osr

gdal.UseExceptions()

print(f"GDAL version: {gdal.__version__}")

drv = gdal.GetDriverByName("EOPFZARR")
if drv:
    print("EOPFZARR driver: Registered")
else:
    print("WARNING: EOPFZARR driver not found!")

GRD_URL = (
    "/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:"
    "notebook-data/tutorial_data/cpm_v264/"
    "S1C_IW_GRDH_1SDV_20260205T120122_20260205T120158_006220_00C7E4_5D6E.zarr"
)

zarr_path = f"EOPFZARR:'{GRD_URL}'"

print("Dataset: S1C_IW_GRDH_1SDV")
print("  Platform:      Sentinel-1C")
print("  Mode:          IW (Interferometric Wide)")
print("  Resolution:    High")
print("  Polarization:  Dual VV/VH")

# Open product in standard mode to see subdatasets
root_ds = gdal.OpenEx(zarr_path, gdal.OF_RASTER, open_options=["GRD_MULTIBAND=NO"])

subdatasets = root_ds.GetMetadata("SUBDATASETS")

# Find the VV GRD measurement subdataset
vv_path = None
for key, value in sorted(subdatasets.items()):
    if "_NAME" in key and "measurements/grd" in value and "_VV" in value:
        vv_path = value
        break

print(f"VV subdataset: {vv_path.split(':/')[1] if vv_path else 'not found'}")

# Open the VV subdataset
ds = gdal.Open(vv_path)
print(f"\nDimensions: {ds.RasterXSize} x {ds.RasterYSize} pixels")
print(f"Data Type:  {gdal.GetDataTypeName(ds.GetRasterBand(1).DataType)}")

root_ds = None

gcp_count = ds.GetGCPCount()
print(f"GCP Count: {gcp_count}")

# Check the GCP spatial reference
gcp_srs = ds.GetGCPSpatialRef()
if gcp_srs:
    print(f"GCP SRS:   {gcp_srs.GetName()} (EPSG:{gcp_srs.GetAuthorityCode(None)})")

# Display the GCPs
gcps = ds.GetGCPs()
print("\nFirst 10 GCPs:")
print(
    f"{'ID':>6}  {'Pixel':>10}  {'Line':>10}  {'Longitude':>12}  {'Latitude':>12}  {'Height':>10}"
)
print("-" * 72)
for gcp in gcps[:10]:
    print(
        f"{gcp.Id:>6}  {gcp.GCPPixel:>10.1f}  {gcp.GCPLine:>10.1f}  "
        f"{gcp.GCPX:>12.4f}  {gcp.GCPY:>12.4f}  {gcp.GCPZ:>10.1f}"
    )

# Compute coverage
lons = [gcp.GCPX for gcp in gcps]
lats = [gcp.GCPY for gcp in gcps]
print("\nGeographic Extent:")
print(f"  Longitude: {min(lons):.4f} to {max(lons):.4f}")
print(f"  Latitude:  {min(lats):.4f} to {max(lats):.4f}")

pixels = [gcp.GCPPixel for gcp in gcps]
lines = [gcp.GCPLine for gcp in gcps]

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Plot 1: GCP positions in pixel space
axes[0].scatter(pixels, lines, c="red", s=15, alpha=0.7)
axes[0].set_xlim(0, ds.RasterXSize)
axes[0].set_ylim(ds.RasterYSize, 0)
axes[0].set_title("GCP Positions in Image Space", fontweight="bold")
axes[0].set_xlabel("Pixel")
axes[0].set_ylabel("Line")
axes[0].set_aspect("equal")

# Plot 2: GCP positions in geographic space
scatter = axes[1].scatter(
    lons, lats, c=range(len(lons)), cmap="viridis", s=15, alpha=0.7
)
axes[1].set_title("GCP Positions in Geographic Space", fontweight="bold")
axes[1].set_xlabel("Longitude")
axes[1].set_ylabel("Latitude")
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Compute grid dimensions
unique_pixels = len(set(int(p) for p in pixels))
unique_lines = len(set(int(ln) for ln in lines))
print(f"GCP grid: {unique_pixels} columns x {unique_lines} rows = {gcp_count} points")
print(
    f"Sampling: ~1 GCP every {ds.RasterXSize // unique_pixels} pixels (range)"
    f" x {ds.RasterYSize // unique_lines} pixels (azimuth)"
)

# Create a GCP-based transformer
transformer = gdal.Transformer(ds, None, ["METHOD=GCP_POLYNOMIAL", "MAX_GCP_ORDER=3"])

# Transform image corners to geographic coordinates
corners = [
    ("Top-Left", 0, 0),
    ("Top-Right", ds.RasterXSize, 0),
    ("Bottom-Left", 0, ds.RasterYSize),
    ("Bottom-Right", ds.RasterXSize, ds.RasterYSize),
    ("Center", ds.RasterXSize // 2, ds.RasterYSize // 2),
]

print("Pixel-to-Geographic Coordinate Transformation:")
print(f"{'Location':<15} {'Pixel':>8} {'Line':>8}   {'Longitude':>12} {'Latitude':>12}")
print("-" * 65)
for name, px, ln in corners:
    success, result = transformer.TransformPoint(False, px, ln)
    if success:
        print(f"{name:<15} {px:>8} {ln:>8}   {result[0]:>12.4f} {result[1]:>12.4f}")
    else:
        print(f"{name:<15} {px:>8} {ln:>8}   transform failed")

# Warp a small subset to WGS 84 using GCPs
# We use a reduced resolution to keep it fast over /vsicurl/
warped = gdal.Warp(
    "",
    ds,
    format="MEM",
    dstSRS="EPSG:4326",
    width=1024,
    height=0,  # auto-compute from aspect ratio
    resampleAlg="near",
)

if warped:
    gt = warped.GetGeoTransform()
    print(f"Warped dataset: {warped.RasterXSize} x {warped.RasterYSize} pixels")
    print(
        f"GeoTransform:   origin=({gt[0]:.4f}, {gt[3]:.4f}), "
        f"pixel size=({gt[1]:.6f}, {gt[5]:.6f})"
    )
    print(f"Projection:     {warped.GetProjection()[:60]}...")

    try:
        data = warped.GetRasterBand(1).ReadAsArray()
        data_float = data.astype(np.float32)
        data_float[data_float == 0] = np.nan
        data_db = 10 * np.log10(data_float + 1)

        # Compute extent from geotransform
        extent = [
            gt[0],
            gt[0] + gt[1] * warped.RasterXSize,
            gt[3] + gt[5] * warped.RasterYSize,
            gt[3],
        ]

        fig, axes = plt.subplots(1, 2, figsize=(14, 6))

        # Original (pixel space)
        orig_sub = ds.GetRasterBand(1).ReadAsArray(
            buf_xsize=1024,
            buf_ysize=int(1024 * ds.RasterYSize / ds.RasterXSize),
        )
        orig_float = orig_sub.astype(np.float32)
        orig_float[orig_float == 0] = np.nan
        orig_db = 10 * np.log10(orig_float + 1)
        vmin, vmax = np.nanpercentile(orig_db, [2, 98])
        axes[0].imshow(orig_db, cmap="gray", vmin=vmin, vmax=vmax)
        axes[0].set_title("Original (Image Space)", fontweight="bold")
        axes[0].set_xlabel("Pixel")
        axes[0].set_ylabel("Line")

        # Warped (geographic space)
        vmin, vmax = np.nanpercentile(data_db, [2, 98])
        axes[1].imshow(data_db, cmap="gray", vmin=vmin, vmax=vmax, extent=extent)
        axes[1].set_title("Warped to WGS 84 (using GCPs)", fontweight="bold")
        axes[1].set_xlabel("Longitude")
        axes[1].set_ylabel("Latitude")

        plt.tight_layout()
        plt.show()
    except ImportError:
        print("ReadAsArray not available (NumPy version mismatch), skipping plot")

    warped = None
else:
    print("Warp failed")

SLC_URL = (
    "/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:"
    "notebook-data/tutorial_data/cpm_v262/"
    "S1C_IW_SLC__1SDV_20251016T165627_20251016T165654_004590_00913B_30C4.zarr"
)

ds_slc = gdal.OpenEx(
    f"EOPFZARR:'{SLC_URL}'", gdal.OF_RASTER, open_options=["BURST=IW1_VV_001"]
)

if ds_slc:
    slc_gcp_count = ds_slc.GetGCPCount()
    slc_srs = ds_slc.GetGCPSpatialRef()
    print("SLC Burst IW1_VV_001:")
    print(f"  Dimensions: {ds_slc.RasterXSize} x {ds_slc.RasterYSize}")
    print(f"  GCP Count:  {slc_gcp_count}")
    if slc_srs:
        print(f"  GCP SRS:    {slc_srs.GetName()}")

    if slc_gcp_count > 0:
        slc_gcps = ds_slc.GetGCPs()
        slc_lons = [g.GCPX for g in slc_gcps]
        slc_lats = [g.GCPY for g in slc_gcps]
        print(
            f"  Extent:     lon [{min(slc_lons):.4f}, {max(slc_lons):.4f}]"
            f"  lat [{min(slc_lats):.4f}, {max(slc_lats):.4f}]"
        )

    ds_slc = None
else:
    print("Failed to open SLC burst")

# Print the actual subdataset path used in this notebook
subds_path = vv_path.split("'")[1] if "'" in vv_path else vv_path

print("Step 1: List subdatasets\n")
print(f"  gdalinfo \"EOPFZARR:'{GRD_URL}'\" -oo GRD_MULTIBAND=NO\n")
print("Step 2: Open a measurement subdataset to see GCPs\n")
print(f"  gdalinfo '{subds_path}'\n")
print("The output will include:")
print('  GCP Projection = GEOGCS["WGS 84", ...]')
print("  GCP[  0]: Id=GCP_1, Info=")
print("            (0.5,0.5) -> (-91.4153, 16.2326, 574.99)")
print("  GCP[  1]: Id=GCP_2, Info=")
print("            (1258.5,0.5) -> (-91.5360, 16.2558, 904.98)")
print("  ...")

ds = None
print("All datasets closed")