Process

Corrections and masks apply per granule, before chaining; re-referencing comes after. See Concepts.

Corrections

correct

correct(ds: Dataset, which: Union[str, Sequence[str], None] = 'product', *, external: Optional[Mapping[str, DataArray]] = None, external_kind: Literal['differential', 'absolute'] = 'differential') -> xr.Dataset

Subtract correction layers from displacement.

Parameters:

Name	Type	Description	Default
ds	Dataset	Unchained DISP data opened with open_disp(..., corrections=True).	required
which	"product", str or sequence of str	Product layers to subtract: "ionospheric_delay" (alias "ionosphere") and "solid_earth_tide" (alias "set"). "product" (default) subtracts every layer the product carries that its processor has not applied already: both for DISP-S1, the solid-earth tide only for DISP-NISAR. None or empty applies only external.	'product'
external	mapping of str to xr.DataArray	Extra corrections in meters of LOS displacement, positive toward the satellite like displacement, for example a troposphere layer: (time, y, x) or static (y, x). Put on the DISP grid with :func:disp_xr.geometry.match_grid; dates must match exactly.	None
external_kind	('differential', 'absolute')	"differential": each time value is already relative to that epoch's reference_time. "absolute": one value per acquisition date; the correction is value(time) - value(reference_time), so the reference dates must be present too.	"differential"

Returns:

Type	Description
Dataset	displacement corrected; short_wavelength_displacement and the correction layers are left as they are. attrs["corrections_applied"] lists every correction applied so far (JSON).

Raises:

Type	Description
CorrectionError	If ds is chained, a layer is missing or unknown, a correction is already applied, or the processor already applied it (DISP-NISAR ionosphere).

Examples:

>>> raw = dxr.open_disp("F11116/", corrections=True)
>>> ds = dxr.correct(raw)
>>> ds = dxr.correct(raw, ["solid_earth_tide"])
>>> ds = dxr.correct(raw, None, external={"troposphere": tropo})

remove_plate_motion

remove_plate_motion(ds: Dataset, static: Dataset, plate: str) -> xr.Dataset

Subtract rigid plate motion from displacement.

Each epoch loses v_los * (time - reference), where the reference is the granule's reference_time before :func:disp_xr.chain and attrs["reference_date"] after it, so it works on either. The result is in a plate-fixed frame, comparable with plate-fixed GNSS.

Raises:

Type	Description
CorrectionError	If plate motion was already removed.

Examples:

>>> static = dxr.open_static(ds)
>>> fixed = dxr.remove_plate_motion(ds, static, "NOAM")

Quality

mask

mask(ds: Dataset, preset: Optional[str] = 'recommended', **criteria: Any) -> xr.Dataset

Set displacement to NaN where a pixel fails the quality criteria.

Only displacement and short_wavelength_displacement are masked; quality layers are kept as they are. See :func:valid_mask for presets and criteria. Every layer is per epoch, so each epoch is masked by its own granule's quality. Works before or after chaining; the zero epoch that :func:disp_xr.chain prepends has no quality layers and stays 0.

Examples:

>>> masked = dxr.mask(ds, "strict", temporal_coherence=0.7)

valid_mask

valid_mask(ds: Dataset, preset: Optional[str] = 'recommended', *, recommended: Optional[bool] = None, temporal_coherence: Optional[float] = None, phase_similarity: Optional[float] = None, connected_components: Optional[bool] = None, water: Optional[bool] = None) -> xr.DataArray

Return True where a pixel passes every chosen quality criterion.

Presets: "recommended" uses the product's recommended_mask; "strict" adds connected components, water and temporal_coherence >= the threshold stored on recommended_mask (0.6 if absent). Keyword arguments add to or override the preset; False turns a criterion off. A pixel at an integer layer's no-data value fails.

Raises:

Type	Description
ValueError	If no criterion is chosen, or a needed layer is not in ds.

Examples:

>>> good = dxr.valid_mask(ds, "recommended")
>>> good = dxr.valid_mask(ds, None, temporal_coherence=0.7)

quality_summary

quality_summary(ds: Dataset, *, reliability_threshold: float = 0.9, by_ministack: bool = False) -> xr.Dataset

Summarize quality per pixel over time (lazy).

Percentages divide by the epochs where the layer has data, so a stacked zero epoch or missing granules do not count as bad.

Returns:

Type	Description
Dataset	pct_masked (% of epochs with recommended_mask == 0), pct_ps (% of epochs flagged persistent scatterer), temporal_coherence_hmean (harmonic mean, which down-weights low coherence periods), reliability (fraction of epochs with recommended_mask == 1) and reliable (reliability >= reliability_threshold). With by_ministack, each is computed per reference_time group (flox groupby).

Examples:

>>> summary = dxr.quality_summary(ds).compute()
>>> summary.reliable.plot.imshow()
>>> per_ministack = dxr.quality_summary(ds, by_ministack=True)

quality_weights

quality_weights(ds: Dataset, layers: Sequence[str] = DEFAULT_LAYERS, *, exponents: Optional[Mapping[str, float]] = None, floors: Optional[Mapping[str, float]] = None, penalty: float = 0.5, f0: float = 0.25, min_residual: float = 0.05) -> xr.Dataset

Combine quality layers into pixel weights between 0 and 1.

Parameters:

Name	Type	Description	Default
ds	Dataset	DISP data with the chosen layers (as from :func:disp_xr.open_disp), chained or not.	required
layers	sequence of str	Any of phase_similarity, temporal_coherence, timeseries_inversion_residuals.	DEFAULT_LAYERS
exponents	mapping	Relative weight of each layer in the geometric mean (default 1 each).	None
floors	mapping	Stretch floor per similarity or coherence layer. Default: the processor's thresholds recorded on recommended_mask, else 0.3 for similarity and 0.6 for coherence.	None
penalty	float	Residual weight factor per 2 pi cycle jump.	0.5
f0	float	Scale of the residual weight's decay with the distance to the nearest whole cycle.	0.25
min_residual	float	Lowest residual weight: cycle jumps down-weight, never exclude.	0.05

Returns:

Type	Description
Dataset	weight (time, y, x) and static_weight (y, x), one weight_<layer> per layer, and for the residuals cycles, ambiguous and jump_fraction (share of dates with a jump). Lazy if the input is.

Examples:

>>> w = dxr.quality_weights(ds)
>>> fit = dxr.fit_trajectory(ds, weights=w.weight)

weights

Quality weights from the DISP quality layers (PRD S2, plan task T26).

Each layer maps to a weight between 0 and 1, and the weights combine as a weighted geometric mean, prod(w_i ** (a_i / sum(a))):

phase_similarity, temporal_coherence (good near 1) Linear stretch clip((v - floor) / (1 - floor), 0, 1): values at or below the floor get weight 0 and exclude the pixel at that epoch. timeseries_inversion_residuals (good near 0) The per-date sum of absolute residuals, in radians, of dolphin's L1 time-series inversion over every interferogram that includes the date. An unwrapping error of one cycle in one interferogram adds about 2 pi to both of its dates, so on real granules the residuals cluster near multiples of 2 pi (residual_cycle_histogram). The residual splits into cycles = round(r / 2 pi) and the distance to that lattice frac = |r / 2 pi - cycles| (0 to 0.5). The weight is penalty ** cycles * exp(-frac / f0), clamped to min_residual so a cycle jump down-weights a pixel but never removes it (PRD decision). A frac near 0.5 is ambiguous between two cycle counts and is flagged.

The per-epoch weight (time, y, x) combines the per-epoch layer weights; the static weight (y, x) combines their time means, with the residual term replaced by one minus the share of dates that have a cycle jump.

Defaults for the floors, penalty, f0 and min_residual are provisional until the trade study (plan task T29).

residual_cycles

residual_cycles(residuals: DataArray) -> xr.Dataset

Split inversion residuals into 2 pi cycles and the distance to them.

Returns:

Type	Description
Dataset	cycles (integer count of 2 pi), frac (0 to 0.5, distance to the nearest whole cycle) and ambiguous (frac within AMBIGUOUS_WIDTH of 0.5: the cycle count could be either).

Examples:

>>> r = xr.DataArray([0.1, 2 * np.pi + 0.2, 3 * np.pi], dims="time")
>>> residual_cycles(r).cycles.values.tolist()
[0, 1, 2]

network_degree

network_degree(ds: Dataset, max_bandwidth: int) -> xr.DataArray

Return the interferogram count of each epoch's date in its ministack.

In a bandwidth-limited network each date pairs with up to max_bandwidth neighbours on each side within its ministack, whose first date is the (compressed) reference. A secondary at position p of N dates has min(p, b) + min(N - 1 - p, b) interferograms.

Examples:

>>> degree = network_degree(ds, max_bandwidth=3)

max_bandwidth_from_granule

max_bandwidth_from_granule(path: Union[str, Path]) -> Optional[int]

Read interferogram_network.max_bandwidth from a granule's dolphin config.

Returns None if the config is missing or the network is not bandwidth-limited.

residual_cycle_histogram

residual_cycle_histogram(ds: Dataset, bins: Optional[ndarray] = None, min_residual: float = 0.05) -> pd.Series

Histogram of the inversion residuals in units of 2 pi.

