Spatial models¶

Estimators with a fit / predict interface for spatial models: trends, the GNSS tie (T33) will build on them. verde is the design reference, not a dependency.

from disp_xr import models

flat = ds.disp.deramp(degree=1)  # a plane removed from every epoch, lazily

coords = models.coordinates(points)  # (easting, northing) of a point table
est = models.Chain([("trend", models.Trend(1)), ("spline", models.Spline(1e-6))])
est.fit(coords, points.velocity.values)
grid = est.grid(like=ds)  # on the frame's grid, with its CRS
scores = models.cross_val_score(est, coords, points.velocity.values, spacing=5000)

Estimator base for spatial models (PRD F9, plan task T19).

Spatial models (trends, GNSS ties, gridders) share one shape, after verde (the design reference only, not a dependency; PRD decision):

fit(coords, data, weights=None) estimates the parameters and returns the estimator; fitted attributes end in _;
predict(coords) evaluates the model at any coordinates;
score(coords, data, weights=None) is the weighted R²;
grid(...) predicts on a regular grid as an xr.DataArray with a CRS;
:class:Chain fits steps in sequence, each on the residual of the last.

Coordinates are a tuple (easting, northing) of arrays of one shape, in the CRS of the data (projected, in meters for DISP frames). :func:coordinates builds them from a raster, a point table or a DataFrame. :func:cross_val_score scores an estimator with spatial block k-fold cross-validation, which keeps neighbors out of each other's folds (Roberts et al. 2017, https://doi.org/10.1111/ecog.02881).

BaseEstimator ¶

Base class of disp-xr spatial estimators.

Subclasses take their parameters as keyword arguments of __init__ (stored under the same names) and implement _fit(easting, northing, data, weights) and _predict(easting, northing) on 1D arrays without NaN; the base class checks inputs, drops NaN and reshapes.

Examples:

>>> class Offset(BaseEstimator):
...     def fit(self, coords, data, weights=None): ...
...     def predict(self, coords): ...

get_params ¶

get_params() -> dict[str, Any]

Return the constructor parameters.

set_params ¶

set_params(**params: Any) -> BaseEstimator

Set constructor parameters; returns the estimator.

fit ¶

fit(coords: Coords, data: Any, weights: Optional[Any] = None) -> BaseEstimator

Fit the model; returns the estimator.

Parameters:

Name	Type	Description	Default
`coords`	`(easting, northing)`	Arrays of one shape.	required
`data`	`array - like`	Values at `coords` (same shape); NaN points are skipped.	required
`weights`	`array - like`	Non-negative weights (same shape), for example `1 / sigma**2`; points with a zero or NaN weight are skipped.	`None`

predict ¶

predict(coords: Coords) -> np.ndarray

Return the model at coords (same shape; NaN where they are NaN).

score ¶

score(coords: Coords, data: Any, weights: Optional[Any] = None) -> float

Return the weighted R² of the predictions at coords.

grid ¶

grid(*, region: Optional[Sequence[float]] = None, spacing: Optional[float] = None, like: Optional[Union[Dataset, DataArray]] = None, crs: Optional[Any] = None, name: str = 'prediction') -> xr.DataArray

Predict on a regular grid.

Parameters:

Name	Type	Description	Default
`region`	`(west, east, south, north)`	Bounds of the pixel centers; default the fitted data's bounds.	`None`
`spacing`	`float`	Pixel size, needed unless `like` is given.	`None`
`like`	`Dataset or DataArray`	Predict on this grid (its `x`, `y` and CRS).	`None`
`crs`	`optional`	CRS to attach (default: the CRS of `like`).	`None`
`name`	`str`	Name of the returned DataArray.	`'prediction'`

Returns:

Type	Description
`DataArray`	`(y, x)`, north up, with the CRS when known.

Chain ¶

Bases: BaseEstimator

Fit estimators in sequence, each on the residual of the previous ones.

predict sums the steps. For example, a planar trend then a gridder of what the trend leaves.

Parameters:

Name	Type	Description	Default
`steps`	`sequence of (name, estimator)`	The estimators, in order.	required

Examples:

>>> from disp_xr.models import Chain, Spline, Trend
>>> chain = Chain([("trend", Trend(1)), ("spline", Spline(damping=1e-2))])
>>> chain.fit(coords, data).predict(coords)

init ¶

__init__(steps: Sequence[tuple[str, BaseEstimator]]) -> None

Store the steps.

Mean ¶

Bases: BaseEstimator

The weighted mean of the data: the simplest model (a baseline).

Examples:

>>> Mean().fit(coords, data).predict(coords)

init ¶

__init__() -> None

Take no parameters.

Trends and de-ramping¶

Polynomial trends and per-epoch de-ramping (PRD F9, plan task T24).

:class:Trend fits a 2D polynomial of the coordinates by weighted least squares, on coordinates centered and scaled for conditioning (the raw UTM values, near 10⁶, would make the normal matrix singular in float64 from degree 2). It is the simplest spatial model and a step of the GNSS tie.

:func:deramp removes a polynomial ramp from every epoch of a stack. In the temporal-denoising study, per-epoch spatial de-ramping cut velocity scatter by 10 to 24%, more than any temporal filter. The fit only needs sums over pixels (the normal equations), so it stays lazy and exact; it runs on every stride-th pixel by default and the ramp is subtracted from the full grid chunk by chunk. Ramps are long-wavelength by construction: remove them only when long-wavelength signal (plate motion, large subsidence bowls) is not wanted, or mask such areas out of the fit.

Trend ¶

Bases: BaseEstimator

A 2D polynomial trend of degree degree, by weighted least squares.

Parameters:

Name	Type	Description	Default
`degree`	`int`	0 is a constant, 1 a plane, 2 a quadratic surface.	`1`

Attributes:

Name	Type	Description
`coef_`	`ndarray`	Coefficients of the terms in :func:`polynomial_terms` order, for the normalized coordinates `(e - center) / scale`.
`center_, scale_`	`tuple of float`	The normalization.

Examples:

>>> plane = Trend(1).fit(coords, data)
>>> ramp = plane.grid(like=ds)

init ¶

__init__(degree: int = 1) -> None

Store the degree.

deramp ¶

deramp(ds: Dataset, degree: int = 1, *, layer: str = 'displacement', mask: Optional[DataArray] = None, weights: Optional[DataArray] = None, stride: int = 4) -> xr.Dataset

Remove a polynomial ramp from every epoch of layer, lazily.

Parameters:

Name	Type	Description	Default
`ds`	`Dataset`	Data with `y`, `x` (and `time`).	required
`degree`	`int`	Degree of the ramp: 1 a plane (default), 2 a quadratic surface.	`1`
`layer`	`str`	Layer to de-ramp.	`'displacement'`
`mask`	`DataArray`	Boolean `(y, x)` or `(time, y, x)`: pixels to fit the ramp on (for example `disp_xr.valid_mask(ds)`, with deforming areas removed). The ramp is subtracted everywhere.	`None`
`weights`	`DataArray`	Non-negative weights for the fit, for example quality weights.	`None`
`stride`	`int`	Fit on every `stride`-th pixel in `y` and `x` (1 uses all).	`4`

Returns:

Type	Description
`Dataset`	`ds` with `layer` de-ramped and `<layer>_ramp_coefficients` `(time, term)`: the coefficients of each epoch for coordinates normalized as `(e - center) / scale` (`ramp_center`, `ramp_scale` attributes).

Examples:

>>> flat = dxr.models.deramp(ds, degree=1)

polynomial_terms ¶

polynomial_terms(degree: int) -> list[tuple[int, int]]

Return the (i, j) exponents of e**i * n**j with i + j <= degree.

Examples:

>>> polynomial_terms(1)
[(0, 0), (1, 0), (0, 1)]

Gridders¶

Gridders for point data (PRD F9, plan task T24.4).

