depsi.arc_estimation

depsi.arc_estimation.periodogram(stm, key_dphase, key_h2ph, key_Btemporal, std_obs=1.0, std_height=50.0, std_vel=0.02, init_height=0.0, init_vel=0.0, init_step_height=1.0, init_step_vel=0.001, min_steps=10)

Periodogram algorithm.

This function performs periodogram unwrapping on arcs.

It uses a deformation model with two parameters: height and velocity to estimate the unwrapped phase.

For computation efficiency, the design matrix is constructed only once for all arcs, utilizing the average height-to-phase conversion factor (h2ph) across all arcs. The effect of using this average is corrected later.

Parameters:

Name	Type	Description	Default
stm	Dataset	Input Space-Time Matrix (STM) containing the wrapped phase, height-to-phase conversion factor, and year-time.	required
key_dphase	str	Key for the wrapped differential phase data variable in the STM.	required
key_h2ph	str	Key for the height-to-phase conversion factor in the STM.	required
key_Btemporal	str	Key for the temporal baseline in the STM. The value should be in decimal years.	required
std_obs	float	A-poriori standard deviation of the observations in rads, by default 1.0. This value is used to construct the stochastic model (Qyy) of the observations.	1.0
std_height	float	A-priori standard deviation of the height in meters, by default 50.0. This value is used to construct the boundaries of the initial search space for the height parameter.	50.0
std_vel	float	A-priori standard deviation of the velocity in meters per year, by default 0.02. This value is used to construct the boundaries of the initial search space for the velocity parameter.	0.02
init_height	float	Initial value for the height parameter in meters, by default 0.0.	0.0
init_vel	float	Initial value for the velocity parameter in meters per year, by default 0.0.	0.0
init_step_height	float	Initial step size for the height parameter in meters, by default 1.0. This value sets the resolution of the initial search space for the height parameter. After every search, the step size will be reduced by a factor of 10.	1.0
init_step_vel	float	Initial step size for the velocity parameter in meters per year, by default 1e-3. This value sets the resolution of the initial search space for the velocity parameter. After every search, the step size will be reduced by a factor of 10.	0.001
min_steps	int	Minimum number of steps in the search space for the height and velocity parameters, by default 10. If the number of steps in the initial search space is smaller than this value, it will be set to this value. After the first search, the number of steps will be set to this value.	10

Returns:

Type

Description

Tuple[DataArray, DataArray, DataArray, DataArray, DataArray]

Returns the unwrapped phase, ambiguities, estimated height, estimated velocity, and temporal coherence. - Unwrapped phase: in rads, shape (n_arcs, n_obs), dtype np.float64. - Ambiguities: unitless, shape (n_arcs, n_obs), dtype np.float64. - Estimated height: in meters, shape (n_arcs,), dtype np.float64. - Estimated velocity: in meters per year, shape (n_arcs,), dtype np.float64. - Temporal coherence: unitless float number, norm of the complex coherence, scalar, dtype np.float64.

Source code in depsi/arc_estimation.py

