Skip to content

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=3.0, init_step_vel=0.002, 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 3.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.

 3.0
init_step_vel float

Initial step size for the velocity parameter in meters per year, by default 2e-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.002
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 
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
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 = 3.0,
    init_step_vel: float = 2e-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 3.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 2e-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)

    # check if the time dimension is not chunked, and unchunk it if necessary
    if "time" in stm.chunks.keys():
        if len(stm.chunks["time"]) != 1:
            stm = stm.chunk({"time": -1})

    # 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
    )

    # Perform one search for all arcs, and get the best initial estimates for height and velocity per arc
    # This is motivated by the fact that the initial search space is the largest, and can be vectorized for all arcs
    # First iteration, candidate modeled phases are identical for all arcs
    # The residuals phase_residual_all_arcs is a large array with n_arcs x n_obs x n_search
    # so use .data to avoid loading it into memory if it is a dask array
    dphase_obs = stm[key_dphase].data  # n_arcs x n_obs x 1

    # If the memory size of phase_residual_all_arcs will exceed the threshold
    # chunk dphase_obs and init_search_space to enable computation
    mem_size_estimation = (
        dphase_obs.shape[0] * dphase_obs.shape[1] * init_search_space.shape[0] * dphase_obs.dtype.itemsize
    ) / (1024**2)  # in MB
    if mem_size_estimation > THRES_TEMP_COH_MEMORY:
        dphase_obs, init_search_space = _chunk_for_temp_coh_compute(dphase_obs, init_search_space)

    # Compute modelled phase for all arcs and all search candidates
    phs_model = B @ init_search_space.T  # n_obs x n_search

    # Expand dimensions and compute the phase residuals for all arcs and all search candidates
    dphase_obs = dphase_obs[:, :, None]  # n_arcs x n_obs x 1
    phs_model = phs_model[None, :, :]  # 1 x n_obs x n_search
    phase_residual_all_arcs = dphase_obs - phs_model

    # Find the best initial height and velocity based on the temporal coherence
    coh_search_space_all_arcs = (
        np.cos(phase_residual_all_arcs).sum(axis=1) + 1j * np.sin(phase_residual_all_arcs).sum(axis=1)
    ) / stm[key_dphase].sizes["time"]  # n_arcs x n_search
    coh_idx_all_arcs = np.argmax(np.abs(coh_search_space_all_arcs), axis=1)  # n_arcs

    # Implicitly compute best coh index if dask array
    coh_idx_all_arcs = coh_idx_all_arcs.compute() if isinstance(coh_idx_all_arcs, da.Array) else coh_idx_all_arcs

    # Build xr.DataArray for the initial height and velocity of all arcs
    da_init_height_all_arcs = xr.DataArray(
        init_search_space[coh_idx_all_arcs, 0],
        dims=["space"],
    )
    da_init_vel_all_arcs = xr.DataArray(
        init_search_space[coh_idx_all_arcs, 1],
        dims=["space"],
    )

    # Apply the _periodogram_arc on stm[key_dphase] along "space" dimension
    # Set up input 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
    # The height and velocity are scalars
    # Therefore, we are only calling it on the "time" dimension for the first two parameters
    # 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"], [], [], []]

    results = xr.apply_ufunc(
        _periodogram_arc,
        stm[key_dphase],
        stm[key_h2ph],
        da_init_height_all_arcs,
        da_init_vel_all_arcs,
        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_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