Sentinel-3 Multi-Product Demo with GDAL

# Import required libraries
from osgeo import gdal
import matplotlib.pyplot as plt
import numpy as np

# Note: xarray import skipped due to version compatibility issues
# We use pure numpy/GDAL approach instead
import warnings

warnings.filterwarnings("ignore")

# Enable GDAL exceptions
gdal.UseExceptions()

# Verify EOPFZARR driver is available
driver = gdal.GetDriverByName("EOPFZARR")
if driver:
    print(f"✅ EOPFZARR driver loaded: {driver.LongName}")
else:
    print("❌ EOPFZARR driver not found!")

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

# Define all Sentinel-3 product URLs
PRODUCTS = {
    "OLCI_L1_EFR": {
        "name": "Sentinel-3 OLCI Level-1 EFR (Full Resolution)",
        "url": "https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/"
        "cpm_v264/S3B_OL_1_EFR____20250819T074058_20250819T074358_20250820T082937_0179_110_106_3420_ESA_O_NT_004.zarr",
        "description": "TOA Radiances - 21 spectral bands (400-1020nm)",
        "sample_band": "measurements/oa01_radiance",
    },
    "OLCI_L1_ERR": {
        "name": "Sentinel-3 OLCI Level-1 ERR (Reduced Resolution)",
        "url": "https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/"
        "cpm_v264/S3A_OL_1_ERR____20251019T145533_20251019T153950_20251020T155912_2657_131_353______PS1_O_NT_004.zarr",
        "description": "TOA Radiances - Reduced resolution (1.2km)",
        "sample_band": "measurements/oa01_radiance",
    },
    "OLCI_L2_LFR": {
        "name": "Sentinel-3 OLCI Level-2 LFR (Land Full Resolution)",
        "url": "https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/"
        "cpm_v264/S3A_OL_2_LFR____20250817T085058_20250817T085358_20250818T125249_0179_129_221_2160_PS1_O_NT_003.zarr",
        "description": "Land products - OGVI, OTCI, RC681, etc.",
        "sample_band": None,  # Will discover
    },
    "SLSTR_L1_RBT": {
        "name": "Sentinel-3 SLSTR Level-1 RBT",
        "url": "https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/"
        "cpm_v264/S3A_SL_1_RBT____20251019T064521_20251019T064821_20251020T175322_0179_131_348_2700_PS1_O_NT_004.zarr",
        "description": "Radiances and Brightness Temperatures",
        "sample_band": None,  # Will discover
    },
    "SLSTR_L2_LST": {
        "name": "Sentinel-3 SLSTR Level-2 LST",
        "url": "https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/"
        "cpm_v264/S3A_SL_2_LST____20251204T100846_20251204T101146_20251205T211618_0179_133_236_2340_PS1_O_NT_004.zarr",
        "description": "Land Surface Temperature",
        "sample_band": None,  # Will discover
    },
}

# Helper function for eopfzarr paths
def eopfzarr_path(url, subdataset=None):
    """Construct EOPFZARR path for remote access

    Note: The subdataset path MUST start with '/' for GEOLOCATION metadata to work.
    This is required by the EOPFZARR driver for proper geolocation array handling.
    """
    base = f'EOPFZARR:"/vsicurl/{url}"'
    if subdataset:
        # Ensure subdataset starts with '/' for geolocation arrays
        if not subdataset.startswith("/"):
            subdataset = "/" + subdataset
        return f"{base}:{subdataset}"
    return base


print(f"Configured {len(PRODUCTS)} Sentinel-3 products for analysis")