On real DISP-S1 granules the counts peak near whole numbers (one 2 pi unwrapping error per interferogram), which the cycle decomposition of :func:quality_weights assumes. Residuals below min_residual radians (clean pixels) are left out so the peaks show.

Returns:

Type	Description
Series	Counts indexed by the lower bin edge, in cycles of 2 pi.

Examples:

>>> counts = residual_cycle_histogram(ds)

Chaining

chain

chain(ds: Dataset, *, start: Optional[Union[str, datetime, Timestamp]] = None, end: Optional[Union[str, datetime, Timestamp]] = None, nan_policy: Union[str, NaNPolicy] = NaNPolicy.propagate, zero_epoch: bool = True) -> xr.Dataset

Chain ministacks so displacement is relative to one reference date.

Parameters:

Name	Type	Description	Default
ds	Dataset	Output of :func:disp_xr.open_disp: time and reference_time coordinates, displacement relative to each granule's reference date.	required
start	date - like	Chain from this ministack reference date; earlier ministacks are dropped. Must be one of the reference dates.	None
end	date - like	Drop epochs after this date.	None
nan_policy	('propagate', 'omit')	See :class:NaNPolicy.	"propagate"
zero_epoch	bool	Prepend the reference date as an epoch with zero displacement. Other layers get their no-data value there (NaN for floats, _FillValue for integers), so their dtype is kept.	True

Returns:

Type	Description
Dataset	Same variables; displacement and short_wavelength_displacement relative to the first reference date, with the whole time axis in one chunk. attrs["reference_date"] records that date.

Raises:

Type	Description
MissingCrossoverError	If a reference date after the first is not an epoch of ds.
ValueError	If start is not a reference date.

Examples:

>>> ds = dxr.chain(dxr.open_disp("F11116/"))
>>> ds = dxr.chain(raw, start="2020-01-01")

NaNPolicy

Bases: str, Enum

How a NaN at a crossover epoch affects later ministacks.

propagate (default): later epochs of that pixel become NaN, since their offset is unknown. omit: the crossover's own displacement is taken as zero, which assumes no motion over the missing interval, so later epochs carry the offset accumulated before it.

Examples:

>>> NaNPolicy("omit")
<NaNPolicy.omit: 'omit'>

Spatial reference

rereference

rereference(ds: Dataset, method: Union[str, ReferenceMethod] = 'auto', *, row: Optional[int] = None, col: Optional[int] = None, lon: Optional[float] = None, lat: Optional[float] = None, bbox: Optional[Sequence[float]] = None, geometry: Optional[Any] = None, border_pixels: int = 3, good_pixels: Optional[DataArray] = None, check_point: bool = True) -> xr.Dataset

Subtract a per-epoch reference value from the displacement layers.

Parameters:

Name	Type	Description	Default
ds	Dataset	DISP dataset, usually after :func:disp_xr.chain.	required
method	str or ReferenceMethod	"point", "median", "border", "high_coherence" (needs good_pixels), "area" or "auto" (:func:find_stable_point).	'auto'
row	int	Pixel for point.	None
col	int	Pixel for point.	None
lon	float	Point in degrees for point, used when row/col are not given.	None
lat	float	Point in degrees for point, used when row/col are not given.	None
bbox	sequence of float	(west, south, east, north) in degrees for area.	None
geometry	shapely geometry	Polygon in lon/lat for area. area needs shapely (aws extra).	None
border_pixels	int	Border width for border.	3
good_pixels	xr.DataArray of bool	Pixels the medians may use, for example disp_xr.valid_mask(ds, "strict").	None
check_point	bool	For point and auto, read the reference pixel's displacement and raise if it is NaN at every epoch.	True

Returns:

Type	Description
Dataset	displacement and short_wavelength_displacement re-referenced; attrs["reference"] records the method and, for a point, its row, column and map coordinates.

Raises:

Type	Description
OutOfBoundsError	If the point is outside the grid or masked, or the area is empty.

Examples:

>>> ds = dxr.rereference(ds, "median")
>>> ds = dxr.rereference(ds, "point", lon=-120.43, lat=35.90)
>>> ds = dxr.rereference(
...     ds, "area", bbox=(-120.5, 35.8, -120.4, 35.9)
... )

Time series

fit_trajectory

fit_trajectory(obj: Union[Dataset, DataArray], terms: Optional[Trajectory] = None, *, layer: str = 'displacement', weights: Optional[DataArray] = None, min_valid: int = 0, allow_denoised: bool = False) -> xr.Dataset

Fit a trajectory model to every pixel by ordinary least squares.

Parameters:

Name	Type	Description	Default
obj	Dataset or DataArray	Chained displacement (time, y, x) (see :func:disp_xr.chain), or a Dataset holding it as layer. NaN epochs are skipped per pixel.	required
terms	Trajectory	Model terms; default offset, velocity and an annual cycle.	None
layer	str	Layer of a Dataset to fit.	'displacement'
weights	DataArray	Per-epoch weights (time, y, x) or (time,), proportional to 1 / variance (for example from quality layers). Zero drops an epoch.	None
min_valid	int	Fewest valid epochs to fit a pixel; at least one more than the number of coefficients is always required.	0
allow_denoised	bool	Fit data marked attrs["denoised"]. Off by default: fits use the raw series (PRD S1).	False

Returns:

Type	Description
Dataset	coefficients and sigma (coefficient, y, x) in the input unit (per year for rates), velocity and velocity_sigma when the model has a linear term, <period>_amplitude and <period>_amplitude_sigma for each sine/cosine pair, rms of the residuals and n_valid epochs. Lazy if the input is.

Examples:

>>> terms = Trajectory(steps=("2023-02-06",))
>>> fit = dxr.fit_trajectory(ds, terms)
>>> fit.velocity.disp.plot_map()

Trajectory dataclass

Terms of a trajectory model.

Parameters:

Name	Type	Description	Default
polynomial	int	Degree: 0 offset only, 1 adds velocity, 2 adds acceleration (the quadratic coefficient is a / 2-normalized, so it is the acceleration).	1
periods	tuple of float	Periods in years of sine/cosine pairs: 1.0 annual, 0.5 semiannual.	(1.0,)
steps	tuple of str	Dates of Heaviside steps (any pandas date string).	()
postseismic_log	tuple of (str, float)	(onset date, tau in days) of log1p(dt / tau) and 1 - exp(-dt / tau) terms.	()
postseismic_exp	tuple of (str, float)	(onset date, tau in days) of log1p(dt / tau) and 1 - exp(-dt / tau) terms.	()

Examples:

>>> Trajectory(steps=("2023-02-06",)).names
['constant', 'linear_trend', 'annual_sin', 'annual_cos', 'step_2023-02-06']

names property

names: list[str]

Coefficient names, in column order of the design matrix.

design

design(times: ndarray) -> np.ndarray

Return the (time, coefficient) design matrix, time 0 at times[0].

predict_trajectory

predict_trajectory(fit: Dataset, times: Union[ndarray, Sequence], dtype: str = 'float32') -> xr.DataArray

Evaluate a fitted trajectory at times (lazily if the fit is).

The model is rebuilt from the coefficient names and the reference time in fit.attrs, so a fit read back from a store can be evaluated. Computed in dtype and never stored unless you write it.

Examples:

>>> fit = dxr.fit_trajectory(ds)
>>> model = dxr.predict_trajectory(fit, ds.time.values)
>>> residual = ds.displacement - model

Velocity for every pixel

Velocity for every pixel: OLS, Whittle colored noise, or MIDAS.

fit_velocity(ds, method=...) runs per pixel on a full frame:

"ols" The linear term of :func:~disp_xr.timeseries.trajectory.fit_trajectory, with white-noise sigmas (too small for real InSAR series). "whittle" The same OLS velocity, with noise estimated per pixel from the periodogram of the OLS residuals (Whittle likelihood, power law plus white): kappa, the white fraction and amplitudes, and the colored-noise sigma of the velocity. Reliable from about 200 epochs; on short records it underestimates |kappa| and the sigma (benchmarks/whittle_bias.py). The sigma uses the covariance of the complete record, so it is approximate for pixels with gaps, and the amplitudes come from the gap-filled grid, which interpolation smooths (about 10% low with 20% gaps). The exact grouped GLS is plan task T23. "midas" The median of slopes between epoch pairs about one year apart, trimmed at two scaled MADs (Blewitt et al., 2016): robust to outliers, seasonal signals and steps, with its own sigma. A step biases it when the pairs that span it are a large share, about one year over the record length, so prefer records of several years. Pairs are chosen from the stack dates; a pair with a NaN end is dropped, where per-series MIDAS would choose a replacement.

Blewitt, G., Kreemer, C., Hammond, W. C., & Gazeaux, J. (2016). MIDAS robust trend estimator for accurate GPS station velocities without step detection. JGR Solid Earth, 121, 2054-2068. https://doi.org/10.1002/2015JB012552

fit_velocity

fit_velocity(obj: Union[Dataset, DataArray], method: Literal['ols', 'whittle', 'midas', 'gls'] = 'whittle', *, terms: Optional[Trajectory] = None, layer: str = 'displacement', min_valid: int = 0, noise: str = 'PLWN', groups: Union[str, DataArray] = 'coherence', n_groups: int = 5, n_samples: int = 16, allow_denoised: bool = False) -> xr.Dataset

