Harmonizing BIOMASS Beta0, Gamma0, and LUT georeferencing¶

This notebook downloads one BIOMASS L1B granule to a local cache, then follows its radar-coordinate arrays to a UTM Gamma0 chip. It deliberately plots each coordinate system rather than making the arrays merely look the same.

gammaNought is not an image whose corners match Beta0. It is a coarse lookup table on relative azimuth time and slant-range time, with padding outside the Beta0 acquisition. The LUT also provides geometry/longitude and geometry/latitude on those same physical axes. For each Beta0 pixel, first derive its physical radar coordinates, locate that coordinate in the LUT vectors, and bilinearly sample the radiometry and geometry arrays. Calculate Gamma0 in radar geometry, then warp it once to the final map grid using the sampled longitude/latitude arrays.

This notebook demonstrates a few things:

Beta0 and the LUT have different array shapes and different physical extents.
The LUT's longitude/latitude arrays supply denser georeferencing than the Beta0 TIFF's GCPs.
The safe alignment maps every Beta0 row and column through the annotation and LUT coordinate vectors before bilinear sampling.

Before running the first code cell, set ESA_MAAP_CLIENT_SECRET and ESA_OFFLINE_TOKEN in the environment. It reuses cached assets when available, otherwise retrieves the item from the ESA STAC API and downloads the required TIFF, NetCDF, and XML files. Set BIOMASS_CACHE_DIR to use a cache location other than /tmp/biomass-cache.

Install the dependencies with the following command:

pip install affine matplotlib netCDF4 numpy obstore pyproj pystac-client rasterio requests scipy

In [1]:

Copied!





import asyncio
import os
from pathlib import Path
from tempfile import NamedTemporaryFile
from urllib.parse import urlparse
from xml.etree import ElementTree

import matplotlib.pyplot as plt
import numpy as np
import obstore as obs
import rasterio
import requests
from affine import Affine
from netCDF4 import Dataset
from obstore.store import HTTPStore
from pyproj import CRS, Transformer
from pystac_client import Client
from rasterio.enums import Resampling
from rasterio.warp import reproject
from scipy.ndimage import map_coordinates

plt.rcParams.update(
    {"figure.dpi": 120, "axes.spines.top": False, "axes.spines.right": False}
)

ESA_STAC_API_URL = "https://catalog.maap.eo.esa.int/catalogue/"
CACHE_DIR = Path(os.environ.get("BIOMASS_CACHE_DIR", "/tmp/biomass-cache"))
NODATA = -9999.0
DISPLAY_MAX_DIMENSION = 900
UTM_RESOLUTION_M = 25.0


async def fetch_assets(urls: dict[str, str], token: str) -> dict[str, bytes]:
    """Download the product's authenticated assets concurrently through obstore."""
    downloads = {}
    for name, url in urls.items():
        parsed_url = urlparse(url)
        store = HTTPStore(
            f"{parsed_url.scheme}://{parsed_url.netloc}",
            client_options={
                "default_headers": {"Authorization": f"Bearer {token}"},
                "timeout": "3m",
            },
        )
        print(f"Starting download of {url}")
        downloads[name] = obs.get_async(store, parsed_url.path.lstrip("/"))

    responses = await asyncio.gather(*downloads.values())
    contents = await asyncio.gather(*(response.bytes_async() for response in responses))
    print(f"Completed {len(contents)} asset downloads")
    return dict(zip(downloads, map(bytes, contents), strict=True))


def cache_paths(item_id: str, cache_dir: Path) -> dict[str, Path]:
    """Return stable paths for an item's source assets."""
    safe_item_id = item_id.replace("/", "_")
    prefix = cache_dir / f"biomass__{safe_item_id}"
    return {
        "beta": prefix.with_name(f"{prefix.name}__beta.tif"),
        "lut": prefix.with_name(f"{prefix.name}__lut.nc"),
        "annotation": prefix.with_name(f"{prefix.name}__annotation.xml"),
    }


def write_cached_asset(path: Path, contents: bytes) -> None:
    """Atomically write one downloaded asset to the local cache."""
    path.parent.mkdir(parents=True, exist_ok=True)
    with NamedTemporaryFile(dir=path.parent, delete=False) as temporary:
        temporary.write(contents)
        temporary_path = Path(temporary.name)

    temporary_path.replace(path)
import asyncio
import os
from pathlib import Path
from tempfile import NamedTemporaryFile
from urllib.parse import urlparse
from xml.etree import ElementTree

import matplotlib.pyplot as plt
import numpy as np
import obstore as obs
import rasterio
import requests
from affine import Affine
from netCDF4 import Dataset
from obstore.store import HTTPStore
from pyproj import CRS, Transformer
from pystac_client import Client
from rasterio.enums import Resampling
from rasterio.warp import reproject
from scipy.ndimage import map_coordinates

plt.rcParams.update(
    {"figure.dpi": 120, "axes.spines.top": False, "axes.spines.right": False}
)

ESA_STAC_API_URL = "https://catalog.maap.eo.esa.int/catalogue/"
CACHE_DIR = Path(os.environ.get("BIOMASS_CACHE_DIR", "/tmp/biomass-cache"))
NODATA = -9999.0
DISPLAY_MAX_DIMENSION = 900
UTM_RESOLUTION_M = 25.0