def explore_product(product_key, max_subdatasets=15):
    """Explore a Sentinel-3 product and list its subdatasets"""
    product = PRODUCTS[product_key]
    print(f"\n{'='*80}")
    print(f"📦 {product['name']}")
    print(f"{'='*80}")
    print(f"Description: {product['description']}")
    print(f"URL: ...{product['url'][-60:]}")

    try:
        ds = gdal.Open(eopfzarr_path(product["url"]))
        if ds is None:
            print("❌ Failed to open dataset")
            return None

        print(f"\nDriver: {ds.GetDriver().ShortName}")
        print(f"Size: {ds.RasterXSize} x {ds.RasterYSize}")

        # Get projection info
        proj = ds.GetProjection()
        if proj:
            # Extract EPSG if available
            if "EPSG" in proj:
                import re

                epsg_match = re.search(r'EPSG","?(\d+)', proj)
                if epsg_match:
                    print(f"CRS: EPSG:{epsg_match.group(1)}")

        # Get subdatasets
        subdatasets = ds.GetMetadata("SUBDATASETS")
        subdataset_list = []

        for key, value in subdatasets.items():
            if key.startswith("SUBDATASET_") and "_NAME" in key:
                path = value.split(":")[-1]
                subdataset_list.append(path)

        print(f"\nSubdatasets: {len(subdataset_list)} total")
        print("\nSample subdatasets:")
        for i, sd in enumerate(subdataset_list[:max_subdatasets]):
            print(f"  {i+1:2d}. {sd}")
        if len(subdataset_list) > max_subdatasets:
            print(f"  ... and {len(subdataset_list) - max_subdatasets} more")

        ds = None
        return subdataset_list

    except Exception as e:
        print(f"❌ Error: {e}")
        return None


# Explore all products
product_subdatasets = {}
for key in PRODUCTS.keys():
    product_subdatasets[key] = explore_product(key)

def read_band_data(url, subdataset, size=None):
    """Read band data from a subdataset"""
    path = eopfzarr_path(url, subdataset)
    ds = gdal.Open(path)
    if ds is None:
        raise ValueError(f"Cannot open {subdataset}")

    band = ds.GetRasterBand(1)

    if size is None:
        data = band.ReadAsArray().astype(np.float32)
    else:
        # Read a subset
        width = min(size, ds.RasterXSize)
        height = min(size, ds.RasterYSize)
        data = band.ReadAsArray(0, 0, width, height).astype(np.float32)

    nodata = band.GetNoDataValue()
    if nodata is not None:
        data[data == nodata] = np.nan

    metadata = {
        "width": ds.RasterXSize,
        "height": ds.RasterYSize,
        "nodata": nodata,
        "dtype": band.DataType,
    }

    ds = None
    return data, metadata


# Read OLCI L1 EFR bands
print("Reading OLCI L1 EFR spectral bands...")

olci_bands = {
    "Oa01 (400nm)": "measurements/oa01_radiance",
    "Oa06 (560nm)": "measurements/oa06_radiance",
    "Oa08 (665nm)": "measurements/oa08_radiance",
    "Oa17 (865nm)": "measurements/oa17_radiance",
}

olci_data = {}
for name, subdataset in olci_bands.items():
    try:
        data, meta = read_band_data(
            PRODUCTS["OLCI_L1_EFR"]["url"], subdataset, size=1000
        )
        olci_data[name] = data
        print(
            f"  ✅ {name}: {data.shape}, range [{np.nanmin(data):.2f}, {np.nanmax(data):.2f}]"
        )
    except Exception as e:
        print(f"  ❌ {name}: {e}")

# Visualize OLCI L1 EFR bands
if olci_data:
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    axes = axes.flatten()

    for idx, (name, data) in enumerate(olci_data.items()):
        # Use robust percentile scaling
        vmin, vmax = np.nanpercentile(data, [2, 98])
        im = axes[idx].imshow(data, cmap="viridis", vmin=vmin, vmax=vmax)
        axes[idx].set_title(f"OLCI L1 EFR - {name}", fontsize=11, fontweight="bold")
        axes[idx].set_xlabel("Column")
        axes[idx].set_ylabel("Row")
        plt.colorbar(im, ax=axes[idx], label="Radiance", fraction=0.046)

    plt.suptitle(
        "Sentinel-3 OLCI L1 EFR - Multi-Spectral Bands", fontsize=14, fontweight="bold"
    )
    plt.tight_layout()
    plt.show()

    print("✅ OLCI L1 EFR visualization complete!")