Estimate a velocity and its sigma for every pixel.

Parameters:

Name	Type	Description	Default
obj	Dataset or DataArray	Chained displacement (time, y, x), or a Dataset holding it.	required
method	('whittle', 'gls', 'ols', 'midas')	"whittle" (default): colored-noise sigmas, reliable from about 200 epochs (6.5 years at 12 days); a warning is logged below. "gls": exact generalized least squares, with the noise shape estimated per group of pixels by exact maximum likelihood and the amplitude per pixel; works on short records too, costs more. "midas": robust to unmodelled steps and outliers, for records of several years. "ols": fast, but its sigma ignores correlated noise and is several times too small. For shorter records use :func:disp_xr.fit_noise on selected pixels. See the guide "Choosing a velocity method" in the documentation.	"whittle"
terms	Trajectory	Deterministic model for "ols" and "whittle" (default offset, velocity and annual cycle). MIDAS needs none.	None
layer	str	Layer of a Dataset.	'displacement'
min_valid	int	Fewest valid epochs to fit a pixel ("ols", "whittle", "gls").	0
noise	('PLWN', 'FNWN', 'RWFNWN', 'WN', 'auto')	Noise model of "gls"; "auto" picks FNWN, PLWN or RWFNWN per group by the median BIC of its sampled pixels.	"PLWN"
groups	'coherence' or DataArray	"gls" groups: n_groups quantile bins of time-mean temporal_coherence (needs the Dataset), or an integer label map (y, x); -1 leaves a pixel out.	'coherence'
n_groups	int	Coherence bins, and pixels per group fitted to estimate its noise.	5
n_samples	int	Coherence bins, and pixels per group fitted to estimate its noise.	5
allow_denoised	bool	Fit data marked attrs["denoised"]. Off by default: velocities and their sigmas come from the raw series (PRD S1).	False

Returns:

Type	Description
Dataset	velocity and velocity_sigma (unit per year) for every method; kappa, sigma_powerlaw, sigma_white and n_valid for Whittle; outlier_fraction and velocity_scatter for MIDAS; noise_amplitude, noise_group, kappa, white_fraction and n_valid for GLS, with the group table in attrs["noise_groups"] (JSON). Lazy if the input is, except the sampled pixels GLS reads to estimate the noise.

Examples:

>>> vel = dxr.fit_velocity(ds, "midas")

gls

Full-frame velocity by generalized least squares with colored noise (T23).

The noise covariance of a pixel depends on the epochs, which of them are valid, and the noise parameters, not on the pixel's values. So the covariance shape (the spectral index and the mix of noise components) is estimated per group of pixels, by exact maximum likelihood (:func:~disp_xr.timeseries.noise.fit_noise) on a sample of each group's pixels, and every pixel is then fitted by exact GLS under its group's shape, with its own noise amplitude profiled from its residuals. Each pixel gets a colored-noise velocity sigma without an n x n factorization per pixel.

Groups are quantile bins of time-averaged temporal_coherence by default, or any integer label map (for example quadtree cells, plan task T28).

Missing epochs. Each group's covariance over the complete record is inverted once, Q = C⁻¹. For a pixel missing the epochs m, the inverse covariance of its valid epochs is the Schur complement P = Q - Q[:, m] Q[m, m]⁻¹ Q[m, :], which is zero on the missing rows and columns. The normal equations and residual quadratic forms follow from Q, Q A and a k x k inverse per pixel (k = len(m)), batched over pixels with the same k.

NoiseShape dataclass

Covariance shape of a group: components (kappa, weight), weights sum to 1.

The covariance of a pixel is sigma2 * sum(w * C(kappa)) over the components, with C(0) the identity and sigma2 the pixel's own amplitude.

kappa property

kappa: float

Spectral index of the power-law component (NaN for mixtures).

white_fraction property

white_fraction: float

Share of the variance in the white component.

covariance

covariance(obs_idx: ndarray) -> np.ndarray

Return the unit-amplitude covariance at the grid indices obs_idx.

coherence_groups

coherence_groups(ds: Dataset, n_groups: int = 5) -> xr.DataArray

Label pixels by quantile bin of their (time-mean) temporal_coherence.

estimate_group_noise

estimate_group_noise(da: DataArray, labels: DataArray, *, terms: Trajectory, noise: str = 'PLWN', n_samples: int = 16, seed: int = 0, max_workers: Optional[int] = None) -> dict[int, NoiseShape]

Estimate a noise shape per group from a random sample of its pixels.

Each sampled pixel is fitted with :func:fit_noise (exact RMLE); the group shape is the median of the samples' shapes. With noise="auto", FNWN, PLWN and RWFNWN are fitted and the model with the lowest median BIC is used for the group.

Model diagnostics

Per-pixel diagnostics of a fitted trajectory model (PRD S4, plan T32.6b-c).

:func:model_diagnostics describes the residuals of every pixel's fit: mean, standard deviation, skewness, excess kurtosis, autocorrelation at chosen lags and the Anderson-Darling statistic for normality. White-noise sigmas (the sigma of :func:disp_xr.fit_trajectory) hold only when the residuals are uncorrelated and Gaussian; autocorrelation above about 0.3 or an Anderson-Darling statistic above about 0.75 (the 5% critical value) says they are optimistic, and :func:disp_xr.fit_velocity with a colored-noise method is the better tool. The model is evaluated on the fly from the coefficients (:func:disp_xr.predict_trajectory), never stored, and every statistic comes from one pass per block of pixels over time-contiguous chunks.

:func:significance turns coefficients and sigmas into z-scores, two-sided p-values and confidence intervals.

The statistics follow the Brisbane notebook (branch _stats_analysis.py), with Anderson-Darling vectorized over pixels in NumPy instead of a per-pixel apply_ufunc(vectorize=True).

model_diagnostics

model_diagnostics(data: Union[Dataset, DataArray], fit: Dataset, *, layer: str = 'displacement', lags: Sequence[int] = (1,)) -> xr.Dataset

Describe the residuals of a trajectory fit at every pixel, lazily.

Parameters:

Name	Type	Description	Default
data	Dataset or DataArray	The data that were fitted, (time, y, x).	required
fit	Dataset	:func:disp_xr.fit_trajectory output for data.	required
layer	str	Layer of a Dataset.	'displacement'
lags	sequence of int	Lags, in epochs, of the residual autocorrelation.	(1,)

Returns:

Type	Description
Dataset	(y, x): n (residuals used), residual_mean, residual_std, skewness, kurtosis (excess), autocorrelation_lag<k> and anderson_darling (normality; above AD_CRITICAL_5 rejects at 5%).

Examples:

>>> fit = dxr.fit_trajectory(ds)
>>> diag = model_diagnostics(ds, fit, lags=(1, 2))
>>> (diag.autocorrelation_lag1 > 0.3).mean()

significance

significance(fit: Dataset, alpha: float = 0.05) -> xr.Dataset

Return z-scores, p-values and confidence intervals of fitted coefficients.

Parameters:

Name	Type	Description	Default
fit	Dataset	:func:disp_xr.fit_trajectory output (coefficients and sigma, (coefficient, y, x)).	required
alpha	float	Confidence intervals cover 1 - alpha.	0.05

Returns:

Type	Description
Dataset	z (coefficient over sigma), p_value (two-sided, Gaussian), ci_low and ci_high, and significant (p_value < alpha).

Examples:

>>> sig = significance(dxr.fit_trajectory(ds))
>>> sig.sel(coefficient="linear_trend").significant.mean()

anderson_darling

anderson_darling(values: ndarray, axis: int = -1) -> np.ndarray

Return the Anderson-Darling normality statistic along axis, ignoring NaN.

Mean and variance are estimated from the data, with Stephens' small-sample correction A² (1 + 0.75 / n + 2.25 / n²). Vectorized: one sort for all series. NaN where fewer than 8 values are finite.

Examples:

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> bool(anderson_darling(rng.standard_normal(500)) < AD_CRITICAL_5)
True

Time grids, gaps and breaks

Regular time grids, refused gaps and declared breaks (PRD S3, plan task T30).

Rules from the temporal-denoising study (Kiritimati DISP-S1, 2026-09):

• Grid step. 15 d where sampling is genuinely 12-day, 30 d for 24-day sampling. A 15-day grid is finer than 24-day sampling supports. The step is chosen from the 75th percentile of the acquisition intervals rather than the median: the study's pre segment has a 12-day median (18 of its 34 intervals) but 14 intervals of 24 days, and the study recommends 30 d for such frames. The step is 15 d per 12 d of that percentile, rounded.
• Gaps are refused, not interpolated. A grid point farther than max_gap_days / 2 from every acquisition lies inside a hole longer than max_gap_days and is marked; a 15-day grid across pre's 288-day hole would otherwise invent about 19 samples.
• Breaks are declared (earthquake catalogs, user input); no automatic detection before 1.0. A break at tb lands on the first epoch at or after it: the first observation that already contains the jump. Filters that know nothing about discontinuities respect them when run per segment (:func:apply_segmented).

regular_time_grid