async def fetch_assets(urls: dict[str, str], token: str) -> dict[str, bytes]:
    """Download the product's authenticated assets concurrently through obstore."""
    downloads = {}
    for name, url in urls.items():
        parsed_url = urlparse(url)
        store = HTTPStore(
            f"{parsed_url.scheme}://{parsed_url.netloc}",
            client_options={
                "default_headers": {"Authorization": f"Bearer {token}"},
                "timeout": "3m",
            },
        )
        print(f"Starting download of {url}")
        downloads[name] = obs.get_async(store, parsed_url.path.lstrip("/"))

    responses = await asyncio.gather(*downloads.values())
    contents = await asyncio.gather(*(response.bytes_async() for response in responses))
    print(f"Completed {len(contents)} asset downloads")
    return dict(zip(downloads, map(bytes, contents), strict=True))


def cache_paths(item_id: str, cache_dir: Path) -> dict[str, Path]:
    """Return stable paths for an item's source assets."""
    safe_item_id = item_id.replace("/", "_")
    prefix = cache_dir / f"biomass__{safe_item_id}"
    return {
        "beta": prefix.with_name(f"{prefix.name}__beta.tif"),
        "lut": prefix.with_name(f"{prefix.name}__lut.nc"),
        "annotation": prefix.with_name(f"{prefix.name}__annotation.xml"),
    }


def write_cached_asset(path: Path, contents: bytes) -> None:
    """Atomically write one downloaded asset to the local cache."""
    path.parent.mkdir(parents=True, exist_ok=True)
    with NamedTemporaryFile(dir=path.parent, delete=False) as temporary:
        temporary.write(contents)
        temporary_path = Path(temporary.name)

    temporary_path.replace(path)

In [2]:

Copied!





ITEM_ID = (
    "BIO_S2_DGM__1S_20260601T121355_20260601T121415_T_G01_M03_C07_T008_F100_02_DSD6MY"
)

item = next(
    Client.open(ESA_STAC_API_URL)
    .search(ids=[ITEM_ID], collections=["BiomassLevel1b"], limit=1)
    .items()
)
if not item:
    raise ValueError(f"No BIOMASS L1B item found with id {ITEM_ID!r}")

try:
    asset_urls = {
        "beta": item.assets["enclosure_tiff"].href,
        "lut": item.assets["enclosure_nc"].href,
        "annotation": item.assets["enclosure_annot_xml"].href,
    }
except KeyError as error:
    raise ValueError(f"{item.id} is missing required asset {error.args[0]}") from error

paths = cache_paths(item.id, CACHE_DIR)
missing_urls = {
    name: url for name, url in asset_urls.items() if not paths[name].exists()
}
if missing_urls:
    client_secret = os.getenv("ESA_MAAP_CLIENT_SECRET")
    offline_token = os.getenv("ESA_OFFLINE_TOKEN")
    if not all((client_secret, offline_token)):
        raise ValueError(
            "Set ESA_MAAP_CLIENT_SECRET and ESA_OFFLINE_TOKEN before downloading source assets."
        )

    response = requests.post(
        "https://iam.maap.eo.esa.int/realms/esa-maap/protocol/openid-connect/token",
        data={
            "client_id": os.getenv("MAAP_CLIENT_ID", "offline-token"),
            "client_secret": client_secret,
            "grant_type": "refresh_token",
            "refresh_token": offline_token,
            "scope": "offline_access openid",
        },
        timeout=30,
    )
    response.raise_for_status()
    access_token = response.json().get("access_token")
    if not access_token:
        raise RuntimeError("The IAM response did not include an access token.")
    for name, contents in (await fetch_assets(missing_urls, access_token)).items():
        write_cached_asset(paths[name], contents)

assert all(path.exists() for path in paths.values())
paths
ITEM_ID = (
    "BIO_S2_DGM__1S_20260601T121355_20260601T121415_T_G01_M03_C07_T008_F100_02_DSD6MY"
)

item = next(
    Client.open(ESA_STAC_API_URL)
    .search(ids=[ITEM_ID], collections=["BiomassLevel1b"], limit=1)
    .items()
)
if not item:
    raise ValueError(f"No BIOMASS L1B item found with id {ITEM_ID!r}")

try:
    asset_urls = {
        "beta": item.assets["enclosure_tiff"].href,
        "lut": item.assets["enclosure_nc"].href,
        "annotation": item.assets["enclosure_annot_xml"].href,
    }
except KeyError as error:
    raise ValueError(f"{item.id} is missing required asset {error.args[0]}") from error

paths = cache_paths(item.id, CACHE_DIR)
missing_urls = {
    name: url for name, url in asset_urls.items() if not paths[name].exists()
}
if missing_urls:
    client_secret = os.getenv("ESA_MAAP_CLIENT_SECRET")
    offline_token = os.getenv("ESA_OFFLINE_TOKEN")
    if not all((client_secret, offline_token)):
        raise ValueError(
            "Set ESA_MAAP_CLIENT_SECRET and ESA_OFFLINE_TOKEN before downloading source assets."
        )

    response = requests.post(
        "https://iam.maap.eo.esa.int/realms/esa-maap/protocol/openid-connect/token",
        data={
            "client_id": os.getenv("MAAP_CLIENT_ID", "offline-token"),
            "client_secret": client_secret,
            "grant_type": "refresh_token",
            "refresh_token": offline_token,
            "scope": "offline_access openid",
        },
        timeout=30,
    )
    response.raise_for_status()
    access_token = response.json().get("access_token")
    if not access_token:
        raise RuntimeError("The IAM response did not include an access token.")
    for name, contents in (await fetch_assets(missing_urls, access_token)).items():
        write_cached_asset(paths[name], contents)

