depsi.network

depsi.network.form_network(stm, key_phase, key_h2ph, key_Btemporal, key_complex='complex', key_xcrds='lon', key_ycrds='lat', network_method='redundant', max_length=None, min_links=16, num_partitions=8, dphase_method='subtract')

Generate an STM of arcs from an STM of points.

Parameters:

Name	Type	Description	Default
stm	Dataset	Space-Time Matrix of scatterers.	required
key_phase	str	Key of the phase values in the STM. This phase will be used to compute the differential arc phase.	required
key_h2ph	str	Key of the h2ph values in the STM. The arc h2ph will be computed as the average between source and target.	required
key_Btemporal	str	Key of the temporal baseline values in the STM.	required
key_complex	str	Key of the complex values, by default "complex"	'complex'
key_xcrds	str	Key of the x coordinates for calulating arc length, by default "lon"	'lon'
key_ycrds	str	Key of the y coordinates for calulating arc length, by default "lat"	'lat'
network_method	Literal['redundant', 'delaunay']	network formation method, by default "redundant"	'redundant'
max_length	float	maximum arc length, by default None	None
min_links	int	minimum links per point, by default 16 only effective when network_method is "redundant"	16
num_partitions	int	number of partitions of searching when forming redundant network, by default 8 only effective when network_method is "redundant"	8
dphase_method	Literal['conjmult', 'subtract']	method of computing phase difference, by default "subtract" "subtract" method subtracts the source phase from the target phase (without re-wrapping); "conjmult" method computes the phase difference by conjugate multiplication: d_phase = np.angle(complex_target * complex_source.conj())	'subtract'

Returns:

Type	Description
Dataset	Space-Time Matrix of arcs, containing the following variables: - d_phase: the arc phase, which is the difference between source and target points - Btemp: the temporal baseline, which is the same for all arcs - h2ph: the arc h2ph, which is the average between source and target points

Source code in depsi/network.py

def form_network(
    stm: xr.Dataset,
    key_phase: str,
    key_h2ph: str,
    key_Btemporal: str,
    key_complex: str = "complex",
    key_xcrds: str = "lon",
    key_ycrds: str = "lat",
    network_method: Literal["redundant", "delaunay"] = "redundant",
    max_length: float = None,
    min_links: int = 16,
    num_partitions: int = 8,
    dphase_method: Literal["conjmult", "subtract"] = "subtract",
) -> xr.Dataset:
    """Generate an STM of arcs from an STM of points.

    Parameters
    ----------
    stm : xr.Dataset
        Space-Time Matrix of scatterers.
    key_phase : str
        Key of the phase values in the STM.
        This phase will be used to compute the differential arc phase.
    key_h2ph : str
        Key of the h2ph values in the STM.
        The arc h2ph will be computed as the average between source and target.
    key_Btemporal : str
        Key of the temporal baseline values in the STM.
    key_complex : str, optional
        Key of the complex values, by default "complex"
    key_xcrds  : str, optional
        Key of the x coordinates for calulating arc length, by default "lon"
    key_ycrds  : str, optional
        Key of the y coordinates for calulating arc length, by default "lat"
    network_method : Literal["redundant", "delaunay"], optional
        network formation method, by default "redundant"
    max_length : float, optional
        maximum arc length, by default None
    min_links : int, optional
        minimum links per point, by default 16
        only effective when network_method is "redundant"
    num_partitions : int, optional
        number of partitions of searching when forming redundant network, by default 8
        only effective when network_method is "redundant"
    dphase_method : Literal["conjmult", "subtract"], optional
        method of computing phase difference, by default "subtract"
        "subtract" method subtracts the source phase from the target phase (without re-wrapping);
        "conjmult" method computes the phase difference by conjugate multiplication:
            d_phase = np.angle(complex_target * complex_source.conj())

    Returns
    -------
    xr.Dataset
        Space-Time Matrix of arcs, containing the following variables:
        - d_phase: the arc phase, which is the difference between source and target points
        - Btemp: the temporal baseline, which is the same for all arcs
        - h2ph: the arc h2ph, which is the average between source and target points
    """
    # Generate the network arcs.
    if network_method == "redundant":
        if min_links <= 0:
            logger.error(f"min_links must be strictly positive (currently: {min_links})")
            return
        if num_partitions <= 0:
            logger.error(f"num_partitions must be strictly positive (currently: {num_partitions})")
            return
    elif network_method != "delaunay":
        raise NotImplementedError(f"Unknown network method {network_method}, known are delaunay and redundant")

    # Collect point coordinates.
    indices = [stm[coord] for coord in [key_xcrds, key_ycrds]]
    coordinates = np.column_stack(indices)

    arcs = None

    # Create network arcs as list of tuples of point ids.
    if network_method == "delaunay":
        arcs = _generate_arcs_delaunay(coordinates, max_length)
    elif network_method == "redundant":
        arcs = _generate_arcs_redundant(coordinates, max_length, min_links, num_partitions)

    # Compute the phase difference.
    arcs_unzipped = list(zip(*arcs, strict=False))
    source_idx = list(arcs_unzipped[0])
    target_idx = list(arcs_unzipped[1])

    if dphase_method not in ["conjmult", "subtract"]:
        raise NotImplementedError(f"Unknown dphase_method '{dphase_method}'.")
    dict_key_method = {"subtract": key_phase, "conjmult": key_complex}  # mapping for selecting the correct key
    d_phase = compute_phase_difference(
        stm.isel(space=source_idx)[dict_key_method[dphase_method]].data,
        stm.isel(space=target_idx)[dict_key_method[dphase_method]].data,
        method=dphase_method,
    )

    # Temporal baseline
    Btemp = stm[key_Btemporal].data

    # Height to phase factor
    h2ph = (stm[key_h2ph].isel(space=source_idx).data + stm[key_h2ph].isel(space=target_idx).data) / 2

    # Generate a unique identifier of arcs based on source and target for easy indexing
    # This is because when updating network, points can be removed and reindexed
    # Therefore we cannot use 2d index (source, target) as uid
    # NOTE: This encoding is safe for typical networks (<10,000 points, <100,000 arcs).
    # For 10,000 points, max UID ≈ 1e9, well within int64 max (2^63-1 ≈ 9.22e18), providing
    # a safety margin of >9 billion times. Even for 1M points, safety margin is >900,000x.
    scale = 10 ** (math.floor(math.log10(stm.sizes["space"])) + 1)  # Scale to ensure no overlap
    uid = scale * (np.array(source_idx) + 1) + (np.array(target_idx) + 1)  # Plus one to avoid zero uid
    uid = uid.astype(np.int64)

    arcs = xr.Dataset(
        data_vars={
            "d_phase": (["space", "time"], d_phase),
            "Btemp": (["time"], Btemp),
            "h2ph": (["space", "time"], h2ph),
        },
        coords={"source": (["space"], source_idx), "target": (["space"], target_idx), "uid": (["space"], uid)},
        attrs=stm.attrs,
    )

    return arcs

