Velocity from a trajectory fit¶

fit_trajectory fits every pixel with a trajectory model: offset and velocity, seasonal sine/cosine pairs, steps and postseismic terms, by least squares over each pixel's valid epochs. The synthetic frame has a known velocity, a 5 mm annual cycle and a 10 mm step, so the fit can be checked.

In [1]:

import tempfile
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

import disp_xr as dxr
from disp_xr.tutorial import make_disp_stack

folder = Path(tempfile.mkdtemp()) / "granules"
truth = make_disp_stack(
    folder, "three_ministacks", shape=(64, 64), epochs_per_ministack=70
)
ds = dxr.open_disp(folder).disp.mask("recommended").disp.chain()
print(
    f"{ds.sizes['time']} epochs over "
    f"{(ds.time[-1] - ds.time[0]).values.astype('timedelta64[D]')}"
)

import tempfile from pathlib import Path import matplotlib.pyplot as plt import numpy as np import disp_xr as dxr from disp_xr.tutorial import make_disp_stack folder = Path(tempfile.mkdtemp()) / "granules" truth = make_disp_stack( folder, "three_ministacks", shape=(64, 64), epochs_per_ministack=70 ) ds = dxr.open_disp(folder).disp.mask("recommended").disp.chain() print( f"{ds.sizes['time']} epochs over " f"{(ds.time[-1] - ds.time[0]).values.astype('timedelta64[D]')}" )

211 epochs over 2520 days

Fit¶

Trajectory lists the terms. The coefficient names encode the model, so a fit written to a store can be evaluated again with predict_trajectory.

In [2]:

step_date = f"{truth.step_date:%Y-%m-%d}"
terms = dxr.Trajectory(polynomial=1, periods=(1.0,), steps=(step_date,))
fit = dxr.fit_trajectory(ds, terms).compute()
print([str(name) for name in fit.coefficient.values])
print(list(fit.data_vars))

step_date = f"{truth.step_date:%Y-%m-%d}" terms = dxr.Trajectory(polynomial=1, periods=(1.0,), steps=(step_date,)) fit = dxr.fit_trajectory(ds, terms).compute() print([str(name) for name in fit.coefficient.values]) print(list(fit.data_vars))

['constant', 'linear_trend', 'annual_sin', 'annual_cos', 'step_2023-06-14']
['coefficients', 'sigma', 'rms', 'n_valid', 'velocity', 'velocity_sigma', 'annual_amplitude', 'annual_amplitude_sigma']

Velocity and its uncertainty¶

These sigmas assume white noise; the next section compares them with colored-noise and robust estimates.

In [3]:

fig, (left, right) = plt.subplots(1, 2, figsize=(11, 4.2), layout="constrained")
(fit.velocity * 1000).plot.imshow(
    ax=left, cmap="RdBu_r", cbar_kwargs={"label": "mm/yr"}
)
left.set_title("velocity", loc="left")
(fit.velocity_sigma * 1000).plot.imshow(ax=right, cbar_kwargs={"label": "mm/yr"})
right.set_title("velocity sigma (white noise)", loc="left")
for ax in (left, right):
    ax.set_aspect("equal")
plt.show()

fig, (left, right) = plt.subplots(1, 2, figsize=(11, 4.2), layout="constrained") (fit.velocity * 1000).plot.imshow( ax=left, cmap="RdBu_r", cbar_kwargs={"label": "mm/yr"} ) left.set_title("velocity", loc="left") (fit.velocity_sigma * 1000).plot.imshow(ax=right, cbar_kwargs={"label": "mm/yr"}) right.set_title("velocity sigma (white noise)", loc="left") for ax in (left, right): ax.set_aspect("equal") plt.show()

Against the truth¶

In [4]:

error = (fit.velocity - truth.velocity) / fit.velocity_sigma
fitted = error.notnull()
print(
    f"velocity within 2 sigma: {float((abs(error) < 2).sum() / fitted.sum()):.0%}"
    " of fitted pixels"
)
print(
    f"annual amplitude: {float(fit.annual_amplitude.median()) * 1000:.2f} mm"
    f" (truth {truth.annual_amplitude * 1000:.0f} mm)"
)
step = fit.coefficients.sel(coefficient=f"step_{step_date}")
print(f"step: {float(step.median()) * 1000:.2f} mm (truth {truth.step * 1000:.0f} mm)")