# Read OLCI L1 ERR bands
print("Reading OLCI L1 ERR (Reduced Resolution) spectral bands...")

err_data = {}
for name, subdataset in olci_bands.items():
    try:
        data, meta = read_band_data(
            PRODUCTS["OLCI_L1_ERR"]["url"], subdataset, size=800
        )
        err_data[name] = data
        print(
            f"  ✅ {name}: {data.shape}, range [{np.nanmin(data):.2f}, {np.nanmax(data):.2f}]"
        )
    except Exception as e:
        print(f"  ❌ {name}: {e}")

# Visualize
if err_data:
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    axes = axes.flatten()

    for idx, (name, data) in enumerate(err_data.items()):
        vmin, vmax = np.nanpercentile(data, [2, 98])
        im = axes[idx].imshow(data, cmap="viridis", vmin=vmin, vmax=vmax)
        axes[idx].set_title(f"OLCI L1 ERR - {name}", fontsize=11, fontweight="bold")
        axes[idx].set_xlabel("Column")
        axes[idx].set_ylabel("Row")
        plt.colorbar(im, ax=axes[idx], label="Radiance", fraction=0.046)

    plt.suptitle(
        "Sentinel-3 OLCI L1 ERR - Reduced Resolution", fontsize=14, fontweight="bold"
    )
    plt.tight_layout()
    plt.show()

    print("✅ OLCI L1 ERR visualization complete!")

# Discover L2 LFR subdatasets
print("Discovering OLCI L2 LFR subdatasets...")

lfr_subdatasets = product_subdatasets.get("OLCI_L2_LFR", [])

# Find measurement subdatasets
lfr_measurements = [
    sd
    for sd in lfr_subdatasets
    if "measurements/" in sd and "latitude" not in sd and "longitude" not in sd
]

print(f"\nAvailable L2 LFR measurements ({len(lfr_measurements)}):")
for sd in lfr_measurements[:15]:
    print(f"  {sd}")

# Read L2 LFR products - try common land products
print("Reading OLCI L2 LFR land products...")

# Common L2 land product names to try
l2_products_to_try = [
    ("GIFAPAR", "measurements/gifapar"),  # OLCI Global Vegetation Index
    ("OTCI", "measurements/otci"),  # OLCI Terrestrial Chlorophyll Index
    ("IWV", "measurements/iwv"),  # Integrated Water Vapour
    ("RC681", "measurements/rc681"),  # Rectified Reflectance 681nm
]

lfr_data = {}
for name, subdataset in l2_products_to_try:
    try:
        data, meta = read_band_data(
            PRODUCTS["OLCI_L2_LFR"]["url"], subdataset, size=800
        )
        lfr_data[name] = data
        print(
            f"  ✅ {name}: {data.shape}, range [{np.nanmin(data):.4f}, {np.nanmax(data):.4f}]"
        )
    except Exception as e:
        print(f"  ⚠️  {name}: Not available or error - {str(e)[:50]}")

# If we found products, visualize them
if lfr_data:
    n_products = len(lfr_data)
    cols = min(2, n_products)
    rows = (n_products + 1) // 2

    fig, axes = plt.subplots(rows, cols, figsize=(7 * cols, 5 * rows))
    if n_products == 1:
        axes = [axes]
    else:
        axes = axes.flatten()

    for idx, (name, data) in enumerate(lfr_data.items()):
        vmin, vmax = np.nanpercentile(data, [2, 98])

        # Choose colormap based on product
        if "OGVI" in name or "OTCI" in name:
            cmap = "RdYlGn"
        else:
            cmap = "viridis"

        im = axes[idx].imshow(data, cmap=cmap, vmin=vmin, vmax=vmax)
        axes[idx].set_title(f"OLCI L2 LFR - {name}", fontsize=11, fontweight="bold")
        axes[idx].set_xlabel("Column")
        axes[idx].set_ylabel("Row")
        plt.colorbar(im, ax=axes[idx], fraction=0.046)

    # Hide unused subplots
    for idx in range(len(lfr_data), len(axes)):
        axes[idx].axis("off")

    plt.suptitle(
        "Sentinel-3 OLCI L2 LFR - Land Products", fontsize=14, fontweight="bold"
    )
    plt.tight_layout()
    plt.show()