regular_time_grid(times: TimesLike, step: Optional[int] = None, *, start: Optional[Any] = None, end: Optional[Any] = None) -> pd.DatetimeIndex

Return a regular grid of dates covering the acquisitions.

Parameters:

Name	Type	Description	Default
times	array of datetimes	Acquisition times, increasing.	required
step	int	Grid step in days; default :func:grid_step_days.	None
start	date - like	Grid bounds; default the first and last acquisition days.	None
end	date - like	Grid bounds; default the first and last acquisition days.	None

Returns:

Type	Description
DatetimeIndex	Midnight dates from start every step days, not beyond end.

Examples:

>>> import pandas as pd
>>> times = pd.date_range("2020-01-01", periods=30, freq="12D")
>>> grid = regular_time_grid(times)
>>> int((grid[1] - grid[0]).days)
15

grid_step_days

grid_step_days(times: TimesLike) -> int

Return the grid step in days for these acquisitions (15 per 12 d of sampling).

Uses the 75th percentile of the intervals, so a frame with a mix of 12- and 24-day intervals gets 30 d (see the module docstring).

Examples:

>>> import pandas as pd
>>> grid_step_days(pd.date_range("2020-01-01", periods=30, freq="12D"))
15

gap_mask

gap_mask(times: TimesLike, grid: TimesLike, max_gap_days: Optional[float] = None) -> np.ndarray

Return True for grid dates inside a refused gap.

A grid date farther than max_gap_days / 2 from every acquisition lies in a hole longer than max_gap_days; resampled values there must be NaN, never interpolated.

Parameters:

Name	Type	Description	Default
times	array of datetimes	Acquisition times.	required
grid	array of datetimes	Grid dates (for example :func:regular_time_grid).	required
max_gap_days	float	Longest hole to bridge; default four times the typical (75th percentile) interval, so a missed acquisition or two is bridged and a season-long hole is not.	None

Returns:

Type	Description
ndarray	Boolean, one per grid date.

Examples:

>>> import pandas as pd
>>> times = pd.to_datetime(["2020-01-01", "2020-01-13", "2020-09-01"])
>>> grid = pd.date_range("2020-01-01", "2020-09-01", freq="30D")
>>> gap_mask(times, grid, max_gap_days=60).sum()

split_at_breaks

split_at_breaks(times: TimesLike, breaks: Optional[Sequence[Any]]) -> list[np.ndarray]

Return the epoch indices of each segment between breaks.

Examples:

>>> import pandas as pd
>>> times = pd.date_range("2020-01-01", periods=6, freq="12D")
>>> [len(s) for s in split_at_breaks(times, ["2020-02-01"])]
[3, 3]

apply_segmented

apply_segmented(fn: Callable, data: Union[DataArray, Dataset, ndarray], breaks: Optional[Sequence[Any]], *, times: Optional[TimesLike] = None, min_len: int = 3) -> Any

Run fn on each segment between breaks independently.

A filter that knows nothing about discontinuities then never smooths across one, at the cost of one-sided support at each segment edge.

Parameters:

Name	Type	Description	Default
fn	callable	For xarray input, fn(segment) -> segment on the time slice (same shape). For NumPy input, fn(t_segment, y_segment) with y of shape (time, ...).	required
data	(DataArray, Dataset or ndarray)	Series with time first (NumPy) or a time dimension (xarray).	required
breaks	sequence of date-like	Break dates.	required
times	array of datetimes	Acquisition times for NumPy input (xarray uses its time).	None
min_len	int	Segments shorter than this are passed through unchanged: one or two epochs carry nothing a smoother could use.	3

Returns:

Type	Description
Same type as ``data``.

Examples:

>>> breaks = ["2023-02-06"]
>>> smooth = apply_segmented(my_filter, ds.displacement, breaks)

Temporal denoising

Temporal denoising: gated Kalman smoother and robust LOWESS (PRD S3, T31).

Both methods share one interface, :func:denoise, and evaluate on the acquisitions or on a regular grid (:func:disp_xr.timeseries.regular_time_grid) with refused gaps. They are for display, resampling, linearity checks and transient hunting: the temporal-denoising study (Kiritimati DISP-S1, 2026-09) found that no filter improves a velocity, because InSAR residuals at 12-day sampling are white and ordinary least squares is already optimal. Fit rates on the raw series; :func:disp_xr.fit_trajectory refuses data marked as denoised.

Kalman (method="kalman"): a local-linear-trend model (displacement and velocity, an integrated random walk) with a Rauch-Tung-Striebel smoother, or a fixed-lag smoother for products that publish as data arrive. The only knobs are the timescale tau_days (20 to 30 d), which sets the process noise q = 3 R / tau³ against the measurement variance R (estimated per pixel from the data), and the χ² gate: an update whose normalized innovation squared exceeds gate (9) is skipped. A one-epoch outlier makes one huge innovation and is rejected; a real transient builds over several epochs that stay inside the gate. In the study it left 0.2 to 0.4% of an injected spike and kept 93% of a 60-day transient. Declared breaks inflate the covariance at the first epoch after them (break_mode: offset, rate or both), and the gate is suspended there, so a real step is kept; the inflated covariance also stops the backward pass at the break.

Unlike the study, each pixel keeps its own 2 x 2 covariance: DISP pixels have their own missing epochs, and the study's production notes ask for a per-pixel gate (unwrapping jumps) rather than its scene-wide one.

LOWESS (method="lowess"): robust local-linear regression (tricube kernel of bandwidth_days, 60 d, and iterations bisquare reweightings, 2, per pixel), defined at any time. Its outlier rejection depends on sampling density (0.1% of a spike survives at 12-day sampling, about 30% at 24-day). Breaks are handled by fitting each segment separately.

Known limitation: both methods were evaluated on one flat, aseismic site (Kiritimati). A deforming site with topography-correlated troposphere is part of the T29 trade-study frames.

denoise

Temporal denoising: gated Kalman smoother and robust LOWESS (PRD S3, T31).

Both methods share one interface, :func:denoise, and evaluate on the acquisitions or on a regular grid (:func:disp_xr.timeseries.regular_time_grid) with refused gaps. They are for display, resampling, linearity checks and transient hunting: the temporal-denoising study (Kiritimati DISP-S1, 2026-09) found that no filter improves a velocity, because InSAR residuals at 12-day sampling are white and ordinary least squares is already optimal. Fit rates on the raw series; :func:disp_xr.fit_trajectory refuses data marked as denoised.

Kalman (method="kalman"): a local-linear-trend model (displacement and velocity, an integrated random walk) with a Rauch-Tung-Striebel smoother, or a fixed-lag smoother for products that publish as data arrive. The only knobs are the timescale tau_days (20 to 30 d), which sets the process noise q = 3 R / tau³ against the measurement variance R (estimated per pixel from the data), and the χ² gate: an update whose normalized innovation squared exceeds gate (9) is skipped. A one-epoch outlier makes one huge innovation and is rejected; a real transient builds over several epochs that stay inside the gate. In the study it left 0.2 to 0.4% of an injected spike and kept 93% of a 60-day transient. Declared breaks inflate the covariance at the first epoch after them (break_mode: offset, rate or both), and the gate is suspended there, so a real step is kept; the inflated covariance also stops the backward pass at the break.

Unlike the study, each pixel keeps its own 2 x 2 covariance: DISP pixels have their own missing epochs, and the study's production notes ask for a per-pixel gate (unwrapping jumps) rather than its scene-wide one.

LOWESS (method="lowess"): robust local-linear regression (tricube kernel of bandwidth_days, 60 d, and iterations bisquare reweightings, 2, per pixel), defined at any time. Its outlier rejection depends on sampling density (0.1% of a spike survives at 12-day sampling, about 30% at 24-day). Breaks are handled by fitting each segment separately.

Known limitation: both methods were evaluated on one flat, aseismic site (Kiritimati). A deforming site with topography-correlated troposphere is part of the T29 trade-study frames.

q_from_timescale

q_from_timescale(tau_days: float, r: Union[float, ndarray]) -> Any

Return the process noise giving a smoothing timescale of tau_days.

For an integrated random walk the timescale where process noise balances measurement noise r is tau = (3 r / q)^(1/3) (in years). Setting q directly is how a filter silently returns a straight line (q too small for r).

kalman_smooth

kalman_smooth(t_years: ndarray, y: ndarray, *, tau_days: float = 20.0, gate: Optional[float] = 9.0, r: Optional[ndarray] = None, p0_velocity: float = 1.0, breaks: Optional[ndarray] = None, break_mode: Literal['offset', 'rate', 'both'] = 'offset', break_var: float = 100.0, lag: Optional[int] = None) -> tuple[np.ndarray, np.ndarray]

Gated local-linear-trend Kalman filter and smoother on (time, npix).

Parameters:

