Tanager-1 Assets: Basic vs Ortho

# ==================
# 🧰 YOUR TOOLKIT
# ==================

# Essential imports for this lesson
import h5py
import numpy as np
import matplotlib.pyplot as plt


# ---- setup -------------------------------------------------------------

def download_tanager_data(url, file_path):
    import os
    import urllib.request

    def _progress(count, block_size, total_size):
        mb_done = count * block_size / 1_000_000
        mb_total = total_size / 1_000_000
        print(f"\rDownloading... {mb_done:.1f} / {mb_total:.1f} MB", end="", flush=True)

    if not os.path.exists(file_path):
        urllib.request.urlretrieve(url, file_path, reporthook=_progress)
        print(f"\nDownload complete: {file_path}")
    else:
        print(f"File already exists, skipping download: {file_path}")
    return file_path

# ---- io ----------------------------------------------------------------

def print_structure(name, obj):
    """Print one HDF5 object (group or dataset) plus its attributes.

    Designed as a callback for ``h5py.File.visititems`` to render the full
    hierarchy of a Tanager-1 file as an indented tree.

    Parameters
    ----------
    name : str
        Full HDF5 path of the object (provided by ``visititems``).
    obj : h5py.Group or h5py.Dataset
        The object itself (provided by ``visititems``).
    """
    level = name.count('/')
    indent = '    ' * level

    if isinstance(obj, h5py.Group):
        print(f"{indent}📂 {name}/")
    elif isinstance(obj, h5py.Dataset):
        print(f"{indent}📄 {name} | {obj.shape} | {obj.dtype}")

    if len(obj.attrs) > 0:
        for key, val in obj.attrs.items():
            val_str = str(val)
            if len(val_str) > 50:
                val_str = val_str[:50] + "..."
            print(f"{indent}    ↳ 🏷️ {key}: {val_str}")

def load_hdf5(file_path,
                      data_path='HDFEOS/SWATHS/HYP/Data Fields/surface_reflectance',
                      band=None):
    """Load a dataset (or one band) from a Tanager-1 HDF5 file, with its properties.

    Works for ANY dataset in the file — reflectance cubes, masks, geometry layers —
    because it never assumes what the values mean. Pass ``band`` to lazily read a
    single 2-D band straight from disk (cheap on RAM); leave it ``None`` for the
    full array.

    Parameters
    ----------
    file_path : str
        Path to the Tanager-1 ``.h5`` file.
    data_path : str, optional
        Internal HDF5 path to the dataset (Basic ``…/SWATHS/…`` by default; switch
        to ``…/GRIDS/…`` for Ortho, or to any mask / geometry layer).
    band : int, optional
        If given, read only that band index (lazy slice); if ``None``, read all.

    Returns
    -------
    data : ndarray
        The full array, or a single 2-D band if ``band`` is given.
    attrs : dict
        All attributes (properties) of the dataset — e.g. ``'wavelengths'``,
        ``'units'``, ``'good_wavelengths'``.
    """
    with h5py.File(file_path, 'r') as f:
        dset = f[data_path]
        data = dset[:] if band is None else dset[band]
        attrs = dict(dset.attrs)
    return data, attrs

# Download the Basic and Ortho SR assets for this scene using your Toolkit's
# download_tanager_data (defined in the Toolkit cell above).
basic_sr_url = "https://storage.googleapis.com/open-cogs/planet-stac/tanager1-release2-core-imagery/basic_sr_hdf5/20250510_112042_16_4001_basic_sr_hdf5.h5"
ortho_sr_url = "https://storage.googleapis.com/open-cogs/planet-stac/tanager1-release2-core-imagery/ortho_sr_hdf5/20250510_112042_16_4001_ortho_sr_hdf5.h5"

# These filenames are reused throughout the lesson:
# Note: You can change the file name to any name you want.
basic_sr_path = "20250510_112042_16_4001_basic_sr_hdf5.h5"
ortho_sr_path = "20250510_112042_16_4001_ortho_sr_hdf5.h5"

download_tanager_data(basic_sr_url, basic_sr_path)
download_tanager_data(ortho_sr_url, ortho_sr_path)