def periodogram(
    stm: xr.Dataset,
    key_dphase: str,
    key_h2ph: str,
    key_Btemporal: str,
    std_obs: float = 1.0,
    std_height: float = 50.0,
    std_vel: float = 0.02,
    init_height: float = 0.0,
    init_vel: float = 0.0,
    init_step_height: float = 1.0,
    init_step_vel: float = 1e-3,
    min_steps: int = 10,
):
    """Periodogram algorithm.

    This function performs periodogram unwrapping on arcs.

    It uses a deformation model with two parameters: height and velocity to estimate the unwrapped phase.

    For computation efficiency, the design matrix is constructed only once for all arcs, utilizing the average
    height-to-phase conversion factor (h2ph) across all arcs. The effect of using this average is corrected later.

    Parameters
    ----------
    stm : xr.Dataset
        Input Space-Time Matrix (STM) containing the wrapped phase, height-to-phase conversion factor, and year-time.
    key_dphase : str
        Key for the wrapped differential phase data variable in the STM.
    key_h2ph : str
        Key for the height-to-phase conversion factor in the STM.
    key_Btemporal : str
        Key for the temporal baseline in the STM.
        The value should be in decimal years.
    std_obs : float, optional
        A-poriori standard deviation of the observations in rads, by default 1.0.
        This value is used to construct the stochastic model (Qyy) of the observations.
    std_height : float, optional
        A-priori standard deviation of the height in meters, by default 50.0.
        This value is used to construct the boundaries of the initial search space for the height parameter.
    std_vel : float, optional
        A-priori standard deviation of the velocity in meters per year, by default 0.02.
        This value is used to construct the boundaries of the initial search space for the velocity parameter.
    init_height : float, optional
        Initial value for the height parameter in meters, by default 0.0.
    init_vel : float, optional
        Initial value for the velocity parameter in meters per year, by default 0.0.
    init_step_height : float, optional
        Initial step size for the height parameter in meters, by default 1.0.
        This value sets the resolution of the initial search space for the height parameter.
        After every search, the step size will be reduced by a factor of 10.
    init_step_vel : float, optional
        Initial step size for the velocity parameter in meters per year, by default 1e-3.
        This value sets the resolution of the initial search space for the velocity parameter.
        After every search, the step size will be reduced by a factor of 10.
    min_steps : int, optional
        Minimum number of steps in the search space for the height and velocity parameters, by default 10.
        If the number of steps in the initial search space is smaller than this value, it will be set to this value.
        After the first search, the number of steps will be set to this value.

    Returns
    -------
    Tuple[xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray]
        Returns the unwrapped phase, ambiguities, estimated height, estimated velocity, and temporal coherence.
        - Unwrapped phase: in rads, shape (n_arcs, n_obs), dtype np.float64.
        - Ambiguities: unitless, shape (n_arcs, n_obs), dtype np.float64.
        - Estimated height: in meters, shape (n_arcs,), dtype np.float64.
        - Estimated velocity: in meters per year, shape (n_arcs,), dtype np.float64.
        - Temporal coherence: unitless float number, norm of the complex coherence, scalar, dtype np.float64.
    """
    # Compute m2ph (meters to phase) conversion factor from wavelength
    if "wavelength" not in stm.attrs:
        raise ValueError(
            "Wavelength is not provided and not found in attributes of STM."
            "Please make sure it is provided."
            "For example: stm = stm.assign_attrs({'wavelength': wavelength})"
        )
    wavelength = stm.attrs["wavelength"]
    m2ph = -4 * np.pi / wavelength

    # Make sure year time only contains the time dimension
    assert (len(stm[key_Btemporal].dims) == 1) and (
        "time" in stm[key_Btemporal].dims
    ), "year time should and only should contain the 'time' dimension."

    # Load year time in memory
    Btemporal = stm[key_Btemporal].values

    # Set up functional and stochastic model for all arcs
    # Here we use the same h2ph (average over all arcs) for all arcs and correct the effect later
    # Doing this avoids perform matrix inversion for each arc
    h2ph_approx = stm[key_h2ph].mean(dim="space").values  # Mean h2ph of all arcs

    # Design matrix B, size n_obs x n_params
    # In B, h2ph should also be multiplied by m2ph since it did not when it was created
    B = np.stack([h2ph_approx * m2ph, Btemporal * m2ph]).T

    # Stochastic model Qyy, size n_obs x n_obs
    # This is the covariance matrix of the observations
    n_obs = stm[key_dphase].sizes["time"]  # number of observations
    Qyy = np.diag(np.repeat(std_obs**2, n_obs))

    # Normal matrix N and rhs for the least squares solution
    N = B.T @ np.linalg.inv(Qyy) @ B  # B.T * Qyy^-1 * B , size n_params x n_params
    # Solve N * x = B.T * Qyy^-1, then rhs =  N^-1 * B.T * Qyy^-1
    rhs = np.linalg.inv(N) @ B.T @ np.linalg.inv(Qyy)

    # Set up core dimensions, which are the dimensions _periodogram_arc will be applied to
    # We are broadcasting _periodogram_arc on stm[key_dphase] and stm[key_h2ph] along the space dimension
    # Therefore, we are calling it on the "time" dimension of every space entry.
    # So we have the input_core_dims as  [["time"], ["time"]]
    input_core_dims = [["time"], ["time"]]
    # There are 5 outputs from _periodogram_arc
    # The first two are np arrays with time dimension
    # The other three are scalars, so they have no dimensions
    output_core_dims = [["time"], ["time"], [], [], []]

    # Build initial search space for height and velocity
    n_steps_height = max(round(2 * std_height / init_step_height), min_steps)
    n_steps_vel = max(round(2 * std_vel / init_step_vel), min_steps)

    init_search_space = _build_periodogram_search_space(
        init_height, init_vel, init_step_height, init_step_vel, n_steps_height, n_steps_vel
    )

    # Apply the _periodogram_arc on stm[key_dphase] along "space" dimension
    # Other parameters are duplicated for each space entry
    # Therefore they can be passed as kwargs
    results = xr.apply_ufunc(
        _periodogram_arc,
        stm[key_dphase],
        stm[key_h2ph],
        input_core_dims=input_core_dims,
        output_core_dims=output_core_dims,
        kwargs={
            "h2ph_approx": h2ph_approx,
            "B": B,
            "Qyy": Qyy,
            "N": N,
            "rhs": rhs,
            "init_search_space": init_search_space,
            "init_step_height": init_step_height,
            "init_step_vel": init_step_vel,
            "min_steps": min_steps,
        },
        vectorize=True,
        dask="parallelized",
        output_dtypes=[np.float64, np.float64, np.float64, np.float64, np.float64],
    )

    return results