Sentinel-1 SLC Burst GCP 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!")

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

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

print("Dataset: S1C_IW_SLC__1SDV")
print("  Platform:      Sentinel-1C")
print("  Mode:          IW (Interferometric Wide)")
print("  Level:         L1 SLC")
print("  Polarizations: VV, VH")
print("  Date:          2025-10-16")

# Open burst IW1_VV_001 using the BURST open option
ds = gdal.OpenEx(
    zarr_path,
    gdal.OF_RASTER | gdal.OF_READONLY,
    open_options=["BURST=IW1_VV_001"],
)

assert ds is not None, "Failed to open SLC burst"

print("Burst:     IW1_VV_001")
print(f"Size:      {ds.RasterXSize} x {ds.RasterYSize} px")
print(f"Bands:     {ds.RasterCount}")
print(f"Data type: {gdal.GetDataTypeName(ds.GetRasterBand(1).DataType)}")
print(f"GCP count: {ds.GetGCPCount()}")

srs = ds.GetGCPSpatialRef()
if srs:
    print(f"GCP SRS:   {srs.GetName()} (EPSG:{srs.GetAuthorityCode(None)})")

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

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

gcps = ds.GetGCPs()
lons = [g.GCPX for g in gcps]
lats = [g.GCPY for g in gcps]
print(f"Lon range: [{min(lons):.4f}, {max(lons):.4f}]")
print(f"Lat range: [{min(lats):.4f}, {max(lats):.4f}]")

print("\nFirst 10 GCPs:")
print(f"  {'ID':<6} {'Pixel':>8} {'Line':>8} {'Lon':>10} {'Lat':>10} {'Height':>8}")
for gcp in gcps[:10]:
    print(
        f"  {gcp.Id:<6} {gcp.GCPPixel:>8.1f} {gcp.GCPLine:>8.1f} "
        f"{gcp.GCPX:>10.4f} {gcp.GCPY:>10.4f} {gcp.GCPZ:>8.1f}"
    )

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

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

# Image space
axes[0].scatter(pixels, lines, c="red", s=20, alpha=0.8)
axes[0].set_xlim(0, ds.RasterXSize)
axes[0].set_ylim(ds.RasterYSize, 0)
axes[0].set_xlabel("Pixel (column)")
axes[0].set_ylabel("Line (row)")
axes[0].set_title("GCP Positions — Image Space")
axes[0].set_aspect("auto")

# Geographic space
axes[1].scatter(lons, lats, c="blue", s=20, alpha=0.8)
axes[1].set_xlabel("Longitude (°)")
axes[1].set_ylabel("Latitude (°)")
axes[1].set_title("GCP Positions — Geographic Space (WGS 84)")
axes[1].grid(True, alpha=0.3)

plt.suptitle(f"SLC Burst IW1_VV_001 — {gcp_count} GCPs", fontweight="bold")
plt.tight_layout()
plt.show()

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

# Transform the four 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),
]

print(f"{'Corner':<14} {'Pixel':>8} {'Line':>8} {'Lon':>10} {'Lat':>10}")
print("-" * 55)
for name, px, ln in corners:
    ok, result = transformer.TransformPoint(0, px, ln, 0)
    if ok:
        print(f"{name:<14} {px:>8} {ln:>8} {result[0]:>10.4f} {result[1]:>10.4f}")

transformer = None

# Read the complex burst data and compute amplitude
raw = ds.GetRasterBand(1).ReadAsArray(buf_xsize=512, buf_ysize=0)
amplitude = np.abs(raw).astype(np.float32)

# Warp to WGS 84 using GCPs (reduced resolution to keep remote I/O fast)
warped = gdal.Warp(
    "",
    ds,
    format="MEM",
    dstSRS="EPSG:4326",
    width=512,
    height=0,
    resampleAlg="near",
    outputType=gdal.GDT_CFloat32,
)

assert warped is not None, "Warp failed"
gt = warped.GetGeoTransform()
print(f"Warped size:  {warped.RasterXSize} x {warped.RasterYSize} px")
print(f"Origin:       lon={gt[0]:.4f}  lat={gt[3]:.4f}")
print(f"Pixel size:   {gt[1]:.6f}° x {abs(gt[5]):.6f}°")

# Compute amplitude in dB for display
raw_warped = warped.GetRasterBand(1).ReadAsArray()
amp = np.abs(raw_warped).astype(np.float32)
amp_db = 20 * np.log10(np.where(amp > 0, amp, np.nan))

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

fig, ax = plt.subplots(figsize=(10, 6))
vmin, vmax = np.nanpercentile(amp_db, [2, 98])
im = ax.imshow(amp_db, cmap="gray", vmin=vmin, vmax=vmax, extent=extent, aspect="auto")
plt.colorbar(im, ax=ax, label="Amplitude (dB)")
ax.set_title(
    "Sentinel-1 SLC Burst IW1_VV_001 — GCP-warped to WGS 84", fontweight="bold"
)
ax.set_xlabel("Longitude (°)")
ax.set_ylabel("Latitude (°)")
ax.grid(True, alpha=0.3, color="yellow", linewidth=0.5)
plt.suptitle(
    f"Extent: lon [{extent[0]:.2f}, {extent[1]:.2f}]  lat [{extent[2]:.2f}, {extent[3]:.2f}]",
    fontsize=9,
)
plt.tight_layout()
plt.show()

