isce3 CSLC baseline computation

baselines.py

import datetime
from pyproj import CRS, Transformer

import h5py
import isce3
import numpy as np

from dolphin._types import Filename


def get_orbit_arrays(
    h5file: Filename,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, datetime.datetime]:
    """Parse orbit info from OPERA S1 CSLC HDF5 file into python types.

    Parameters
    ----------
    h5file : Filename
        Path to OPERA S1 CSLC HDF5 file.

    Returns
    -------
    times : np.ndarray
        Array of times in seconds since reference epoch.
    positions : np.ndarray
        Array of positions in meters.
    velocities : np.ndarray
        Array of velocities in meters per second.
    reference_epoch : datetime.datetime
        Reference epoch of orbit.
    """
    with h5py.File(h5file) as hf:
        orbit_group = hf["science/SENTINEL1/CSLC/metadata/orbit"]
        times = orbit_group["time"][:]
        positions = np.stack([orbit_group[f"position_{p}"] for p in ["x", "y", "z"]]).T
        velocities = np.stack([orbit_group[f"velocity_{p}"] for p in ["x", "y", "z"]]).T
        # take out ":26" once COMPASS fixes format
        reference_epoch = datetime.datetime.fromisoformat(
            orbit_group["reference_epoch"][()][:26].decode()
        )

    return times, positions, velocities, reference_epoch


def get_cslc_orbit(h5file: Filename) -> isce3.core.Orbit:
    """Parse orbit info from OPERA S1 CSLC HDF5 file into an isce3.core.Orbit.

    Parameters
    ----------

    h5file : Filename
        Path to OPERA S1 CSLC HDF5 file.

    Returns
    -------
    orbit : isce3.core.Orbit
        Orbit object.
    """
    times, positions, velocities, reference_epoch = get_orbit_arrays(h5file)
    orbit_svs = []
    for t, x, v in zip(times, positions, velocities):
        orbit_svs.append(
            isce3.core.StateVector(
                isce3.core.DateTime(reference_epoch + datetime.timedelta(seconds=t)),
                x,
                v,
            )
        )

    return isce3.core.Orbit(orbit_svs)


def _get_xy_grid(h5file: Filename, subsample: int = 100):
    """Get x and y grid from OPERA S1 CSLC HDF5 file.

    Parameters
    ----------
    h5file : Filename
        Path to OPERA S1 CSLC HDF5 file.
    subsample : int, optional
        Subsampling factor, by default 100

    Returns
    -------
    x : np.ndarray
        Array of x coordinates in meters.
    y : np.ndarray
        Array of y coordinates in meters.
    projection : int
        EPSG code of projection.
    """
    with h5py.File(h5file) as hf:
        x = hf["science/SENTINEL1/CSLC/grids/x_coordinates"][:]
        y = hf["science/SENTINEL1/CSLC/grids/y_coordinates"][:]
        projection = hf["science/SENTINEL1/CSLC/grids/projection"][()]

    return x[::subsample], y[::subsample], projection


def get_lonlat_grid(h5file: Filename, subsample: int = 100):
    """Get 2D latitude and longitude grid from OPERA S1 CSLC HDF5 file.

    Parameters
    ----------
    h5file : Filename
        Path to OPERA S1 CSLC HDF5 file.
    subsample : int, optional
        Subsampling factor, by default 100

    Returns
    -------
    lat : np.ndarray
        2D Array of latitude coordinates in degrees.
    lon : np.ndarray
        2D Array of longitude coordinates in degrees.
    projection : int
        EPSG code of projection.
    """

    x, y, projection = _get_xy_grid(h5file, subsample)
    X, Y = np.meshgrid(x, y)
    xx = X.flatten()
    yy = Y.flatten()
    crs = CRS.from_epsg(projection)
    transformer = Transformer.from_crs(crs, CRS.from_epsg(4326), always_xy=True)
    lon, lat = transformer.transform(xx=xx, yy=yy, radians=True)
    lon = lon.reshape(X.shape)
    lat = lat.reshape(Y.shape)
    return lon, lat