Name	Type	Description	Default
t_years	ndarray	Increasing times in years; epochs where y is NaN are predicted without an update (use them for grid times).	required
y	ndarray	(time, npix) series.	required
tau_days	float	Smoothing timescale and χ² gate (None turns the gate off).	20.0
gate	float	Smoothing timescale and χ² gate (None turns the gate off).	20.0
r	ndarray	Measurement variance per pixel; default estimated from y.	None
p0_velocity	float	Prior velocity variance; must stay diffuse (a tight prior shrinks every rate toward zero), in data units squared per year squared.	1.0
breaks	ndarray	Boolean per epoch: the covariance is inflated there.	None
break_mode	('offset', 'rate', 'both')	What may jump at a break.	"offset"
break_var	float	Variance added at a break, in data units squared (and per year squared for the rate).	100.0
lag	int	Fixed-lag smoothing: epoch j uses data through j + lag (about 3 tau converges). None is the full RTS smoother.	None

Returns:

Type	Description
(ndarray, ndarray)	Smoothed displacement and velocity, (time, npix).

lowess

lowess(t_years: ndarray, y: ndarray, t_out: ndarray, *, bandwidth_days: float = 60.0, iterations: int = 2, min_points: int = 3) -> np.ndarray

Robust local-linear regression of (time, npix) series at t_out.

The window at each output time is bandwidth_days wide on each side, widened to the min_points nearest acquisitions where sampling is sparse; pixels with fewer than min_points valid values in it are NaN. Each of iterations passes reweights points by the bisquare of their residual over six times the pixel's median absolute residual (Cleveland 1979, https://doi.org/10.1080/01621459.1979.10481038).

denoise

denoise(obj: Union[Dataset, DataArray], method: Literal['kalman', 'lowess'] = 'kalman', *, grid: Union[None, str, Sequence, DatetimeIndex] = None, layer: str = 'displacement', tau_days: float = 20.0, gate: Optional[float] = 9.0, lag: Optional[int] = None, bandwidth_days: float = 60.0, iterations: int = 2, breaks: Optional[Sequence[Any]] = None, break_mode: Literal['offset', 'rate', 'both'] = 'offset', max_gap_days: Optional[float] = None) -> Any

Denoise displacement series in time, on the acquisitions or a grid.

Parameters:

Name	Type	Description	Default
obj	Dataset or DataArray	(time, y, x) (or any dims with time) series.	required
method	('kalman', 'lowess')	See the module docstring.	"kalman"
grid	(None, 'auto' or dates)	Output times: the acquisitions (None), a regular grid ("auto": :func:disp_xr.timeseries.regular_time_grid), or given dates. Grid dates inside gaps longer than max_gap_days are NaN.	None
layer	str	Layer of a Dataset.	'displacement'
tau_days	float	Kalman timescale, χ² gate (None: off) and fixed lag (None: RTS).	20.0
gate	float	Kalman timescale, χ² gate (None: off) and fixed lag (None: RTS).	20.0
lag	float	Kalman timescale, χ² gate (None: off) and fixed lag (None: RTS).	20.0
bandwidth_days	float	LOWESS half-window and robust passes.	60.0
iterations	float	LOWESS half-window and robust passes.	60.0
breaks	sequence of date-like	Declared discontinuities.	None
break_mode	('offset', 'rate', 'both')	Kalman: what may jump at a break.	"offset"
max_gap_days	float	Longest hole a grid may bridge (default from :func:disp_xr.timeseries.gap_mask).	None

Returns:

Type	Description
Same type as ``obj``	Denoised layer (a Dataset keeps only that layer) with time the output times; attrs["denoised"] names the method, so fits refuse it. Lazy if the input is.

Examples:

>>> smooth = dxr.denoise(ds, "kalman", tau_days=20)
>>> on_grid = dxr.denoise(ds, "lowess", grid="auto")
>>> stepped = dxr.denoise(ds, breaks=["2023-02-06"])

kalman_smooth

kalman_smooth(t_years: ndarray, y: ndarray, *, tau_days: float = 20.0, gate: Optional[float] = 9.0, r: Optional[ndarray] = None, p0_velocity: float = 1.0, breaks: Optional[ndarray] = None, break_mode: Literal['offset', 'rate', 'both'] = 'offset', break_var: float = 100.0, lag: Optional[int] = None) -> tuple[np.ndarray, np.ndarray]

Gated local-linear-trend Kalman filter and smoother on (time, npix).

Parameters:

Name	Type	Description	Default
t_years	ndarray	Increasing times in years; epochs where y is NaN are predicted without an update (use them for grid times).	required
y	ndarray	(time, npix) series.	required
tau_days	float	Smoothing timescale and χ² gate (None turns the gate off).	20.0
gate	float	Smoothing timescale and χ² gate (None turns the gate off).	20.0
r	ndarray	Measurement variance per pixel; default estimated from y.	None
p0_velocity	float	Prior velocity variance; must stay diffuse (a tight prior shrinks every rate toward zero), in data units squared per year squared.	1.0
breaks	ndarray	Boolean per epoch: the covariance is inflated there.	None
break_mode	('offset', 'rate', 'both')	What may jump at a break.	"offset"
break_var	float	Variance added at a break, in data units squared (and per year squared for the rate).	100.0
lag	int	Fixed-lag smoothing: epoch j uses data through j + lag (about 3 tau converges). None is the full RTS smoother.	None

Returns:

Type	Description
(ndarray, ndarray)	Smoothed displacement and velocity, (time, npix).

lowess

lowess(t_years: ndarray, y: ndarray, t_out: ndarray, *, bandwidth_days: float = 60.0, iterations: int = 2, min_points: int = 3) -> np.ndarray

Robust local-linear regression of (time, npix) series at t_out.

The window at each output time is bandwidth_days wide on each side, widened to the min_points nearest acquisitions where sampling is sparse; pixels with fewer than min_points valid values in it are NaN. Each of iterations passes reweights points by the bisquare of their residual over six times the pixel's median absolute residual (Cleveland 1979, https://doi.org/10.1080/01621459.1979.10481038).

q_from_timescale

q_from_timescale(tau_days: float, r: Union[float, ndarray]) -> Any

Return the process noise giving a smoothing timescale of tau_days.

For an integrated random walk the timescale where process noise balances measurement noise r is tau = (3 r / q)^(1/3) (in years). Setting q directly is how a filter silently returns a straight line (q too small for r).

Colored noise (one series)

Colored-noise trend estimation for one time series (maximum likelihood).

The deterministic model (a :class:~disp_xr.timeseries.trajectory.Trajectory or any design matrix) is fitted by generalized least squares inside the likelihood of a noise model, and the noise parameters are estimated by (restricted) maximum likelihood. The resulting sigmas account for temporal correlation, which ordinary least squares ignores: for GNSS and InSAR series OLS velocity sigmas are often 5 to 10 times too small.

Noise models:

"WN" White noise only (equivalent to OLS). "FNWN" Flicker (spectral index -1) plus white. "PLWN" Power law with a free spectral index kappa plus white. "RWFNWN" Random walk (-2), flicker and white.

Each likelihood evaluation factorizes an n x n covariance (O(n³)), so this is for one series or a group of pixels; full frames use the vectorized estimators of plan task T22. method="whittle" estimates the noise from the periodogram of the OLS residuals (O(n log n)) and makes one exact GLS solve.

The method follows Hector / HectorP (Bos et al.), written from the papers:

• Bos, M. S., Fernandes, R. M. S., Williams, S. D. P., & Bastos, L. (2013). Fast error analysis of continuous GNSS observations with missing data. Journal of Geodesy, 87(4), 351-360. https://doi.org/10.1007/s00190-012-0605-0
• Williams, S. D. P. (2003). The effect of coloured noise on the uncertainties of rates estimated from geodetic time series. Journal of Geodesy, 76(9-10), 483-494. https://doi.org/10.1007/s00190-002-0283-4
• Hosking, J. R. M. (1981). Fractional differencing. Biometrika, 68(1), 165-176. https://doi.org/10.1093/biomet/68.1.165

fit_noise

fit_noise(times: Sequence, y: Sequence[float], design: Optional[ndarray] = None, *, terms: Optional[Trajectory] = None, names: Optional[Sequence[str]] = None, model: NoiseModel = 'PLWN', rmle: bool = True, method: Literal['exact', 'whittle'] = 'exact', sampling_days: Optional[float] = None) -> NoiseResult

Fit a trajectory to one series under a colored-noise model.

Parameters:

Name	Type	Description	Default
times	sequence of datetime-like	Epochs; need not be regular. NaN values in y are dropped.	required
y	sequence of float	The series (any unit; outputs use the same unit).	required
design	ndarray	(time, coefficient) design matrix. By default built from terms.	None
terms	Trajectory	Model terms when design is not given; default offset, velocity and an annual cycle.	None
names	sequence of str	Coefficient names for a given design.	None
model	('PLWN', 'FNWN', 'RWFNWN', 'WN')	Noise model; see the module docstring.	"PLWN"
rmle	bool	Restricted maximum likelihood, which corrects the noise estimate for the degrees of freedom used by the deterministic model.	True
method	('exact', 'whittle')	Time-domain (R)MLE (O(n³) per evaluation), or the noise from the Whittle likelihood of the OLS residuals (O(n log n)) followed by one exact GLS solve. Whittle supports PLWN and WN.	"exact"
sampling_days	float	Nominal sampling interval of the regular grid the noise lives on. Default: the smallest spacing between epochs, in whole days.	None

Returns:

Type	Description
NoiseResult

Raises:

Type	Description
ValueError	If fewer than 10 epochs are valid, the model is rank deficient, or two epochs fall in one sampling interval.

Examples:

>>> series = ds.displacement[:, 10, 20]
>>> res = fit_noise(series.time.values, series.values)
>>> res.velocity, res.velocity_sigma, res.kappa

NoiseResult dataclass

A trajectory fitted under a colored-noise model.

Amplitudes are in the unit of the series; rates per year.

Examples:

>>> result = dxr.fit_noise(series, noise="PLWN")
>>> result.velocity, result.velocity_sigma, result.kappa

parameters property

parameters: dict[str, tuple[float, float]]

name -> (value, sigma) for every coefficient.

velocity property

velocity: float

Linear rate per year (NaN without a linear_trend term).

velocity_sigma property

velocity_sigma: float

Colored-noise sigma of the velocity.

Downsampling

downsample

Quality-weighted downsampling to a regular coarse grid (PRD S2, plan task T27).

Each output cell aggregates a factor x factor block of 30 m pixels:

• float layers: the weighted mean (stat="wmean") or weighted median (stat="wmedian"), with the weighted spread, the effective sample count n_eff = (Σw)² / Σw² and sigma = spread / sqrt(n_eff);
• integer mask and label layers: the most common value ("mode") or whether any pixel is set ("any"), never a mean.

With spread, each float layer also gets its spatial variation within the cell, all weighted with the same weights: <layer>_std (weighted standard deviation of the pixels), <layer>_rmse (weighted RMS of the pixels about the cell value; equal to the std for "wmean") and <layer>_mad (weighted median of the absolute deviations from the weighted median; times 1.4826 it estimates the std of Gaussian data).

Every cell also gets n_valid, valid_fraction (of the full block, so edge cells padded beyond the frame show how much they hold) and quality (mean static weight). Weighted-mean levels keep sum_w and sum_w2 so a coarser mean can be rebuilt from them; a weighted median cannot, so every median level is computed from 30 m. The sigma treats pixels as independent, which is far too small at 1 to 10 km; pass a variogram model (:mod:disp_xr.residuals) to also get <layer>_sigma_corr, which accounts for spatially correlated errors.

Odd factors (33, 67, 333 for 990 m, 2,010 m, 9,990 m) put each cell centre on a 30 m pixel centre, so the coarse transform is exact. Run after chaining and re-referencing: averaging series with different reference dates is wrong.

downsample

downsample(ds: Dataset, *, factor: Optional[int] = None, spacing: Optional[float] = None, stat: Stat = 'wmean', weights: Union[str, DataArray, None] = 'quality', layers: Optional[list[str]] = None, categorical: Optional[Mapping[str, Categorical]] = None, nan_threshold: float = 1.0, variogram: Optional[Mapping[str, VariogramModel]] = None, spread: Union[Spread, Sequence[Spread], None] = None) -> xr.Dataset

Aggregate a DISP dataset to coarse cells with quality weights.

Parameters:

Name	Type	Description	Default
ds	Dataset	Chained, re-referenced data with y and x (and optionally time) dimensions.	required
factor	int	Pixels per cell side. Odd factors keep cell centres on pixel centres (33, 67, 333 for about 1, 2 and 10 km at 30 m).	None
spacing	float	Cell size in the CRS unit instead of factor; rounded to the nearest odd factor.	None
stat	('wmean', 'wmedian')	Weighted mean or weighted median of float layers.	"wmean"
weights	('quality', DataArray or None)	"quality" (default): :func:disp_xr.quality_weights from the quality layers present, per epoch when the data have time. A DataArray (time, y, x) or (y, x) of weights, or None for equal weights.	'quality'
layers	list of str	Float layers to aggregate; default every float layer with y and x except quality layers used for the weights.	None
categorical	mapping	Mask and label layers and their reduction, "mode" or "any"; defaults in CATEGORICAL_DEFAULTS (mode for masks and labels, any for persistent scatterers and layover/shadow). Other integer layers (counts) are averaged like floats, ignoring their fill value.	None
nan_threshold	float	Largest share of missing pixels (of the full block) a cell may have; cells above it become NaN. 1.0 keeps any cell with data. Follows xarray-regrid's nan_threshold.	1.0
variogram	mapping of str to VariogramModel	A fitted variogram (:func:disp_xr.residuals.fit_variogram) per layer, estimated on that layer's residual field. Adds <layer>_sigma_corr = sqrt(sill / n_eff_corr) and n_eff_corr, the effective number of independent pixels of the cell mean given the spatial correlation (Hugonnet et al. 2022). The model is the same for every epoch.	None
spread	('std', 'rmse', 'mad')	Spatial variation of each float layer within the cell, weighted with the same weights: <layer>_std (weighted standard deviation), <layer>_rmse (weighted RMS about the cell value; equals the std for "wmean", and adds the mean-median offset for "wmedian") and <layer>_mad (weighted median absolute deviation from the weighted median; robust to outlier pixels). Per epoch when the data have time. Unlike <layer>_sigma, these describe the pixels, not the uncertainty of the cell value.	"std"

Returns:

Type	Description
Dataset	One variable per aggregated layer on the coarse grid, with <layer>_sigma and, shared by all, n_eff, n_valid, valid_fraction, quality, centroid_x, centroid_y, and for "wmean" sum_w and sum_w2; with variogram, <layer>_sigma_corr and n_eff_corr; with spread, <layer>_std, <layer>_rmse or <layer>_mad. Integer layers are reduced with mode or any. Lazy if the input is.

Examples:

>>> coarse = dxr.downsample(ds, factor=33)
>>> points = dxr.to_points(coarse)
>>> robust = dxr.downsample(
...     ds, factor=67, stat="wmedian", spread="mad"
... )
>>> robust.displacement_mad

to_points

to_points(coarse: Dataset, min_valid_fraction: float = 0.0) -> xr.Dataset

Return the occupied cells of a downsampled dataset as points.

The result has dimension point (and time if present), ordered by cell row then column, with the weighted centroid (centroid_x, centroid_y), the cell centre (x, y), the cell id (row * n_cols + col) and every cell variable. .to_dataframe() gives the table.

Examples:

>>> coarse = dxr.downsample(ds, factor=33)
>>> points = dxr.to_points(coarse, min_valid_fraction=0.5)
>>> points.to_dataframe()

to_points

to_points(coarse: Dataset, min_valid_fraction: float = 0.0) -> xr.Dataset

Return the occupied cells of a downsampled dataset as points.

The result has dimension point (and time if present), ordered by cell row then column, with the weighted centroid (centroid_x, centroid_y), the cell centre (x, y), the cell id (row * n_cols + col) and every cell variable. .to_dataframe() gives the table.

Examples:

>>> coarse = dxr.downsample(ds, factor=33)
>>> points = dxr.to_points(coarse, min_valid_fraction=0.5)
>>> points.to_dataframe()

downsample_hex

downsample_hex(ds: Dataset, *, spacing: Optional[float] = None, resolution: Optional[int] = None, grid: Literal['planar', 'h3'] = 'planar', stat: Stat = 'wmean', weights: Union[str, DataArray, None] = 'quality', layers: Optional[list[str]] = None, categorical: Optional[Mapping[str, Categorical]] = None, nan_threshold: float = 1.0, variogram: Optional[Mapping[str, VariogramModel]] = None, spread: Union[Spread, Sequence[Spread], None] = None) -> xr.Dataset

Aggregate a DISP dataset to hexagonal cells with quality weights.

Parameters:

Name	Type	Description	Default
ds	Dataset	Chained, re-referenced data with y and x (and optionally time) dimensions, in a projected CRS.	required
spacing	float	grid="planar": distance between neighboring cell centers, in the CRS unit (a hexagon covers 0.866 * spacing²).	None
resolution	int	grid="h3": H3 resolution, for example 7 (about 5.2 km²), 8 (0.74 km²) or 9 (0.11 km²).	None
grid	('planar', 'h3')	Planar hexagons in the frame's CRS anchored at its origin, or H3 global cells (the hex extra), which line up across frames and UTM zones. H3 cells are hexagons on the sphere, so their pixel counts vary slightly; valid_fraction and quality use each cell's own area.	"planar"
stat	Stat	As in :func:disp_xr.downsample. valid_fraction is relative to the expected pixel count of a full cell (its area over the pixel area), so it can slightly exceed 1.	'wmean'
weights	Stat	As in :func:disp_xr.downsample. valid_fraction is relative to the expected pixel count of a full cell (its area over the pixel area), so it can slightly exceed 1.	'wmean'
layers	Stat	As in :func:disp_xr.downsample. valid_fraction is relative to the expected pixel count of a full cell (its area over the pixel area), so it can slightly exceed 1.	'wmean'
categorical	Stat	As in :func:disp_xr.downsample. valid_fraction is relative to the expected pixel count of a full cell (its area over the pixel area), so it can slightly exceed 1.	'wmean'
nan_threshold	Stat	As in :func:disp_xr.downsample. valid_fraction is relative to the expected pixel count of a full cell (its area over the pixel area), so it can slightly exceed 1.	'wmean'
variogram	Stat	As in :func:disp_xr.downsample. valid_fraction is relative to the expected pixel count of a full cell (its area over the pixel area), so it can slightly exceed 1.	'wmean'
spread	Stat	As in :func:disp_xr.downsample. valid_fraction is relative to the expected pixel count of a full cell (its area over the pixel area), so it can slightly exceed 1.	'wmean'