assert all(path.exists() for path in paths.values())
paths

Out[2]:

{'beta': PosixPath('/tmp/biomass-cache/biomass__BIO_S2_DGM__1S_20260601T121355_20260601T121415_T_G01_M03_C07_T008_F100_02_DSD6MY__beta.tif'),
 'lut': PosixPath('/tmp/biomass-cache/biomass__BIO_S2_DGM__1S_20260601T121355_20260601T121415_T_G01_M03_C07_T008_F100_02_DSD6MY__lut.nc'),
 'annotation': PosixPath('/tmp/biomass-cache/biomass__BIO_S2_DGM__1S_20260601T121355_20260601T121415_T_G01_M03_C07_T008_F100_02_DSD6MY__annotation.xml')}

1. Load the native arrays and their metadata¶

The TIFF's four bands are Beta0 amplitudes in radar pixels (line, sample). Its GCPs link selected radar pixels to longitude/latitude; the first plot overlays them for reference. The NetCDF LUT has its own coarse shape and coordinate vectors. Its first dimension is azimuth time and its second is slant-range time, matching the Beta0 line/sample order. gammaNought, geometry/longitude, and geometry/latitude must share that LUT grid.

In [3]:

Copied!





with rasterio.open(paths["beta"]) as source:
    beta = source.read(masked=True).filled(np.nan).astype("float32")
    beta_height, beta_width = source.height, source.width
    beta_nodata = source.nodata
    gcps, _ = source.gcps

gcp_rows = np.asarray([gcp.row for gcp in gcps])
gcp_cols = np.asarray([gcp.col for gcp in gcps])

with Dataset(paths["lut"]) as dataset:
    gamma_lut = np.asarray(dataset["radiometry/gammaNought"][:], dtype="float32")
    longitude_lut = np.asarray(dataset["geometry/longitude"][:], dtype="float64")
    latitude_lut = np.asarray(dataset["geometry/latitude"][:], dtype="float64")
    lut_azimuth = np.asarray(dataset["relativeAzimuthTimeRGC"][:], dtype="float64")
    lut_range = np.asarray(dataset["slantRangeTimeRGC"][:], dtype="float64")

for name, values in {
    "gammaNought": gamma_lut,
    "geometry/longitude": longitude_lut,
    "geometry/latitude": latitude_lut,
}.items():
    if values.shape != gamma_lut.shape:
        raise ValueError(f"{name} does not match the gammaNought LUT shape")

gamma_lut = np.where(gamma_lut == NODATA, np.nan, gamma_lut)
longitude_lut = np.where(longitude_lut == NODATA, np.nan, longitude_lut)
latitude_lut = np.where(latitude_lut == NODATA, np.nan, latitude_lut)

sar_image = ElementTree.parse(paths["annotation"]).getroot().find("sarImage")
if sar_image is None:
    raise ValueError("The product annotation has no sarImage element.")

azimuth_step = float(sar_image.findtext("azimuthTimeInterval", ""))
range_spacing = float(sar_image.findtext("rangePixelSpacing", ""))

ground_to_slant = np.fromstring(
    sar_image.findtext(
        "rangeCoordinateConversion/coordinateConversion/groundToSlantCoefficients", ""
    ),
    sep=" ",
)

if not ground_to_slant.size:
    raise ValueError("The product annotation has no ground-to-slant coefficients.")

print(f"Beta0: {beta.shape} (band, line, sample); {len(gcps)} GCPs")
print(f"LUT:   {gamma_lut.shape} (relativeAzimuthTimeRGC, slantRangeTimeRGC)")
print(f"Beta0 nodata: {beta_nodata}")
with rasterio.open(paths["beta"]) as source:
    beta = source.read(masked=True).filled(np.nan).astype("float32")
    beta_height, beta_width = source.height, source.width
    beta_nodata = source.nodata
    gcps, _ = source.gcps

gcp_rows = np.asarray([gcp.row for gcp in gcps])
gcp_cols = np.asarray([gcp.col for gcp in gcps])

with Dataset(paths["lut"]) as dataset:
    gamma_lut = np.asarray(dataset["radiometry/gammaNought"][:], dtype="float32")
    longitude_lut = np.asarray(dataset["geometry/longitude"][:], dtype="float64")
    latitude_lut = np.asarray(dataset["geometry/latitude"][:], dtype="float64")
    lut_azimuth = np.asarray(dataset["relativeAzimuthTimeRGC"][:], dtype="float64")
    lut_range = np.asarray(dataset["slantRangeTimeRGC"][:], dtype="float64")

for name, values in {
    "gammaNought": gamma_lut,
    "geometry/longitude": longitude_lut,
    "geometry/latitude": latitude_lut,
}.items():
    if values.shape != gamma_lut.shape:
        raise ValueError(f"{name} does not match the gammaNought LUT shape")