depsi.network.spatial_integration(stm_pnts, stm_arcs, key_sdphase='sd_phase', key_arc_quality='temp_coh', threshold_arc_quality=0.5, idx_refpnt=None, min_arc_connections=3, parallel=False, sparse_mode=False, ensure_network_while_mht=False, arc_estimation_method='periodogram')

Spatially integrate the ambiguities of network arcs to points.

This function estimates the integer ambiguities of the points from arc ambiguities. It assumes a network has been formed from the points by calling "network.form_network", and arc ambiguities have been estimated by calling relevant functions in the "depsi.arc_estimation" module.

The function returns an updated STM of arcs which contains the adjusted arc ambiguities after spatial integration, as well as an updated STM of points which contains the estimated point ambiguities.

The following steps are performed: 1. Validate the network arcs (stm_arcs) against the point STM (stm_pnts) 2. Select arcs based on quality threshold and ensure minimum connections for all points 3. Select a reference point which assumes zero phase (hence zero ambiguity) 4. Adjust the network by removing arcs/points which potentially cause errors using Multi-Hypothesis Testing (MHT) 5. Adjust the ambiguities per time epoch to make sure spatial solutions give zero residuals. 6. Calculate point ambiguities and unwrapped phases w.r.t. the reference point.

Parameters:

Name	Type	Description	Default
stm_pnts	Dataset	Space-Time Matrix of points.	required
stm_arcs	Dataset	Space-Time Matrix of arcs. The arcs should be formed from stm_pnts, using the "network.form_network" function. This function generates coordinates "source" and "target" in stm_arcs which refer to the indices in stm_pnts. An "uid" data variable is also generated by "form_network" for easy indexing. An arc estimation should have been applied on stm_arcs before calling "spatial_integration". One can use relevant functions in "depsi.arc_estimation" module for this purpose. Arc estimation adds the variable "ambiguities" to stm_arcs, which are the estimated arc ambiguities. It also adds quality variables such as "temp_coh" (ensemble coherence), which are used to select arcs for spatial integration.	required
key_sdphase	str	Key of the single difference phase variable in stm_pnts, by default "sd_phase" This phase is used to compute unwrapped phases after ambiguity estimation.	'sd_phase'
key_arc_quality	str	Key of the arc quality variable in stm_arcs, by default "temp_coh"	'temp_coh'
threshold_arc_quality	float	Threshold for arc quality, by default 0.5	0.5
idx_refpnt	int | None	Index of the reference point in stm_pnts. If None, the source point of the arc with highest quality is selected as the reference point.	None
min_arc_connections	int	Minimum number of connections for arcs, by default 3	3
parallel	bool	Whether to use parallel processing, by default False	False
sparse_mode	bool	Whether to use sparse matrix format for large networks, by default False	False
ensure_network_while_mht	bool	Whether to ensure minimum connections in MHT network adjustment, by default False	False
arc_estimation_method	Literal['periodogram']	Method used for arc estimation, by default "periodogram". This constrains the method used for VCM computation.	'periodogram'