else:
    print("⚠️  No L2 LFR products loaded. Checking available subdatasets...")
    print(f"Available: {lfr_measurements[:5]}")

# Discover SLSTR L1 RBT subdatasets
print("Discovering SLSTR L1 RBT subdatasets...")

rbt_subdatasets = product_subdatasets.get("SLSTR_L1_RBT", [])

# Find measurement subdatasets
rbt_measurements = [sd for sd in rbt_subdatasets if "measurements/" in sd]

print(f"\nAvailable SLSTR L1 RBT measurements ({len(rbt_measurements)}):")
for sd in rbt_measurements[:20]:
    print(f"  {sd}")

# Read SLSTR L1 RBT bands
print("Reading SLSTR L1 RBT bands...")

# SLSTR bands to try (common naming patterns)
slstr_bands_to_try = [
    ("S1 (555nm)", "measurements/anadir/s1_radiance_an"),
    ("S2 (659nm)", "measurements/anadir/s2_radiance_an"),
    ("S3 (865nm)", "measurements/anadir/s3_radiance_an"),
    ("S7 (3.7μm BT)", "measurements/inadir/s7_bt_in"),
    ("S8 (10.8μm BT)", "measurements/inadir/s8_bt_in"),
    ("S9 (12μm BT)", "measurements/inadir/s9_bt_in"),
]

rbt_data = {}
for name, subdataset in slstr_bands_to_try:
    try:
        data, meta = read_band_data(
            PRODUCTS["SLSTR_L1_RBT"]["url"], subdataset, size=600
        )
        rbt_data[name] = data
        print(
            f"  ✅ {name}: {data.shape}, range [{np.nanmin(data):.2f}, {np.nanmax(data):.2f}]"
        )
    except Exception as e:
        print(f"  ⚠️  {name}: {str(e)[:60]}")

# Visualize
if rbt_data:
    n_bands = len(rbt_data)
    cols = min(3, n_bands)
    rows = (n_bands + 2) // 3

    fig, axes = plt.subplots(rows, cols, figsize=(5 * cols, 4 * rows))
    if n_bands == 1:
        axes = [axes]
    else:
        axes = axes.flatten()

    for idx, (name, data) in enumerate(rbt_data.items()):
        vmin, vmax = np.nanpercentile(data, [2, 98])

        # Use thermal colormap for brightness temperature bands
        if "BT" in name or "μm" in name:
            cmap = "inferno"
            label = "Brightness Temp (K)"
        else:
            cmap = "viridis"
            label = "Radiance"

        im = axes[idx].imshow(data, cmap=cmap, vmin=vmin, vmax=vmax)
        axes[idx].set_title(f"{name}", fontsize=10, fontweight="bold")
        plt.colorbar(im, ax=axes[idx], label=label, fraction=0.046)

    # Hide unused
    for idx in range(len(rbt_data), len(axes)):
        axes[idx].axis("off")

    plt.suptitle(
        "Sentinel-3 SLSTR L1 RBT - Radiances & Brightness Temperatures",
        fontsize=13,
        fontweight="bold",
    )
    plt.tight_layout()
    plt.show()
else:
    print("⚠️  No SLSTR L1 RBT bands loaded")

# Discover SLSTR L2 LST subdatasets
print("Discovering SLSTR L2 LST subdatasets...")

lst_subdatasets = product_subdatasets.get("SLSTR_L2_LST", [])