gamma_lut = np.where(gamma_lut == NODATA, np.nan, gamma_lut)
longitude_lut = np.where(longitude_lut == NODATA, np.nan, longitude_lut)
latitude_lut = np.where(latitude_lut == NODATA, np.nan, latitude_lut)

sar_image = ElementTree.parse(paths["annotation"]).getroot().find("sarImage")
if sar_image is None:
    raise ValueError("The product annotation has no sarImage element.")

azimuth_step = float(sar_image.findtext("azimuthTimeInterval", ""))
range_spacing = float(sar_image.findtext("rangePixelSpacing", ""))

ground_to_slant = np.fromstring(
    sar_image.findtext(
        "rangeCoordinateConversion/coordinateConversion/groundToSlantCoefficients", ""
    ),
    sep=" ",
)

if not ground_to_slant.size:
    raise ValueError("The product annotation has no ground-to-slant coefficients.")

print(f"Beta0: {beta.shape} (band, line, sample); {len(gcps)} GCPs")
print(f"LUT:   {gamma_lut.shape} (relativeAzimuthTimeRGC, slantRangeTimeRGC)")
print(f"Beta0 nodata: {beta_nodata}")

Beta0: (4, 6191, 2627) (band, line, sample); 35 GCPs
LUT:   (2120, 806) (relativeAzimuthTimeRGC, slantRangeTimeRGC)
Beta0 nodata: -9999.0

In [8]:

Copied!





def display_stride(shape: tuple[int, int], maximum: int = DISPLAY_MAX_DIMENSION) -> int:
    """Return a stride that limits a two-dimensional display's longest side."""
    return max(1, int(np.ceil(max(shape) / maximum)))


def finite_display_limits(array: np.ndarray) -> tuple[float, float]:
    """Return 2nd-to-98th percentile limits from finite values for plotting."""
    values = array[np.isfinite(array)]
    if not values.size:
        raise ValueError("Cannot plot an array with no finite values.")
    vmin, vmax = np.percentile(values, (2, 98))
    if vmin == vmax:
        padding = max(abs(vmin) * 0.01, 1.0)
        return float(vmin - padding), float(vmax + padding)
    return float(vmin), float(vmax)


def plot_radar_array(
    array: np.ndarray,
    title: str,
    *,
    cmap: str = "viridis",
    points: tuple[np.ndarray, np.ndarray] | None = None,
) -> None:
    """Plot a radar-grid array in line/sample coordinates."""
    height, width = array.shape
    stride = display_stride(array.shape)
    vmin, vmax = finite_display_limits(array)
    figure, axis = plt.subplots(figsize=(10, 5))
    image = axis.imshow(
        array[::stride, ::stride],
        extent=(0, width, height, 0),
        cmap=cmap,
        interpolation="nearest",
        vmin=vmin,
        vmax=vmax,
    )
    if points is not None:
        axis.scatter(
            points[1],
            points[0],
            s=22,
            facecolors="none",
            edgecolors="crimson",
            label="GCP",
        )
        axis.legend(loc="lower right")
    axis.set(title=title, xlabel="column / sample", ylabel="row / line")
    figure.colorbar(image, ax=axis, label="array value")
    plt.show()


gcp_rows = np.asarray([gcp.row for gcp in gcps])
gcp_cols = np.asarray([gcp.col for gcp in gcps])

plot_radar_array(
    beta[0], "Beta0 HH amplitude in the native radar grid", points=(gcp_rows, gcp_cols)
)
def display_stride(shape: tuple[int, int], maximum: int = DISPLAY_MAX_DIMENSION) -> int:
    """Return a stride that limits a two-dimensional display's longest side."""
    return max(1, int(np.ceil(max(shape) / maximum)))


def finite_display_limits(array: np.ndarray) -> tuple[float, float]:
    """Return 2nd-to-98th percentile limits from finite values for plotting."""
    values = array[np.isfinite(array)]
    if not values.size:
        raise ValueError("Cannot plot an array with no finite values.")
    vmin, vmax = np.percentile(values, (2, 98))
    if vmin == vmax:
        padding = max(abs(vmin) * 0.01, 1.0)
        return float(vmin - padding), float(vmax + padding)
    return float(vmin), float(vmax)


def plot_radar_array(
    array: np.ndarray,
    title: str,
    *,
    cmap: str = "viridis",
    points: tuple[np.ndarray, np.ndarray] | None = None,
) -> None:
    """Plot a radar-grid array in line/sample coordinates."""
    height, width = array.shape
    stride = display_stride(array.shape)
    vmin, vmax = finite_display_limits(array)
    figure, axis = plt.subplots(figsize=(10, 5))
    image = axis.imshow(
        array[::stride, ::stride],
        extent=(0, width, height, 0),
        cmap=cmap,
        interpolation="nearest",
        vmin=vmin,
        vmax=vmax,
    )
    if points is not None:
        axis.scatter(
            points[1],
            points[0],
            s=22,
            facecolors="none",
            edgecolors="crimson",
            label="GCP",
        )
        axis.legend(loc="lower right")
    axis.set(title=title, xlabel="column / sample", ylabel="row / line")
    figure.colorbar(image, ax=axis, label="array value")
    plt.show()