print("✅ Assets ready")

# Basic Surface Reflectance
with h5py.File(basic_sr_path, 'r') as f: 
    f.visititems(print_structure)

# Ortho Surface Reflectance
with h5py.File(ortho_sr_path, 'r') as f:
    f.visititems(print_structure)

# Define the dataset paths (Basic = SWATHS, Ortho = GRIDS)
basic_data_path = 'HDFEOS/SWATHS/HYP/Data Fields/surface_reflectance'
ortho_data_path = 'HDFEOS/GRIDS/HYP/Data Fields/surface_reflectance'

# Define the band
band_number = 100  # Feel free to change the band (0–425)

# Lazy single-band read via your Toolkit loader.
# Passing band=... reads just that 2-D slice straight from disk — no full-cube load.
basic_slice, basic_attrs = load_hdf5(basic_sr_path, basic_data_path, band=band_number)
ortho_slice, ortho_attrs = load_hdf5(ortho_sr_path, ortho_data_path, band=band_number)

# Full cube dimensions without loading the cube: the band count is simply the
# number of wavelengths reported in the dataset's attributes.
basic_bands = len(basic_attrs['wavelengths'])
ortho_bands = len(ortho_attrs['wavelengths'])
print(f"Basic Data Shape (Bands, Rows, Cols): ({basic_bands}, {basic_slice.shape[0]}, {basic_slice.shape[1]})")
print(f"Loaded Basic Slice: {basic_slice.shape}")
print(f"Ortho Data Shape (Bands, Rows, Cols): ({ortho_bands}, {ortho_slice.shape[0]}, {ortho_slice.shape[1]})")
print(f"Loaded Ortho Slice: {ortho_slice.shape}")

# Create side-by-side figure: unpack axes as (basic, ortho) for clarity
fig, (basic, ortho) = plt.subplots(1, 2, figsize=(14, 8))

# Left: Basic (sensor geometry)
basic.imshow(basic_slice, cmap='gray')
basic.set_title("Basic (Sensor Geometry)\nRaw Perspective")
basic.set_xlabel("Sample (Column)")
basic.set_ylabel("Line (Row)")

# Right: Ortho (map geometry)
ortho.imshow(ortho_slice, cmap='gray')
ortho.set_title("Ortho (Map Geometry)\nProjected North-Up")
ortho.set_xlabel("Sample (Column)")
ortho.set_ylabel("Line (Row)")

plt.tight_layout()
plt.show()

# plot histogram for Ortho Slice

plt.figure(figsize=(10, 5))
plt.hist(ortho_slice.flatten(), bins=100, edgecolor='black', alpha=0.7)
plt.xlabel('Surface Reflectance')
plt.ylabel('Pixel Count')
plt.title(f'Ortho Band {band_number} Distribution')
plt.grid(True, alpha=0.3)
plt.show()

# replace -9999 with np.nan
ortho_slice[ortho_slice == -9999] = np.nan

# Replot the histogram after removing `-9999` values
plt.figure(figsize=(10, 5))
plt.hist(ortho_slice.flatten(), bins=100, edgecolor='black', alpha=0.7)
plt.xlabel('Surface Reflectance')
plt.ylabel('Pixel Count')
plt.title(f'Ortho Band {band_number} Distribution')
plt.grid(True, alpha=0.3)
plt.show()

# Create side-by-side figure: unpack axes as (Basic, Ortho) for clarity
fig, (basic, ortho) = plt.subplots(1, 2, figsize=(14, 8))

# Left: Basic (sensor geometry)
basic.imshow(basic_slice, cmap='gray')
basic.set_title("Basic (Sensor Geometry)\nRaw Perspective")
basic.set_xlabel("Sample (Column)")
basic.set_ylabel("Line (Row)")

# Right: Ortho (map geometry) 
ortho.imshow(ortho_slice, cmap='gray')
ortho.set_title("Ortho (Map Geometry)\nProjected North-Up")
ortho.set_xlabel("Sample (Column)")
ortho.set_ylabel("Line (Row)")

plt.tight_layout()
plt.show()