# Find measurement subdatasets
lst_measurements = [sd for sd in lst_subdatasets if "measurements/" in sd]

print(f"\nAvailable SLSTR L2 LST measurements ({len(lst_measurements)}):")
for sd in lst_measurements[:15]:
    print(f"  {sd}")

lst_conditions = [sd for sd in lst_subdatasets if "conditions/" in sd]
print(f"\nAvailable SLSTR L2 LST conditions ({len(lst_conditions)}):")

for sd in lst_conditions[:15]:
    print(f"  {sd}")

# Read SLSTR L2 LST data
print("Reading SLSTR L2 LST products...")

lst_products_to_try = [
    ("LST", "measurements/lst"),
    ("LST Uncertainty", "quality/lst_uncertainty"),
    ("NDVI", "conditions/auxiliary/ndvi"),
    ("Biome", "conditions/auxiliary/biome"),
]

lst_data = {}
for name, subdataset in lst_products_to_try:
    try:
        data, meta = read_band_data(
            PRODUCTS["SLSTR_L2_LST"]["url"], subdataset, size=600
        )
        lst_data[name] = data
        print(
            f"  ✅ {name}: {data.shape}, range [{np.nanmin(data):.2f}, {np.nanmax(data):.2f}]"
        )
    except Exception as e:
        print(f"  ⚠️  {name}: {str(e)[:60]}")

# Visualize LST
if lst_data:
    n_products = len(lst_data)
    cols = min(2, n_products)
    rows = (n_products + 1) // 2

    fig, axes = plt.subplots(rows, cols, figsize=(7 * cols, 5 * rows))
    if n_products == 1:
        axes = [axes]
    else:
        axes = axes.flatten()

    for idx, (name, data) in enumerate(lst_data.items()):
        vmin, vmax = np.nanpercentile(data, [2, 98])

        # Choose appropriate colormap
        if "LST" in name and "Uncertainty" not in name:
            cmap = "RdYlBu_r"  # Temperature colormap
            label = "Temperature (K)"
        elif "NDVI" in name:
            cmap = "RdYlGn"
            label = "NDVI"
        else:
            cmap = "viridis"
            label = "Value"

        im = axes[idx].imshow(data, cmap=cmap, vmin=vmin, vmax=vmax)
        axes[idx].set_title(f"SLSTR L2 LST - {name}", fontsize=11, fontweight="bold")
        axes[idx].set_xlabel("Column")
        axes[idx].set_ylabel("Row")
        plt.colorbar(im, ax=axes[idx], label=label, fraction=0.046)

    # Hide unused
    for idx in range(len(lst_data), len(axes)):
        axes[idx].axis("off")

    plt.suptitle(
        "Sentinel-3 SLSTR L2 LST - Land Surface Temperature",
        fontsize=14,
        fontweight="bold",
    )
    plt.tight_layout()
    plt.show()
else:
    print("⚠️  No SLSTR L2 LST products loaded")

# Demonstrate EOPFZARR Geolocation Arrays Support for ALL Products
print("=" * 80)
print("EOPFZARR GEOLOCATION ARRAYS DEMO - All Sentinel-3 Products")
print("=" * 80)

# Define sample subdatasets for each product type
product_samples = {
    "OLCI_L1_EFR": "measurements/oa01_radiance",
    "OLCI_L1_ERR": "measurements/oa01_radiance",
    "OLCI_L2_LFR": "measurements/gifapar",
    "SLSTR_L1_RBT": "measurements/inadir/s8_bt_in",
    "SLSTR_L2_LST": "measurements/lst",
}

geoloc_results = {}