gcp_rows = np.asarray([gcp.row for gcp in gcps])
gcp_cols = np.asarray([gcp.col for gcp in gcps])

plot_radar_array(
    beta[0], "Beta0 HH amplitude in the native radar grid", points=(gcp_rows, gcp_cols)
)

No description has been provided for this image

2. Put both arrays on their physical radar axes¶

The annotation converts a Beta0 line to relative azimuth time. It converts a Beta0 sample to ground range and then to slant-range time through the product polynomial. These arrays, not shape ratios, tell us where each Beta0 pixel belongs in the LUT.

The table below shows the important mismatch: the LUT extends before and after the Beta0 interval in both directions. That padding is why mapping Beta0's corners to the LUT's corners is unsafe.

In [9]:

Copied!





beta_azimuth = np.arange(beta_height, dtype="float64") * azimuth_step
beta_range = np.polynomial.polynomial.polyval(
    np.arange(beta_width, dtype="float64") * range_spacing,
    ground_to_slant,
)

coverage = {
    "Beta0 azimuth time": (beta_azimuth[0], beta_azimuth[-1]),
    "LUT azimuth time": (lut_azimuth[0], lut_azimuth[-1]),
    "Beta0 slant-range time": (beta_range[0], beta_range[-1]),
    "LUT slant-range time": (lut_range[0], lut_range[-1]),
}
for name, (start, stop) in coverage.items():
    print(f"{name:24}: {start:.12g} to {stop:.12g}")

assert lut_azimuth[0] < beta_azimuth[0] < beta_azimuth[-1] < lut_azimuth[-1]
assert lut_range[0] < beta_range[0] < beta_range[-1] < lut_range[-1]
beta_azimuth = np.arange(beta_height, dtype="float64") * azimuth_step
beta_range = np.polynomial.polynomial.polyval(
    np.arange(beta_width, dtype="float64") * range_spacing,
    ground_to_slant,
)

coverage = {
    "Beta0 azimuth time": (beta_azimuth[0], beta_azimuth[-1]),
    "LUT azimuth time": (lut_azimuth[0], lut_azimuth[-1]),
    "Beta0 slant-range time": (beta_range[0], beta_range[-1]),
    "LUT slant-range time": (lut_range[0], lut_range[-1]),
}
for name, (start, stop) in coverage.items():
    print(f"{name:24}: {start:.12g} to {stop:.12g}")

assert lut_azimuth[0] < beta_azimuth[0] < beta_azimuth[-1] < lut_azimuth[-1]
assert lut_range[0] < beta_range[0] < beta_range[-1] < lut_range[-1]

Beta0 azimuth time      : 0 to 20.2525849829
LUT azimuth time        : -2.35039921771 to 22.6083761108
Beta0 slant-range time  : 0.0049582390078 to 0.00516271734917
LUT slant-range time    : 0.0049540091194 to 0.00516682521111

In [10]:

Copied!





def physical_radar_plot(
    values: np.ndarray, x: np.ndarray, y: np.ndarray, title: str
) -> None:
    """Plot a downsampled radar array against its physical coordinate vectors."""
    stride = display_stride(values.shape)
    vmin, vmax = finite_display_limits(values)
    figure, axis = plt.subplots(figsize=(9, 5))
    mesh = axis.pcolormesh(
        x[::stride],
        y[::stride],
        values[::stride, ::stride],
        shading="auto",
        cmap="viridis",
        vmin=vmin,
        vmax=vmax,
    )
    axis.set(
        title=title, xlabel="slant-range time (s)", ylabel="relative azimuth time (s)"
    )
    figure.colorbar(mesh, ax=axis, label="array value")
    plt.show()


physical_radar_plot(
    beta[0], beta_range, beta_azimuth, "Beta0 HH on its physical radar coordinates"
)
physical_radar_plot(
    gamma_lut,
    lut_range,
    lut_azimuth,
    "Native gammaNought LUT on its physical radar coordinates",
)
def physical_radar_plot(
    values: np.ndarray, x: np.ndarray, y: np.ndarray, title: str
) -> None:
    """Plot a downsampled radar array against its physical coordinate vectors."""
    stride = display_stride(values.shape)
    vmin, vmax = finite_display_limits(values)
    figure, axis = plt.subplots(figsize=(9, 5))
    mesh = axis.pcolormesh(
        x[::stride],
        y[::stride],
        values[::stride, ::stride],
        shading="auto",
        cmap="viridis",
        vmin=vmin,
        vmax=vmax,
    )
    axis.set(
        title=title, xlabel="slant-range time (s)", ylabel="relative azimuth time (s)"
    )
    figure.colorbar(mesh, ax=axis, label="array value")
    plt.show()


physical_radar_plot(
    beta[0], beta_range, beta_azimuth, "Beta0 HH on its physical radar coordinates"
)
physical_radar_plot(
    gamma_lut,
    lut_range,
    lut_azimuth,
    "Native gammaNought LUT on its physical radar coordinates",
)

3. Sample radiometry and georeferencing from the native LUT grid¶

The LUT maps a Beta0 line/sample to fractional row/column coordinates through the same physical radar axes. Sample all three LUT arrays at those coordinates: gammaNought corrects radiometry, while longitude and latitude are the source geolocation for the final warp. No TIFF GCPs are used for that warp.