Returns:

Type	Description
(Dataset, Dataset)	Updated Space-Time Matrix of arcs and updated Space-Time Matrix of points. For arcs, the "ambiguities" variable contains the adjusted arc ambiguities after spatial integration. For points, the "ambiguities" variable contains the estimated point ambiguities, and "unwrapped_phase" contains the unwrapped phase w.r.t. the reference point.

References

Van Leijen, F.J.. "Persistent scatterer interferometry based on geodetic estimation theory." (2014).

Source code in depsi/network.py

def spatial_integration(
    stm_pnts: xr.Dataset,
    stm_arcs: xr.Dataset,
    key_sdphase: str = "sd_phase",
    key_arc_quality: str = "temp_coh",
    threshold_arc_quality: float = 0.5,
    idx_refpnt: int | None = None,
    min_arc_connections: int = 3,
    parallel: bool = False,
    sparse_mode: bool = False,
    ensure_network_while_mht: bool = False,
    arc_estimation_method: Literal["periodogram"] = "periodogram",
) -> tuple[xr.Dataset, xr.Dataset]:
    """Spatially integrate the ambiguities of network arcs to points.

    This function estimates the integer ambiguities of the points from arc ambiguities. It assumes a network
    has been formed from the points by calling "network.form_network", and arc ambiguities have been estimated
    by calling relevant functions in the "depsi.arc_estimation" module.

    The function returns an updated STM of arcs which contains the adjusted arc ambiguities after spatial integration,
    as well as an updated STM of points which contains the estimated point ambiguities.

    The following steps are performed:
    1. Validate the network arcs (`stm_arcs`) against the point STM (`stm_pnts`)
    2. Select arcs based on quality threshold and ensure minimum connections for all points
    3. Select a reference point which assumes zero phase (hence zero ambiguity)
    4. Adjust the network by removing arcs/points which potentially cause errors using Multi-Hypothesis Testing (MHT)
    5. Adjust the ambiguities per time epoch to make sure spatial solutions give zero residuals.
    6. Calculate point ambiguities and unwrapped phases w.r.t. the reference point.

    Parameters
    ----------
    stm_pnts : xr.Dataset
        Space-Time Matrix of points.
    stm_arcs : xr.Dataset
        Space-Time Matrix of arcs.
        The arcs should be formed from stm_pnts, using the "network.form_network" function. This function generates
        coordinates "source" and "target" in stm_arcs which refer to the indices in stm_pnts. An "uid" data variable
        is also generated by "form_network" for easy indexing.
        An arc estimation should have been applied on stm_arcs before calling "spatial_integration". One can use
        relevant functions in "depsi.arc_estimation" module for this purpose. Arc estimation adds the variable
        "ambiguities" to stm_arcs, which are the estimated arc ambiguities. It also adds quality variables such as
        "temp_coh" (ensemble coherence), which are used to select arcs for spatial integration.
    key_sdphase : str, optional
        Key of the single difference phase variable in stm_pnts, by default "sd_phase"
        This phase is used to compute unwrapped phases after ambiguity estimation.
    key_arc_quality : str, optional
        Key of the arc quality variable in stm_arcs, by default "temp_coh"
    threshold_arc_quality : float, optional
        Threshold for arc quality, by default 0.5
    idx_refpnt : int | None, optional
        Index of the reference point in stm_pnts. If None, the source point of the arc with highest quality is selected
        as the reference point.
    min_arc_connections : int, optional
        Minimum number of connections for arcs, by default 3
    parallel : bool, optional
        Whether to use parallel processing, by default False
    sparse_mode : bool, optional
        Whether to use sparse matrix format for large networks, by default False
    ensure_network_while_mht : bool, optional
        Whether to ensure minimum connections in MHT network adjustment, by default False
    arc_estimation_method : Literal["periodogram"], optional
        Method used for arc estimation, by default "periodogram".
        This constrains the method used for VCM computation.

    Returns
    -------
    xr.Dataset, xr.Dataset
        Updated Space-Time Matrix of arcs and updated Space-Time Matrix of points.
        For arcs, the "ambiguities" variable contains the adjusted arc ambiguities after spatial integration.
        For points, the "ambiguities" variable contains the estimated point ambiguities, and "unwrapped_phase"
        contains the unwrapped phase w.r.t. the reference point.

    References
    ----------
    Van Leijen, F.J.. "Persistent scatterer interferometry based on geodetic estimation theory." (2014).
    """
    # Check parallelization behavior
    if parallel:
        raise NotImplementedError("Dask support is not implemented yet for spatial_integration.")
    else:
        # Compute all data into memory
        stm_pnts = stm_pnts.compute()
        stm_arcs = stm_arcs.compute()

    # Validate that stm_arcs are formed from stm_pnts
    if (stm_arcs["source"].max().values >= stm_pnts.sizes["space"]) or (
        stm_arcs["target"].max().values >= stm_pnts.sizes["space"]
    ):
        raise ValueError("stm_arcs contains source/target indices that exceed the number of points in stm_pnts.")
    if "ambiguities" not in stm_arcs:
        raise ValueError("stm_arcs does not contain 'ambiguities' variable. Please estimate arc ambiguities first.")
    if key_arc_quality not in stm_arcs:
        raise ValueError(f"stm_arcs does not contain '{key_arc_quality}' variable for arc quality assessment.")

    # Check arc estimation method, this constrains VCM computation method
    if arc_estimation_method not in ["periodogram"]:
        raise NotImplementedError(f"Unknown arc estimation method {arc_estimation_method}.")

    if sparse_mode:
        raise NotImplementedError("Sparse mode is not implemented yet for spatial_integration.")

    # If idx_refpnt is specified
    # Get radar coordinates of the reference point before any shape change
    if idx_refpnt is not None:
        azimuth_refpnt = stm_pnts["azimuth"].isel(space=idx_refpnt).values
        range_refpnt = stm_pnts["range"].isel(space=idx_refpnt).values

    # Select arcs with quality > threshold_arc_quality
    # Then ensure all points have at least min_arc_connections connections
    mask = (np.abs(stm_arcs[key_arc_quality]) > threshold_arc_quality).compute()
    stm_arcs = stm_arcs.where(mask, drop=True)
    stm_arcs, stm_pnts = _ensure_network_min_connections(stm_arcs, stm_pnts, min_arc_connections)

    # Select reference point as the source pnt of arcs with highest temp_coh
    if idx_refpnt is None:
        idx_arc_max_coh = stm_arcs[key_arc_quality].argmax().values
        idx_refpnt = stm_arcs["source"].isel(space=idx_arc_max_coh).values
        azimuth_refpnt = stm_pnts["azimuth"].isel(space=idx_refpnt).values
        range_refpnt = stm_pnts["range"].isel(space=idx_refpnt).values
    else:
        # Make sure idx_refpnt is still valid after arc selection and point removal
        mask_refpnt = (stm_pnts["azimuth"].values == azimuth_refpnt) & (stm_pnts["range"].values == range_refpnt)
        if not np.any(mask_refpnt):
            raise ValueError(
                f"Reference point ({azimuth_refpnt}, {range_refpnt}) removed after arc selection. "
                f"Please choose another reference point."
            )

        # Get reference point from radar coordinates
        idx_refpnt = np.where(mask_refpnt)[0][0]

    # Adjust the network by removing bad arcs/points using MHT
    stm_arcs_adjusted, stm_pnts_adjusted = _mht_network_adjustment(
        stm_arcs,
        stm_pnts,
        idx_refpnt,
        azimuth_refpnt,
        range_refpnt,
        ensure_network_while_mht,
        sparse_mode,
        arc_estimation_method,
    )

    # Update idx_refpnt after MHT adjustment
    idx_refpnt = np.where(
        (stm_pnts_adjusted["azimuth"].values == azimuth_refpnt) & (stm_pnts_adjusted["range"].values == range_refpnt)
    )[0][0]

    # Adjust ambiguities to fix unwrapping errors
    stm_arcs_output, stm_pnts_output, idx_refpnt = _ambiguities_adjustment(
        stm_arcs_adjusted,
        stm_pnts_adjusted,
        idx_refpnt,
        sparse_mode,
        arc_estimation_method,
    )

    # Assign idx_refpnt as attribute to stm_pnts_output
    stm_pnts_output = stm_pnts_output.assign_attrs({"idx_refpnt": idx_refpnt})

    # Add unwrapped phase to stm_pnts_output
    # Unwrapped phase is w.r.t. the reference point
    # Therefore the sd_phase of the reference point is subtracted
    unwrapped_phase_pnts = (
        stm_pnts_output[key_sdphase].data
        + stm_pnts_output["ambiguities"].data * 2 * np.pi
        - np.tile(
            stm_pnts_output[key_sdphase].isel(space=idx_refpnt).data,
            (stm_pnts_output.sizes["space"], 1),
        )
    )
    stm_pnts_output["unwrapped_phase"] = (("space", "time"), unwrapped_phase_pnts)

    return stm_arcs_output, stm_pnts_output