for product_key, subdataset in product_samples.items():
    url = PRODUCTS[product_key]["url"]
    path = eopfzarr_path(url, subdataset)

    print(f"\n{'─'*70}")
    print(f"📦 {product_key}")
    print(f"   Subdataset: {subdataset}")

    try:
        ds = gdal.Open(path)
        if ds:
            geoloc = ds.GetMetadata("GEOLOCATION")

            if geoloc:
                x_path = geoloc.get("X_DATASET", "N/A")
                y_path = geoloc.get("Y_DATASET", "N/A")
                srs = geoloc.get("SRS", "N/A")

                # Extract just the measurement path
                x_meas = x_path.split(":")[-1] if x_path else "N/A"
                y_meas = y_path.split(":")[-1] if y_path else "N/A"

                print("   ✅ GEOLOCATION metadata found!")
                print(f"      X_DATASET (Lon): {x_meas}")
                print(f"      Y_DATASET (Lat): {y_meas}")
                print(f"      PIXEL_STEP: {geoloc.get('PIXEL_STEP', 'N/A')}")
                print(f"      LINE_STEP: {geoloc.get('LINE_STEP', 'N/A')}")

                # Check if coordinates are WGS84
                if "WGS" in srs or "EPSG" in srs:
                    print("      SRS: WGS 84 (EPSG:4326)")

                geoloc_results[product_key] = {
                    "status": "success",
                    "x_path": x_path,
                    "y_path": y_path,
                    "size": f"{ds.RasterXSize}x{ds.RasterYSize}",
                    "subdataset": subdataset,
                }
            else:
                print("   ⚠️ No GEOLOCATION metadata")
                geoloc_results[product_key] = {"status": "no_geoloc"}

            ds = None
        else:
            print("   ❌ Failed to open dataset")
            geoloc_results[product_key] = {"status": "failed"}

    except Exception as e:
        print(f"   ❌ Error: {str(e)[:50]}")
        geoloc_results[product_key] = {"status": "error", "msg": str(e)}

# Summary
print(f"\n{'='*80}")
print("GEOLOCATION ARRAYS SUMMARY")
print("=" * 80)
success_count = sum(1 for r in geoloc_results.values() if r["status"] == "success")
print(f"\n✅ {success_count}/{len(geoloc_results)} products have GEOLOCATION arrays\n")

print(f"{'Product':<20} {'Status':<12} {'Lon/Lat Arrays'}")
print("-" * 70)
for product, result in geoloc_results.items():
    if result["status"] == "success":
        lon_arr = result["x_path"].split("/")[-1]
        lat_arr = result["y_path"].split("/")[-1]
        print(f"{product:<20} {'✅ Found':<12} {lon_arr} / {lat_arr}")
    else:
        print(f"{product:<20} {'❌ Missing':<12} -")

# Visualize geolocation for ALL 5 products
print("\n" + "=" * 80)
print("GEOGRAPHIC VISUALIZATION - All 5 Sentinel-3 Products")
print("=" * 80)

# Create a 5x2 grid for all products (Array View | Geographic View)
fig, axes = plt.subplots(5, 2, figsize=(14, 25))

# All products to visualize
viz_products = [
    (
        "OLCI_L1_EFR",
        "measurements/oa01_radiance",
        "OLCI L1 EFR (Oa01 Radiance)",
        "viridis",
    ),
    (
        "OLCI_L1_ERR",
        "measurements/oa01_radiance",
        "OLCI L1 ERR (Oa01 Radiance)",
        "viridis",
    ),
    ("OLCI_L2_LFR", "measurements/gifapar", "OLCI L2 LFR (GIFAPAR)", "RdYlGn"),
    ("SLSTR_L1_RBT", "measurements/inadir/s8_bt_in", "SLSTR L1 RBT (S8 BT)", "inferno"),
    (
        "SLSTR_L2_LST",
        "measurements/lst",
        "SLSTR L2 LST (Land Surface Temp)",
        "RdYlBu_r",
    ),
]