In [11]:

Copied!





lut_row_for_beta_row = np.interp(
    beta_azimuth, lut_azimuth, np.arange(lut_azimuth.size, dtype="float64")
)
lut_col_for_beta_col = np.interp(
    beta_range, lut_range, np.arange(lut_range.size, dtype="float64")
)
lut_rows, lut_cols = np.broadcast_arrays(
    lut_row_for_beta_row[:, None], lut_col_for_beta_col[None, :]
)

assert np.all(np.diff(lut_row_for_beta_row) > 0)
assert np.all(np.diff(lut_col_for_beta_col) > 0)


def sample_lut(values: np.ndarray) -> np.ndarray:
    """Bilinearly sample a LUT array onto the native Beta0 radar grid."""
    if values.shape != gamma_lut.shape:
        raise ValueError("LUT values do not match the gammaNought LUT shape")
    return map_coordinates(
        values,
        np.stack((lut_rows, lut_cols)),
        order=1,
        mode="nearest",
        prefilter=False,
    ).astype("float32")


longitude = sample_lut(longitude_lut)
latitude = sample_lut(latitude_lut)
if not np.isfinite(longitude).any() or not np.isfinite(latitude).any():
    raise ValueError("The sampled geometry LUT has no valid longitude/latitude pairs.")

plot_radar_array(longitude, "Longitude sampled from the LUT onto the Beta0 grid")
plot_radar_array(latitude, "Latitude sampled from the LUT onto the Beta0 grid")
lut_row_for_beta_row = np.interp(
    beta_azimuth, lut_azimuth, np.arange(lut_azimuth.size, dtype="float64")
)
lut_col_for_beta_col = np.interp(
    beta_range, lut_range, np.arange(lut_range.size, dtype="float64")
)
lut_rows, lut_cols = np.broadcast_arrays(
    lut_row_for_beta_row[:, None], lut_col_for_beta_col[None, :]
)

assert np.all(np.diff(lut_row_for_beta_row) > 0)
assert np.all(np.diff(lut_col_for_beta_col) > 0)


def sample_lut(values: np.ndarray) -> np.ndarray:
    """Bilinearly sample a LUT array onto the native Beta0 radar grid."""
    if values.shape != gamma_lut.shape:
        raise ValueError("LUT values do not match the gammaNought LUT shape")
    return map_coordinates(
        values,
        np.stack((lut_rows, lut_cols)),
        order=1,
        mode="nearest",
        prefilter=False,
    ).astype("float32")


longitude = sample_lut(longitude_lut)
latitude = sample_lut(latitude_lut)
if not np.isfinite(longitude).any() or not np.isfinite(latitude).any():
    raise ValueError("The sampled geometry LUT has no valid longitude/latitude pairs.")

plot_radar_array(longitude, "Longitude sampled from the LUT onto the Beta0 grid")
plot_radar_array(latitude, "Latitude sampled from the LUT onto the Beta0 grid")

Beta0 contains amplitude, so its linear intensity is amplitude². The LUT is already an intensity-domain correction factor. Nodata becomes NaN before the calculation so it cannot leak into valid output values.

At this point beta, gamma_nought, gamma0, longitude, and latitude have the same native radar grid. This is the last safe point for the radiometric calculation. The three panels below separate the Beta0 intensity, sampled correction factor, and corrected Gamma0. Every plot uses the finite 2nd-to-98th percentile range for display contrast only; it does not clip the calculation or output arrays.

In [12]:

Copied!





gamma_nought = sample_lut(gamma_lut)
beta = np.where(beta == beta_nodata, np.nan, beta)
gamma0 = beta**2 * gamma_nought[None, :, :]

assert gamma0.shape == beta.shape
assert np.isnan(gamma0[np.isnan(beta)]).all()

stride = display_stride((beta_height, beta_width))
figure, axes = plt.subplots(1, 3, figsize=(18, 5), constrained_layout=True)
for axis, values, title, cmap in zip(
    axes,
    (beta[0] ** 2, gamma_nought, gamma0[0]),
    (
        "Beta0 HH intensity in the native radar grid",
        "Sampled gammaNought correction factor",
        "Gamma0 HH in the native Beta0 radar grid",
    ),
    ("magma", "viridis", "magma"),
    strict=True,
):
    vmin, vmax = finite_display_limits(values)
    image = axis.imshow(
        values[::stride, ::stride],
        extent=(0, beta_width, beta_height, 0),
        cmap=cmap,
        interpolation="nearest",
        vmin=vmin,
        vmax=vmax,
    )
    axis.set(title=title, xlabel="column / sample", ylabel="row / line")
    figure.colorbar(image, ax=axis, label="array value")
plt.show()
gamma_nought = sample_lut(gamma_lut)
beta = np.where(beta == beta_nodata, np.nan, beta)
gamma0 = beta**2 * gamma_nought[None, :, :]

assert gamma0.shape == beta.shape
assert np.isnan(gamma0[np.isnan(beta)]).all()