:class:Spline is the biharmonic spline of Sandwell (1987, https://doi.org/10.1029/GL014i002p00139): the field is a sum of Green's functions g(r) = r² (ln r - 1) centered on the data points, with forces from weighted least squares and optional damping (verde's Spline is the design reference, not a dependency). It suits smooth fields from scattered points (GNSS velocities, downsampled cells). The system is dense, n x n for n points, so keep n to a few thousand (downsample first).

:class:Nearest copies the value of the nearest point with positive weight: a fallback with no smoothing.

Spline ¶

Bases: BaseEstimator

Biharmonic spline gridder (Sandwell 1987) with optional damping.

Parameters:

Name	Type	Description	Default
`damping`	`float`	Regularization of the forces (`damping * \|\|f\|\|²`); larger values give a smoother surface that no longer passes through the data. None or 0 interpolates.	`None`

Attributes:

Name	Type	Description
`force_coords_`	`tuple of numpy.ndarray`	The data points, where the forces sit.
`forces_`	`ndarray`	The force at each point.

Examples:

>>> field = Spline(damping=1e-2).fit(coords, data).grid(like=ds)

init ¶

__init__(damping: Optional[float] = None) -> None

Store the damping.

Nearest ¶

Bases: BaseEstimator

Nearest-neighbor gridder: the value of the closest data point.

Examples:

>>> filled = Nearest().fit(coords, data).grid(like=ds)

init ¶

__init__() -> None

Take no parameters.

GNSS tie¶

GNSS tie: remove long-wavelength DISP error with GNSS (PRD F9, plan task T33).

The residual DISP minus GNSS at the stations is fitted with a smooth surface (a polynomial trend, optionally plus a damped spline) and the predicted field is subtracted from DISP, either for the velocity or for every epoch.

What the tie field may hold. Only long-wavelength error: orbits, long-wavelength troposphere and ionosphere, reference-frame and plate-motion differences. Local deformation (groundwater subsidence, coseismic steps) must stay in the calibrated product, so stations in deforming areas are left out of the fit (exclude) and their signal is restored by construction: the smooth field does not reach it. A plane (degree=1) is the default; add the spline only with damping strong enough that its scale stays well above the station spacing. Leave-one-station-out scores improve as the field gets more local, partly by construction (a field that bends to each station predicts its neighbors well), so a better score alone does not show that the signal was preserved.

fit_tie ¶

fit_tie(comparison: Any, *, mode: Mode = 'velocity', degree: int = 1, damping: Optional[float] = None, exclude: Any = None) -> TieResult

Fit a GNSS tie to the DISP-minus-GNSS residual at the stations.

Parameters:

Name	Type	Description	Default
`comparison`	`Comparison`	:func:`disp_xr.residuals.compare` output.	required
`mode`	`('velocity', 'epoch')`	Fit the velocity residual (`insar_velocity - gps_velocity`, weights `1 / velocity_sigma²`) or the displacement residual of every epoch (weights `1 / sigma_gnss²`).	`"velocity"`
`degree`	`int`	:class:`GNSSTie` parameters.	`1`
`damping`	`int`	:class:`GNSSTie` parameters.	`1`
`exclude`	`list of station ids, or a lon/lat polygon`	Stations left out of the fit, for example those in deforming areas (subsidence, coseismic): their signal must stay in the product.	`None`

Returns:

Type	Description
`TieResult`

Examples:

>>> comparison = dxr.residuals.compare(ds, gps, static)
>>> tie = fit_tie(comparison, mode="velocity")
>>> tie = fit_tie(comparison, exclude=["P539"])

apply_tie ¶

apply_tie(ds: Dataset, tie: TieResult, *, layer: str = 'displacement') -> xr.Dataset

Subtract a GNSS tie from layer.

Epoch mode subtracts each epoch's correction (zero where no tie was fitted). Velocity mode subtracts the correction from a velocity layer, or from displacement as the corrected rate times the time since the first epoch.

Examples:

>>> corrected = apply_tie(ds, tie)

cross_validate_tie ¶

cross_validate_tie(comparison: Any, *, mode: Mode = 'velocity', degree: int = 1, damping: Optional[float] = None, exclude: Any = None) -> pd.DataFrame

Leave-one-station-out: residual RMSE at each station before and after.

Each station used in the fit is left out in turn, the tie is fitted on the others and predicted at it. Excluded stations are scored with the tie of all used stations (they never enter the fit).

Returns:

Type	Description
`DataFrame`	Per station: `rmse_before` (of its residual) and `rmse_after` (after subtracting the prediction), `excluded`; `attrs` holds the network RMSE before and after (`rmse_`), and the median of the per-station values (`median_`). The RMSE is driven by single stations (local subsidence, a short record); the median describes the typical station.

Examples:

>>> scores = cross_validate_tie(comparison, mode="epoch")
>>> scores.attrs["rmse_before"], scores.attrs["rmse_after"]
>>> scores.attrs["median_before"], scores.attrs["median_after"]

GNSSTie ¶

Bases: BaseEstimator

A polynomial trend, optionally plus a damped spline of what it leaves.

Parameters:

Name	Type	Description	Default
`degree`	`int`	Degree of the trend (1, a plane, by default).	`1`
`damping`	`float`	Damping of a biharmonic spline fitted to the trend's residual; None (default) fits the trend alone. Keep it large: the tie field should hold long-wavelength error only.	`None`

Examples:

>>> tie = GNSSTie(degree=1, damping=None)
>>> tie.fit(coords, residuals).grid(like=ds)

init ¶

__init__(degree: int = 1, damping: Optional[float] = None) -> None

Store the parameters.

TieResult `dataclass` ¶

A fitted GNSS tie: one estimator (velocity) or one per epoch.

Examples:

>>> tie = fit_tie(comparison, mode="epoch")
>>> tie.correction(ds).isel(time=-1).plot()
>>> tie.stations

estimators `instance-attribute` ¶

estimators: Any

The fitted tie (velocity mode) or {time: tie} (epoch mode).

stations `instance-attribute` ¶

stations: DataFrame

Stations used and left out, with coordinates and residuals.

correction ¶

correction(like: Union[Dataset, DataArray]) -> xr.DataArray

Return the correction field on the grid of like (lazy in time).

Helpers¶

coordinates ¶

coordinates(obj: Union[Dataset, DataArray, DataFrame]) -> Coords

Return (easting, northing) of a raster, a point table or a DataFrame.

A raster with x and y dimensions gives 2D arrays ((y, x) order); a point table (x and y coordinates along one dimension, as from :func:disp_xr.to_points) or a DataFrame with x and y columns gives 1D arrays.

Examples:

>>> coords = coordinates(ds)

cross_val_score ¶

cross_val_score(estimator: BaseEstimator, coords: Coords, data: Any, *, spacing: float, weights: Optional[Any] = None, n_splits: int = 5, seed: int = 0) -> np.ndarray

Score an estimator with spatial block k-fold cross-validation.

Points are grouped into square blocks of side spacing; the blocks are shuffled and split into n_splits folds. Each fold is predicted by a copy fitted on the other folds, and scored with the weighted R².

Returns:

Type	Description
`ndarray`	One score per fold.

Examples:

>>> from disp_xr.models import Trend
>>> scores = cross_val_score(Trend(1), coords, data, spacing=5e3)

clone ¶

clone(estimator: BaseEstimator) -> BaseEstimator

Return an unfitted copy of estimator with the same parameters.

Examples:

>>> from disp_xr.models import Trend
>>> fresh = clone(Trend(2))
>>> fresh.degree
2

NotFittedError ¶

Bases: RuntimeError

Raised when predict runs before fit.

Spatial models¶

BaseEstimator ¶

get_params ¶

set_params ¶

fit ¶

predict ¶

score ¶

grid ¶

Chain ¶

__init__ ¶

Mean ¶

__init__ ¶

Trends and de-ramping¶

Trend ¶

__init__ ¶

deramp ¶

polynomial_terms ¶

Gridders¶

Spline ¶

__init__ ¶

Nearest ¶

__init__ ¶

GNSS tie¶

fit_tie ¶

apply_tie ¶

cross_validate_tie ¶

GNSSTie ¶

__init__ ¶

TieResult dataclass ¶

estimators instance-attribute ¶

stations instance-attribute ¶

correction ¶

Helpers¶

coordinates ¶

cross_val_score ¶

clone ¶

NotFittedError ¶

init ¶

init ¶

init ¶

init ¶

init ¶

init ¶

TieResult `dataclass` ¶

estimators `instance-attribute` ¶

stations `instance-attribute` ¶