error = (fit.velocity - truth.velocity) / fit.velocity_sigma fitted = error.notnull() print( f"velocity within 2 sigma: {float((abs(error) < 2).sum() / fitted.sum()):.0%}" " of fitted pixels" ) print( f"annual amplitude: {float(fit.annual_amplitude.median()) * 1000:.2f} mm" f" (truth {truth.annual_amplitude * 1000:.0f} mm)" ) step = fit.coefficients.sel(coefficient=f"step_{step_date}") print(f"step: {float(step.median()) * 1000:.2f} mm (truth {truth.step * 1000:.0f} mm)")

velocity within 2 sigma: 95% of fitted pixels

annual amplitude: 5.00 mm (truth 5 mm)
step: 10.01 mm (truth 10 mm)

Model and residuals at a pixel¶

predict_trajectory evaluates the fit at any dates, lazily and in float32; the model is never stored unless you write it.

In [5]:

model = dxr.predict_trajectory(fit, ds.time.values)
row, col = 30, 50
fig, ax = plt.subplots(figsize=(8, 3.2), layout="constrained")
ax.plot(
    ds.time,
    ds.displacement[:, row, col] * 1000,
    "o",
    ms=3,
    color="#2a78d6",
    label="data",
)
ax.plot(ds.time, model[:, row, col] * 1000, color="#e34948", lw=2, label="trajectory")
ax.set_ylabel("displacement (mm)")
ax.legend(frameon=False)
ax.set_title(f"row {row}, col {col}", loc="left")
plt.show()
residual = ds.displacement[:, row, col] - model[:, row, col]
print(f"residual RMS: {float(np.sqrt((residual**2).mean())) * 1000:.2f} mm")

model = dxr.predict_trajectory(fit, ds.time.values) row, col = 30, 50 fig, ax = plt.subplots(figsize=(8, 3.2), layout="constrained") ax.plot( ds.time, ds.displacement[:, row, col] * 1000, "o", ms=3, color="#2a78d6", label="data", ) ax.plot(ds.time, model[:, row, col] * 1000, color="#e34948", lw=2, label="trajectory") ax.set_ylabel("displacement (mm)") ax.legend(frameon=False) ax.set_title(f"row {row}, col {col}", loc="left") plt.show() residual = ds.displacement[:, row, col] - model[:, row, col] print(f"residual RMS: {float(np.sqrt((residual**2).mean())) * 1000:.2f} mm")

residual RMS: 0.96 mm

Colored-noise and robust sigmas¶

fit_velocity estimates the noise of every pixel (see "Choosing a velocity method" for when to use which). "whittle" keeps the OLS velocity and estimates power-law plus white noise from the periodogram of the residuals (reliable from about 200 epochs); "midas" takes the median of one-year slopes, with a deliberately conservative sigma. The synthetic noise is white, so OLS and Whittle should agree and Whittle should find kappa near 0; on real InSAR series the white-noise sigma is several times too small.

In [6]:

for method in ("ols", "whittle", "midas"):
    result = dxr.fit_velocity(ds, method, terms=terms).compute()
    sigma = float(result.velocity_sigma.median()) * 1000
    print(f"{method:>8}: median velocity sigma {sigma:.2f} mm/yr")
whittle = dxr.fit_velocity(ds, "whittle", terms=terms).compute()
print(f"median kappa: {float(whittle.kappa.median()):.2f} (white noise is 0)")

for method in ("ols", "whittle", "midas"): result = dxr.fit_velocity(ds, method, terms=terms).compute() sigma = float(result.velocity_sigma.median()) * 1000 print(f"{method:>8}: median velocity sigma {sigma:.2f} mm/yr") whittle = dxr.fit_velocity(ds, "whittle", terms=terms).compute() print(f"median kappa: {float(whittle.kappa.median()):.2f} (white noise is 0)")

     ols: median velocity sigma 0.07 mm/yr

 whittle: median velocity sigma 0.07 mm/yr

   midas: median velocity sigma 0.62 mm/yr

median kappa: 0.00 (white noise is 0)