"""
geo2rdr(lon_lat_height: numpy.ndarray[numpy.float64[3, 1]], ellipsoid: isce3.ext.isce3.core.Ellipsoid = <isce3.ext.isce3.core.Ellipsoid object at 0x7fb6d83bbbf0>, orbit: isce3.ext.isce3.core.Orbit, doppler: isce3.ext.isce3.core.LUT2d, wavelength: float, side: object, threshold: float = 1e-08, maxiter: int = 50, delta_range: float = 10.0) -> Tuple[float, float]


This is the elementary transformation from map geometry to radar geometry.
The transformation is applicable for a single lon/lat/h coordinate (i.e., a
single point target).

Arguments:
    lon_lat_height  Lon/Lat/Hae of target of interest
    ellipsoid       Ellipsoid object
    orbit           Orbit object
    doppler         LUT2d Doppler model
    wavelength      Radar wavelength
    side            Left or Right
    threshold       azimuth time convergence threshold in meters
    maxiter         Maximum number of Newton-Raphson iterations
    delta_range     Step size used for computing derivative of doppler
    """


def _get_wavelength(h5file) -> float:
    with h5py.File(h5file) as hf:
        return hf[
            "science/SENTINEL1/CSLC/metadata/processing_information/s1_burst_metadata/wavelength"
        ][()]


def geo2rdr_cslc(
    h5file: Filename,
    threshold: float = 1e-08,
    maxiter: int = 50,
    delta_range: float = 10.0,
):
    """Get 2D latitude and longitude grid from OPERA S1 CSLC HDF5 file.

    Parameters
    ----------
    h5file : Filename
        Path to OPERA S1 CSLC HDF5 file.

    Returns
    -------
    lat : np.ndarray
        2D Array of latitude coordinates in degrees.
    lon : np.ndarray
        2D Array of longitude coordinates in degrees.
    projection : int
        EPSG code of projection.
    """
    lon_grid, lat_grid = get_lonlat_grid(h5file)
    lon_arr = lon_grid.ravel()
    lat_arr = lat_grid.ravel()
    height = 2000.0
    ellipsoid = isce3.core.Ellipsoid()
    orbit = get_cslc_orbit(h5file)
    zero_doppler = isce3.core.LUT2d()
    wavelength = _get_wavelength(h5file)
    side = isce3.core.LookSide.Right

    az_times, ranges = [], []
    for lon, lat in zip(lon_arr, lat_arr):
        llh_rad = np.array([lon, lat, height]).reshape((3, 1))
        az_time, range = isce3.geometry.geo2rdr(
            llh_rad,
            ellipsoid,
            orbit,
            zero_doppler,
            wavelength,
            side,
            threshold=threshold,
            maxiter=maxiter,
            delta_range=delta_range,
        )
        az_times.append(az_time)
        ranges.append(range)

    az_times = np.array(az_times).reshape(lon_grid.shape)
    ranges = np.array(ranges).reshape(lon_grid.shape)
    return az_times, ranges


def get_baselines(ref_h5file: Filename, sec_h5file: Filename) -> np.ndarray:
    """Get 2D baseline grid from OPERA S1 CSLC HDF5 file.

    Parameters
    ----------
    ref_h5file : Filename
        Path to reference OPERA S1 CSLC HDF5 file.
    sec_h5file : Filename
        Path to secondary OPERA S1 CSLC HDF5 file.

    Returns
    -------
    baseline : np.ndarray
        3D Array of baselines (X, Y, Z) in meters.
    """
    az_times_ref, ranges_ref = geo2rdr_cslc(ref_h5file)
    az_times_sec, ranges_sec = geo2rdr_cslc(sec_h5file)
    assert az_times_ref.shape == az_times_sec.shape

    orbit_ref = get_cslc_orbit(ref_h5file)
    orbit_sec = get_cslc_orbit(sec_h5file)
    baselines = np.zeros((3, *az_times_ref.shape))
    for i in range(az_times_ref.shape[0]):
        for j in range(az_times_ref.shape[1]):
            pos_ref = orbit_ref.interpolate(az_times_ref[i, j])[0]
            pos_sec = orbit_sec.interpolate(az_times_sec[i, j])[0]
            baselines[:, i, j] = pos_ref - pos_sec
    return baselines