stride = display_stride((beta_height, beta_width))
figure, axes = plt.subplots(1, 3, figsize=(18, 5), constrained_layout=True)
for axis, values, title, cmap in zip(
    axes,
    (beta[0] ** 2, gamma_nought, gamma0[0]),
    (
        "Beta0 HH intensity in the native radar grid",
        "Sampled gammaNought correction factor",
        "Gamma0 HH in the native Beta0 radar grid",
    ),
    ("magma", "viridis", "magma"),
    strict=True,
):
    vmin, vmax = finite_display_limits(values)
    image = axis.imshow(
        values[::stride, ::stride],
        extent=(0, beta_width, beta_height, 0),
        cmap=cmap,
        interpolation="nearest",
        vmin=vmin,
        vmax=vmax,
    )
    axis.set(title=title, xlabel="column / sample", ylabel="row / line")
    figure.colorbar(image, ax=axis, label="array value")
plt.show()

5. Warp calibrated Gamma0 to UTM with LUT longitude/latitude¶

For an illustrative chip, the next cell builds a 100 m north-up UTM grid from every valid LUT geolocation pixel in the chip. Production processing should instead supply its already-defined target grid, for example exact MGRS tile bounds at 100 m. The target_transform, target_height, and target_width must be identical for every acquisition intended for a time series. Do not let each acquisition choose its own final grid.

reproject receives src_geoloc_array=(longitude, latitude) and not the TIFF's GCPs or an approximate src_transform. It uses the sampled LUT geography directly, keeps NaN nodata throughout, and is the single map-space interpolation.

In [13]:

Copied!





chip_size = 1_024
chip_row_start = (beta_height - chip_size) // 2
chip_col_start = (beta_width - chip_size) // 2
chip_slice = np.s_[
    chip_row_start : chip_row_start + chip_size,
    chip_col_start : chip_col_start + chip_size,
]
chip_gamma0 = gamma0[0, chip_slice[0], chip_slice[1]]
chip_longitude = longitude[chip_slice]
chip_latitude = latitude[chip_slice]

valid_geolocation = np.isfinite(chip_longitude) & np.isfinite(chip_latitude)
if not valid_geolocation.any():
    raise ValueError("The selected chip has no valid LUT geolocation.")

centre_lon = float(np.mean(chip_longitude[valid_geolocation]))
centre_lat = float(np.mean(chip_latitude[valid_geolocation]))
zone = int(np.floor((centre_lon + 180) / 6) + 1)
utm_epsg = (32600 if centre_lat >= 0 else 32700) + zone
utm_crs = CRS.from_epsg(utm_epsg)
to_utm = Transformer.from_crs("EPSG:4326", utm_crs, always_xy=True)
chip_x, chip_y = to_utm.transform(
    chip_longitude[valid_geolocation], chip_latitude[valid_geolocation]
)

xmin = np.floor(np.min(chip_x) / UTM_RESOLUTION_M) * UTM_RESOLUTION_M
ymin = np.floor(np.min(chip_y) / UTM_RESOLUTION_M) * UTM_RESOLUTION_M
xmax = np.ceil(np.max(chip_x) / UTM_RESOLUTION_M) * UTM_RESOLUTION_M
ymax = np.ceil(np.max(chip_y) / UTM_RESOLUTION_M) * UTM_RESOLUTION_M
target_width = int(round((xmax - xmin) / UTM_RESOLUTION_M))
target_height = int(round((ymax - ymin) / UTM_RESOLUTION_M))
target_transform = Affine.translation(xmin, ymax) * Affine.scale(
    UTM_RESOLUTION_M, -UTM_RESOLUTION_M
)

print(
    f"UTM target: EPSG:{utm_epsg}, {target_width} × {target_height} pixels at {UTM_RESOLUTION_M:g} m"
)
assert target_width > 0 and target_height > 0
assert (
    target_transform.a == UTM_RESOLUTION_M and target_transform.e == -UTM_RESOLUTION_M
)
chip_size = 1_024
chip_row_start = (beta_height - chip_size) // 2
chip_col_start = (beta_width - chip_size) // 2
chip_slice = np.s_[
    chip_row_start : chip_row_start + chip_size,
    chip_col_start : chip_col_start + chip_size,
]
chip_gamma0 = gamma0[0, chip_slice[0], chip_slice[1]]
chip_longitude = longitude[chip_slice]
chip_latitude = latitude[chip_slice]

valid_geolocation = np.isfinite(chip_longitude) & np.isfinite(chip_latitude)
if not valid_geolocation.any():
    raise ValueError("The selected chip has no valid LUT geolocation.")

centre_lon = float(np.mean(chip_longitude[valid_geolocation]))
centre_lat = float(np.mean(chip_latitude[valid_geolocation]))
zone = int(np.floor((centre_lon + 180) / 6) + 1)
utm_epsg = (32600 if centre_lat >= 0 else 32700) + zone
utm_crs = CRS.from_epsg(utm_epsg)
to_utm = Transformer.from_crs("EPSG:4326", utm_crs, always_xy=True)
chip_x, chip_y = to_utm.transform(
    chip_longitude[valid_geolocation], chip_latitude[valid_geolocation]
)