for row, (product_key, subdataset, title, cmap) in enumerate(viz_products):
    url = PRODUCTS[product_key]["url"]
    path = eopfzarr_path(url, subdataset)

    try:
        ds = gdal.Open(path)
        geoloc = ds.GetMetadata("GEOLOCATION")

        if geoloc and ds:
            # Read coordinate arrays
            lon_ds = gdal.Open(geoloc["X_DATASET"])
            lat_ds = gdal.Open(geoloc["Y_DATASET"])

            # Read subset
            size = 300
            lon = lon_ds.GetRasterBand(1).ReadAsArray(
                0, 0, min(size, lon_ds.RasterXSize), min(size, lon_ds.RasterYSize)
            )
            lat = lat_ds.GetRasterBand(1).ReadAsArray(
                0, 0, min(size, lat_ds.RasterXSize), min(size, lat_ds.RasterYSize)
            )
            data = (
                ds.GetRasterBand(1)
                .ReadAsArray(0, 0, min(size, ds.RasterXSize), min(size, ds.RasterYSize))
                .astype(np.float32)
            )

            # Handle nodata
            nodata = ds.GetRasterBand(1).GetNoDataValue()
            if nodata is not None:
                data[data == nodata] = np.nan

            # Match dimensions
            min_r = min(data.shape[0], lat.shape[0], lon.shape[0])
            min_c = min(data.shape[1], lat.shape[1], lon.shape[1])
            data, lat, lon = (
                data[:min_r, :min_c],
                lat[:min_r, :min_c],
                lon[:min_r, :min_c],
            )

            vmin, vmax = np.nanpercentile(data, [2, 98])

            # Array view
            im1 = axes[row, 0].imshow(data, cmap=cmap, vmin=vmin, vmax=vmax)
            axes[row, 0].set_title(
                f"{title}\n(Array View)", fontweight="bold", fontsize=10
            )
            axes[row, 0].set_xlabel("Column")
            axes[row, 0].set_ylabel("Row")
            plt.colorbar(im1, ax=axes[row, 0], fraction=0.046)

            # Geographic view
            step = 2
            im2 = axes[row, 1].pcolormesh(
                lon[::step, ::step],
                lat[::step, ::step],
                data[::step, ::step],
                cmap=cmap,
                vmin=vmin,
                vmax=vmax,
                shading="auto",
            )
            axes[row, 1].set_title(
                f"{title}\n(Geographic View)", fontweight="bold", fontsize=10
            )
            axes[row, 1].set_xlabel("Longitude (°E)")
            axes[row, 1].set_ylabel("Latitude (°N)")
            axes[row, 1].grid(True, alpha=0.3)
            plt.colorbar(im2, ax=axes[row, 1], fraction=0.046)

            # Print coordinate range
            print(
                f"{product_key}: Lat [{np.nanmin(lat):.2f}° to {np.nanmax(lat):.2f}°], Lon [{np.nanmin(lon):.2f}° to {np.nanmax(lon):.2f}°]"
            )

            lon_ds = lat_ds = ds = None

    except Exception as e:
        print(f"{product_key}: Error - {str(e)[:50]}")
        axes[row, 0].text(
            0.5,
            0.5,
            f"Error: {str(e)[:30]}",
            ha="center",
            va="center",
            transform=axes[row, 0].transAxes,
        )
        axes[row, 1].text(
            0.5,
            0.5,
            f"Error: {str(e)[:30]}",
            ha="center",
            va="center",
            transform=axes[row, 1].transAxes,
        )

plt.suptitle(
    "EOPFZARR Geolocation Arrays - All 5 Sentinel-3 Products\nSwath (Array) → Geographic Transformation",
    fontsize=14,
    fontweight="bold",
    y=1.01,
)
plt.tight_layout()
plt.show()

# Summary table
print("\n" + "=" * 100)
print("SENTINEL-3 PRODUCT SUMMARY")
print("=" * 100)
print(f"{'Product':<20} {'Subdatasets':<15} {'Description'}")
print("-" * 100)

for key, subdatasets in product_subdatasets.items():
    count = len(subdatasets) if subdatasets else 0
    desc = PRODUCTS[key]["description"][:50]
    print(f"{key:<20} {count:<15} {desc}")