Returns:

Type	Description
Dataset	Dimension point (and time), one point per cell holding at least one pixel of the grid; coordinates cell, the center x, y, and q, r (planar) or h3_cell, lon, lat (H3); the variables of :func:disp_xr.downsample and pixels. Lazy if the input is.

Examples:

>>> hexes = dxr.downsample_hex(ds, spacing=1075)
>>> hexes = dxr.downsample_hex(ds, grid="h3", resolution=8)
>>> dxr.plot.plot_hex(hexes)

quadtree

Quality-weighted quadtree downsampling (PRD S2 quadtree mode, plan task T28).

Cells of level k cover 2**k x 2**k pixels of the 30 m grid, anchored at its first pixel, from min_level (finest) to max_level (coarsest). A cell splits into its four children while its split score exceeds a threshold and every child keeps enough reliable pixels; the cells that do not split are the leaves. Smooth areas end up in large cells and areas with structure (deformation, unwrapping errors) in small ones.

The per-cell sums of :func:disp_xr.downsample (Σw, Σw², Σwv, Σwv², valid count) compose exactly: a level is the 2 x 2 block sum of the level below. So the data are read twice: once for the sums of the finest level (and from them every coarser level and every score), once for the values of the leaves. Leaf selection works on the small score maps, level by level, as boolean masks, with no recursion per cell.

Split metrics (metric):

• "displacement": RMS over epochs of the weighted spread of the pixels in the cell, in the unit of the data, so a threshold reads as "split where pixels in a cell differ by more than 5 mm";
• "velocity": weighted spread of a per-pixel velocity (pass velocity, for example from :func:disp_xr.fit_velocity), in its unit;
• "structure": the spread of the four children's means divided by their noise (the mean of their sigmas), both RMS over epochs: signal structure at the cell scale against noise, after the structure function of GPS Imaging (Hammond et al. 2016, https://doi.org/10.1002/2016JB013458). Unitless; split above about 2;
• or a callable metric(level, children) -> (y, x) DataArray that gets the level's sums as a Dataset and those of its children (None at min_level).

Plan task T29 compares the three and sets defaults.

Quadtree

Bases: NamedTuple

Leaves of a quadtree.

Examples:

>>> tree = dxr.quadtree(ds, n_points=5000)
>>> tree.points.level.to_series().value_counts()
>>> tree.to_raster("displacement", time=-1).plot.imshow()

points instance-attribute

points: Dataset

One point per leaf (dimension point, and time).

leaf_id instance-attribute

leaf_id: DataArray

(y, x) index into points of the leaf holding each pixel, -1 for pixels in no leaf (no reliable pixels in the cell).

to_raster

to_raster(layer: str = 'displacement', time=-1) -> xr.DataArray

Return a leaf variable painted back on the (y, x) grid.

time is an epoch index (default the last) or a date (nearest); pixels in no leaf are NaN.

quadtree

quadtree(ds: Dataset, *, threshold: Optional[float] = None, n_points: Optional[int] = None, metric: Metric = 'structure', min_level: int = 3, max_level: int = 8, min_reliable: float = 0.5, layer: str = 'displacement', weights: Union[str, DataArray, None] = 'quality', velocity: Optional[DataArray] = None) -> Quadtree

Downsample to quadtree leaves: small cells where the data vary.

Parameters:

Name	Type	Description	Default
ds	Dataset	Chained, re-referenced data with y, x (and time).	required
threshold	float	Split score above which a cell splits (see the module docstring for the unit of each metric).	None
n_points	int	Number of leaves wanted instead of a threshold: the threshold whose leaf count is closest is chosen (recorded as quadtree_threshold). Counts outside what min_level and max_level allow give the nearest extreme. Pass threshold or n_points, not both.	None
metric	('displacement', 'velocity', 'structure')	Split metric. Default "structure", the best of the three in the T29 trade study (docs/studies/split_metric.md); at the same number of points the regular grid of :func:downsample still reconstructs smooth signals better.	"displacement"
min_level	int	Finest and coarsest levels: cells of 2**min_level to 2**max_level pixels on a side (3 and 8 are 240 m and 7.7 km at 30 m).	3
max_level	int	Finest and coarsest levels: cells of 2**min_level to 2**max_level pixels on a side (3 and 8 are 240 m and 7.7 km at 30 m).	3
min_reliable	float	Share of a child's pixels inside the frame that must be reliable (static weight above 0 and data at some epoch) in each of the four children for a cell to split. Children outside the frame are ignored, so cells on the frame edge split like the others.	0.5
layer	str	Layer to aggregate and score.	'displacement'
weights	('quality', DataArray or None)	As in :func:disp_xr.downsample.	'quality'
velocity	DataArray	(y, x) velocity for metric="velocity".	None

Returns:

Type	Description
Quadtree	points: one point per leaf with level, size (meters), the cell center x, y, row, col (at its level), the weighted centroid, score, the weighted mean <layer> and its <layer>_sigma, n_eff, n_valid, valid_fraction and quality; lazy in the data. leaf_id: the (y, x) leaf of each pixel.

Examples:

>>> tree = dxr.quadtree(ds, threshold=0.005)
>>> tree.points.level.to_series().value_counts()
>>> tree = dxr.quadtree(ds, n_points=5000)

Quadtree

Bases: NamedTuple

Leaves of a quadtree.

Examples:

>>> tree = dxr.quadtree(ds, n_points=5000)
>>> tree.points.level.to_series().value_counts()
>>> tree.to_raster("displacement", time=-1).plot.imshow()

points instance-attribute

points: Dataset

One point per leaf (dimension point, and time).

leaf_id instance-attribute

leaf_id: DataArray

(y, x) index into points of the leaf holding each pixel, -1 for pixels in no leaf (no reliable pixels in the cell).

to_raster

to_raster(layer: str = 'displacement', time=-1) -> xr.DataArray

Return a leaf variable painted back on the (y, x) grid.

time is an epoch index (default the last) or a date (nearest); pixels in no leaf are NaN.

Quality-weighted downsampling to hexagonal cells (PRD S2 hexagon mode, T27.8).

Pixels are grouped into pointy-top hexagons of a planar lattice in the frame's CRS, anchored at the CRS origin, so overlapping frames of one UTM zone share cell ids. spacing is the distance between neighboring cell centers (the flat-to-flat width); a hexagon covers sqrt(3) / 2 * spacing², so spacing = 1.0746 * side matches the area of a square of that side.

Each pixel goes to the hexagon whose center is nearest, found with axial coordinates and cube rounding (Patel, "Hexagonal Grids", Red Blob Games). The per-cell statistics are those of :func:disp_xr.downsample (weighted mean or median, sigma, n_eff, valid fraction, quality, weighted centroid, mode or "any" for masks and labels), computed with grouped sums (one np.bincount per dask block) instead of block sums. Hexagons are not a raster, so the result is a (time, point) table in the :func:disp_xr.to_points layout.