xmin = np.floor(np.min(chip_x) / UTM_RESOLUTION_M) * UTM_RESOLUTION_M
ymin = np.floor(np.min(chip_y) / UTM_RESOLUTION_M) * UTM_RESOLUTION_M
xmax = np.ceil(np.max(chip_x) / UTM_RESOLUTION_M) * UTM_RESOLUTION_M
ymax = np.ceil(np.max(chip_y) / UTM_RESOLUTION_M) * UTM_RESOLUTION_M
target_width = int(round((xmax - xmin) / UTM_RESOLUTION_M))
target_height = int(round((ymax - ymin) / UTM_RESOLUTION_M))
target_transform = Affine.translation(xmin, ymax) * Affine.scale(
    UTM_RESOLUTION_M, -UTM_RESOLUTION_M
)

print(
    f"UTM target: EPSG:{utm_epsg}, {target_width} × {target_height} pixels at {UTM_RESOLUTION_M:g} m"
)
assert target_width > 0 and target_height > 0
assert (
    target_transform.a == UTM_RESOLUTION_M and target_transform.e == -UTM_RESOLUTION_M
)

UTM target: EPSG:32649, 1327 × 1289 pixels at 25 m

Checklist for the production workflow¶

Convert every Beta0 line/sample to the physical azimuth/slant-range coordinates specified by the annotation.
Locate those coordinates in the LUT vectors and bilinearly sample gammaNought, longitude, and latitude. Do not transpose or flip any LUT array.
Use NaN for internal nodata; calculate Gamma0 = Beta0_amplitude² × gammaNought before any map warp.
Use the sampled geometry/longitude and geometry/latitude arrays as the source geolocation. Do not fall back to the TIFF GCPs for this warp.
Build one explicit UTM target grid per output tile: CRS, snapped bounds, north-up transform, shape, and resolution. Reuse it across dates.
Warp Gamma0 directly from the LUT geolocation to the fixed UTM grid once with bilinear resampling. Write -9999.0 only at the COG boundary.

All processing helpers used above are defined in this notebook.

In [14]:

Copied!





gamma0_utm = np.full((target_height, target_width), np.nan, dtype="float32")
reproject(
    source=chip_gamma0,
    destination=gamma0_utm,
    src_geoloc_array=(chip_longitude, chip_latitude),
    src_crs=CRS.from_epsg(4326),
    dst_crs=utm_crs,
    dst_transform=target_transform,
    src_nodata=np.nan,
    dst_nodata=np.nan,
    resampling=Resampling.bilinear,
)

assert gamma0_utm.shape == (target_height, target_width)
assert np.isfinite(gamma0_utm).any(), "The UTM warp produced only nodata."

figure, axes = plt.subplots(1, 2, figsize=(14, 6), constrained_layout=True)
source_vmin, source_vmax = finite_display_limits(chip_gamma0)
source_image = axes[0].imshow(
    chip_gamma0, cmap="magma", vmin=source_vmin, vmax=source_vmax
)
axes[0].set(
    title="Calibrated Gamma0 chip in radar pixels", xlabel="sample", ylabel="line"
)
figure.colorbar(source_image, ax=axes[0])

utm_vmin, utm_vmax = finite_display_limits(gamma0_utm)
utm_image = axes[1].imshow(
    gamma0_utm,
    extent=(xmin, xmax, ymin, ymax),
    origin="upper",
    cmap="magma",
    vmin=utm_vmin,
    vmax=utm_vmax,
)
axes[1].set(
    title=f"The same chip after one LUT-geolocation warp to EPSG:{utm_epsg}",
    xlabel="UTM easting (m)",
    ylabel="UTM northing (m)",
)
axes[1].set_aspect("equal")
figure.colorbar(utm_image, ax=axes[1])
plt.show()
gamma0_utm = np.full((target_height, target_width), np.nan, dtype="float32")
reproject(
    source=chip_gamma0,
    destination=gamma0_utm,
    src_geoloc_array=(chip_longitude, chip_latitude),
    src_crs=CRS.from_epsg(4326),
    dst_crs=utm_crs,
    dst_transform=target_transform,
    src_nodata=np.nan,
    dst_nodata=np.nan,
    resampling=Resampling.bilinear,
)

assert gamma0_utm.shape == (target_height, target_width)
assert np.isfinite(gamma0_utm).any(), "The UTM warp produced only nodata."

figure, axes = plt.subplots(1, 2, figsize=(14, 6), constrained_layout=True)
source_vmin, source_vmax = finite_display_limits(chip_gamma0)
source_image = axes[0].imshow(
    chip_gamma0, cmap="magma", vmin=source_vmin, vmax=source_vmax
)
axes[0].set(
    title="Calibrated Gamma0 chip in radar pixels", xlabel="sample", ylabel="line"
)
figure.colorbar(source_image, ax=axes[0])

utm_vmin, utm_vmax = finite_display_limits(gamma0_utm)
utm_image = axes[1].imshow(
    gamma0_utm,
    extent=(xmin, xmax, ymin, ymax),
    origin="upper",
    cmap="magma",
    vmin=utm_vmin,
    vmax=utm_vmax,
)
axes[1].set(
    title=f"The same chip after one LUT-geolocation warp to EPSG:{utm_epsg}",
    xlabel="UTM easting (m)",
    ylabel="UTM northing (m)",
)
axes[1].set_aspect("equal")
figure.colorbar(utm_image, ax=axes[1])
plt.show()