warped = None

bursts = ["IW1_VV_001", "IW2_VV_001", "IW3_VV_001"]

print(f"{'Burst':<14} {'Size':^18} {'GCPs':>5}  {'Lon range':^20}  {'Lat range':^20}")
print("-" * 82)

burst_extents = {}
for burst_name in bursts:
    ds_b = gdal.OpenEx(zarr_path, gdal.OF_RASTER, open_options=[f"BURST={burst_name}"])
    if ds_b is None:
        print(f"{burst_name:<14}  NOT FOUND")
        continue
    gcps_b = ds_b.GetGCPs()
    n = len(gcps_b)
    lons_b = [g.GCPX for g in gcps_b]
    lats_b = [g.GCPY for g in gcps_b]
    size_str = f"{ds_b.RasterXSize}x{ds_b.RasterYSize}"
    lon_str = f"[{min(lons_b):.2f}, {max(lons_b):.2f}]"
    lat_str = f"[{min(lats_b):.2f}, {max(lats_b):.2f}]"
    print(f"{burst_name:<14} {size_str:^18} {n:>5}  {lon_str:^20}  {lat_str:^20}")
    burst_extents[burst_name] = (lons_b, lats_b)
    ds_b = None

# Plot geographic footprints of all three bursts
fig, ax = plt.subplots(figsize=(10, 6))
colors = {"IW1_VV_001": "blue", "IW2_VV_001": "orange", "IW3_VV_001": "green"}
for burst_name, (lons_b, lats_b) in burst_extents.items():
    ax.scatter(
        lons_b,
        lats_b,
        s=8,
        alpha=0.5,
        label=burst_name,
        color=colors.get(burst_name, "gray"),
    )

ax.set_xlabel("Longitude (°)")
ax.set_ylabel("Latitude (°)")
ax.set_title(
    "GCP Geographic Footprints — IW1, IW2, IW3 (VV burst 001)", fontweight="bold"
)
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

from scipy.spatial import ConvexHull

fig, ax = plt.subplots(figsize=(8, 7))
ax.set_facecolor("#c8d8e8")  # ocean fallback

try:
    import geopandas as gpd

    world = gpd.read_file(
        "https://naciscdn.org/naturalearth/110m/cultural/ne_110m_admin_0_countries.zip"
    )
    world.plot(ax=ax, color="#e8e0d0", edgecolor="gray", linewidth=0.5)
except Exception:
    pass  # plain ocean background is fine

colors = {
    "IW1_VV_001": "tab:blue",
    "IW2_VV_001": "tab:orange",
    "IW3_VV_001": "tab:green",
}

for burst_name, (lons_b, lats_b) in burst_extents.items():
    pts = np.column_stack([lons_b, lats_b])
    hull = ConvexHull(pts)
    hull_pts = pts[hull.vertices]
    hull_pts = np.vstack([hull_pts, hull_pts[0]])
    color = colors.get(burst_name, "gray")
    ax.fill(hull_pts[:, 0], hull_pts[:, 1], alpha=0.4, color=color, label=burst_name)
    ax.plot(hull_pts[:, 0], hull_pts[:, 1], color=color, linewidth=2)

# Auto-zoom to scene region with margin
all_lons = [lon for lons_b, _ in burst_extents.values() for lon in lons_b]
all_lats = [lat for _, lats_b in burst_extents.values() for lat in lats_b]
margin = 3
ax.set_xlim(min(all_lons) - margin, max(all_lons) + margin)
ax.set_ylim(min(all_lats) - margin, max(all_lats) + margin)

ax.set_xlabel("Longitude (°)")
ax.set_ylabel("Latitude (°)")
ax.set_title("Burst Footprints — IW1, IW2, IW3 (VV burst 001)", fontweight="bold")
ax.legend(loc="upper right")
ax.grid(True, alpha=0.3)

plt.suptitle(
    f"lon [{min(all_lons):.2f}, {max(all_lons):.2f}]  "
    f"lat [{min(all_lats):.2f}, {max(all_lats):.2f}]",
    fontsize=9,
)
plt.tight_layout()
plt.show()

print("Inspect GCPs on an SLC burst:\n")
print(f"  gdalinfo \"EOPFZARR:'{SLC_URL}'\" -oo BURST=IW1_VV_001\n")

print("Warp burst to WGS 84:\n")
print(
    f"  gdalwarp -t_srs EPSG:4326 -ts 512 0 \\\n"
    f"    \"EOPFZARR:'{SLC_URL}' -oo BURST=IW1_VV_001\" \\\n"
    f"    output_burst_wgs84.tif\n"
)

ds = None
print("All datasets closed")