Why hexagons: all six neighbors are equidistant and a hexagon is more compact than a square of equal area, so the sampling is more isotropic (Birch et al. 2007, https://doi.org/10.1016/j.ecolmodel.2007.03.041). The lattice is the one matplotlib.hexbin and PostGIS ST_HexagonGrid build (Carr et al. 1987, https://doi.org/10.1080/01621459.1987.10478445); here the cells carry quality weights and sigmas, lazily.

hex_axial

hex_axial(x: ndarray, y: ndarray, spacing: float) -> tuple[np.ndarray, np.ndarray]

Return the axial coordinates (q, r) of the hexagon holding each point.

The lattice is pointy-top, anchored at (0, 0), with centers at x = spacing * (q + r / 2), y = spacing * sqrt(3) / 2 * r.

hex_centers

hex_centers(q: ndarray, r: ndarray, spacing: float) -> tuple[np.ndarray, np.ndarray]

Return the center (x, y) of hexagons (q, r).

hexagon_vertices

hexagon_vertices(x: ndarray, y: ndarray, spacing: float) -> np.ndarray

Return the six vertices of pointy-top hexagons centered at (x, y).

Shape (n, 6, 2), for matplotlib.collections.PolyCollection or shapely polygons.

cell_polygons

cell_polygons(hexes: Dataset) -> list[np.ndarray]

Return the outline of every cell of downsample_hex in the data CRS.

One (k, 2) array of (x, y) vertices per point.

to_geodataframe

to_geodataframe(hexes: Dataset, time: Optional[Any] = -1) -> Any

Return the cells of downsample_hex as a GeoDataFrame of polygons.

Parameters:

Name	Type	Description	Default
hexes	Dataset	Output of :func:downsample_hex.	required
time	(int, date - like or None)	Epoch to take for time-dependent variables: an index (default the last), a date (nearest), or None for one row per cell and epoch.	-1

Returns:

Type	Description
GeoDataFrame	One row per cell (and epoch), the cell variables, the coordinates and the hexagon polygon, in the data CRS (needs the hex extra).

Examples:

>>> hexes = dxr.downsample_hex(ds, spacing=1000.0)
>>> gdf = dxr.to_geodataframe(hexes)
>>> gdf.to_file("hexes.gpkg")

Correlated-error sigma

empirical_variogram

empirical_variogram(field: DataArray, *, n_pairs: int = 200000, n_bins: int = 30, max_distance: Optional[float] = None, seed: int = 0) -> xr.Dataset

Estimate the empirical semivariogram of a 2D field from random pairs.

Each pair starts at a random valid pixel; its partner lies at a distance drawn log-uniformly between one pixel and max_distance, in a random direction. Every logarithmic distance bin then gets about the same number of pairs, unlike all pairs of a random sample, which leave short lags (where the cell sigma is decided) nearly empty.

Parameters:

Name	Type	Description	Default
field	DataArray	A (y, x) residual field, ideally zero-mean with no deformation signal (stable ground, or data minus a fitted model). NaN pixels are skipped. It is loaded into memory.	required
n_pairs	int	Number of pairs drawn (pairs whose partner is NaN or off the grid are dropped).	200000
n_bins	int	Number of logarithmic distance bins.	30
max_distance	float	Largest distance, in the CRS unit; default half the shorter side of the field.	None
seed	int	Random seed.	0

Returns:

Type	Description
Dataset	Dimension lag, with distance (mean pair distance in the bin, a coordinate), gamma (semivariance, half the mean squared difference) and count (pairs per bin). Empty bins are dropped.

Examples:

>>> residual = ds.displacement.isel(time=-1) - trend
>>> empirical = dxr.empirical_variogram(
...     residual, max_distance=20e3
... )

fit_variogram

fit_variogram(empirical: Dataset, models: Sequence[ModelKind] = ('spherical', 'spherical'), *, nugget: bool = True) -> VariogramModel

Fit a nugget and a sum of models to an empirical variogram.

Weighted least squares with the approximate standard error of each bin, gamma * sqrt(2 / count), so every bin counts by its relative misfit and short lags are not swamped by the large semivariances of long ones. Two spherical models (a short and a long range) are xDEM's default for elevation residuals and fit InSAR residuals (turbulent troposphere plus long-wavelength errors) well.

Parameters:

Name	Type	Description	Default
empirical	Dataset	Output of :func:empirical_variogram.	required
models	sequence of {"spherical", "exponential", "gaussian"}	One entry per component, from short to long range.	('spherical', 'spherical')
nugget	bool	Fit a nugget (uncorrelated pixel noise); otherwise it is 0.	True

Returns:

Type	Description
VariogramModel

Examples:

>>> model = dxr.fit_variogram(empirical, ("exponential",))
>>> coarse = dxr.downsample(
...     ds, factor=33, variogram={"displacement": model}
... )

VariogramModel dataclass

A nugget plus a sum of bounded variogram models.

gamma(h) = nugget + sum(sill_k * f_k(h / range_k)) for h > 0 and 0 at h = 0. Ranges are practical ranges (the semivariance reaches 95% of the partial sill for exponential and Gaussian models, all of it for spherical ones), in the CRS unit of the field.

Examples:

>>> model = dxr.fit_variogram(empirical)
>>> model.sill, model.ranges

components instance-attribute

components: tuple[tuple[str, float, float], ...]

(kind, partial sill, range) per model.

sill property

sill: float

Total variance of a pixel: nugget plus partial sills.

gamma

gamma(h: Union[float, ndarray]) -> np.ndarray

Return the semivariance at distance h.

correlation

correlation(h: Union[float, ndarray]) -> np.ndarray

Return the correlation of two pixels at distance h.

to_json

to_json() -> str

Return the model as a JSON string (for attributes).

Residuals: spatial variograms (T27.3b) and DISP against GNSS (T32).

:func:compare compares DISP with GNSS at the stations, with geepers computing the metrics (see its docstring). The rest of this docstring covers the variograms.

The sigma of a cell mean computed as spread / sqrt(n) assumes independent pixels. InSAR errors (troposphere, unwrapping, orbits) are correlated over kilometres, so at 1 to 10 km cells that sigma is far too small. Following xDEM (Hugonnet et al. 2022), this module

1. estimates an empirical variogram of a residual field (a field that should be zero: stable ground, or data minus a model) from a random sample of pixel pairs, with distances drawn log-uniformly so short lags are as well sampled as long ones (:func:empirical_variogram);
2. fits a sum of spherical, exponential or Gaussian models and a nugget (:func:fit_variogram);
3. turns the model into the effective number of independent samples of a block mean (:func:mean_block_correlation, :func:n_eff_correlated), which :func:disp_xr.downsample uses for <layer>_sigma_corr.

The residual semivariogram of plan task T32.1 reuses this code.

References

Hugonnet, R., et al. (2022). Uncertainty analysis of digital elevation models by spatial inference from stable terrain. IEEE JSTARS 15, 6456-6472. https://doi.org/10.1109/JSTARS.2022.3188922

Rolstad, C., Haug, T. and Denby, B. (2009). Spatially integrated geodetic glacier mass balance and its uncertainty based on geostatistical analysis. Journal of Glaciology 55(192), 666-680. https://doi.org/10.3189/002214309789470950

mean_block_correlation

mean_block_correlation(model: VariogramModel, factor: int, spacing: Union[float, tuple[float, float]]) -> float

Return the mean correlation between two distinct pixels of a block.

Exact for a full factor x factor block: the double sum over pixel pairs is a sum over the (2 factor - 1)² pixel offsets, each counted (factor - |dy|) (factor - |dx|) times.

Parameters:

Name	Type	Description	Default
model	VariogramModel	The fitted variogram.	required
factor	int	Pixels per block side.	required
spacing	float or (float, float)	Pixel spacing, or (dy, dx), in the unit of the model ranges.	required

Examples:

>>> rho = mean_block_correlation(model, factor=33, spacing=30.0)

mean_mask_correlation

mean_mask_correlation(model: VariogramModel, mask: ndarray, spacing: Union[float, tuple[float, float]]) -> float

Return the mean correlation between two distinct pixels of a cell shape.

mask marks the pixels of one full cell (a hexagon, for example). The number of pixel pairs at each offset is the autocorrelation of the mask, computed with an FFT, so the double sum over pairs costs one transform. :func:mean_block_correlation is the square-block case.

Examples:

>>> rho = mean_mask_correlation(model, cell_mask, spacing=30.0)

n_eff_correlated

n_eff_correlated(sum_w, sum_w2, rho_bar: float)

Return the effective sample count of a weighted mean of correlated pixels.

The variance of a weighted mean is sill * Σᵢⱼ wᵢ wⱼ ρᵢⱼ / (Σw)². Replacing each off-diagonal ρᵢⱼ by the block mean rho_bar gives n_eff = (Σw)² / (Σw² + ((Σw)² - Σw²) rho_bar): (Σw)² / Σw² for independent pixels, 1 for fully correlated ones. Exact for a full block with equal weights; it assumes valid pixels and weights are spread evenly over the block, so it overestimates n_eff (underestimates the sigma) for blocks whose valid pixels cluster, such as edge cells.

Examples:

>>> n_eff_correlated(sum_w=100.0, sum_w2=100.0, rho_bar=0.0)
100.0

Units, geometry and checks

to_units

to_units(obj: XarrayObj, unit: str) -> XarrayObj

Convert every length variable of obj to unit (m, cm or mm).

Variables whose units attribute is not a length (unitless, radians) are returned unchanged.

Raises:

Type	Description
ValueError	If unit is not a supported length unit.

Examples:

>>> import numpy as np, xarray as xr
>>> da = xr.DataArray(np.array([0.01]), attrs={"units": "meters"})
>>> float(to_units(da, "mm")[0])
10.0

enu_to_los

enu_to_los(east: ArrayT, north: ArrayT, up: ArrayT, los: Dataset | tuple[ArrayT, ArrayT, ArrayT]) -> ArrayT

Project east, north and up motion onto the line of sight.

Parameters:

Name	Type	Description	Default
east	(float, array or DataArray)	Motion components, in any unit; the result has the same unit.	required
north	(float, array or DataArray)	Motion components, in any unit; the result has the same unit.	required
up	(float, array or DataArray)	Motion components, in any unit; the result has the same unit.	required
los	Dataset or tuple	A dataset with los_east, los_north and los_up (as returned by :func:disp_xr.open_static) or the three components.	required

Returns:

Type	Description
Same type as the inputs	LOS motion, positive toward the satellite.

Examples:

>>> float(enu_to_los(0.0, 0.0, 1.0, (0.5, 0.0, 0.75**0.5)))  # pure uplift
0.8660254037844386

validate_disp

validate_disp(obj: Union[Dataset, DataArray], *, variables: Sequence[str] = (), require_reference_time: bool = False, per_pixel: bool = False) -> None

Raise if obj is not a usable DISP stack; list every problem found.

Checks dims (time, y, x), monotonic x and y coordinates, a CRS, the requested variables and, optionally, the reference_time coordinate. With per_pixel=True, warns when time is split across dask chunks, since per-pixel operations then rechunk (slow on big frames).

Raises:

Type	Description
DispValidationError	If any check fails.

Examples:

>>> dxr.validate_disp(ds, variables=["displacement"])