GNSS¶

GNSS stations next to DISP, with geepers as the engine (install the gnss extra). geepers loads, caches and reads the station series; disp_xr.gnss returns them as xarray and aligns them with DISP.

gps = dxr.gnss.load(ds, start="2017-01-01")            # (station, time) E/N/U
gps = dxr.gnss.to_los(gps, dxr.open_static(ds))        # + los, sigma_los
insar = dxr.gnss.sample(ds, gps, buffer_m=100)         # DISP at the stations
gps_on_disp = dxr.gnss.align(gps["los"], insar.time)   # DISP epochs and reference
residual = insar - gps_on_disp

GNSS stations next to DISP, with geepers as the engine (PRD F10, task T25).

geepers (https://github.com/opera-adt/geepers, the gnss extra) loads the GNSS series: its sources (UNR, UNR grid, sideshow) download, cache and read station files, and its sampler reads DISP around each station. This module adapts them to xarray and to DISP:

:func:get_source picks a geepers source; :func:stations lists the stations over a frame;
:func:load returns a (station, time) Dataset of east, north and up motion with sigmas and correlations, and station coordinates in lon/lat and in the frame's CRS;
:func:to_los projects them onto the line of sight with the static layers (sigma from geepers get_sigma_los);
:func:sample reads DISP at the stations (median in a window, geepers sample_insar);
:func:align puts GNSS on the DISP epochs and on the DISP reference date.
:func:centred stamps daily solutions at the centre of their UTC day.
:func:drop_relocations keeps each station's longest record between monument moves (applied by :func:load).

Typical use::

gps = dxr.gnss.load(ds, start="2017-01-01")
gps_los = dxr.gnss.to_los(gps, dxr.open_static(ds))
insar = dxr.gnss.sample(ds, gps, buffer_m=100)
gps_aligned = dxr.gnss.align(gps_los["los"], insar.time)

load ¶

load(target: Union[Dataset, DataArray, Bbox, list[str]], *, source: Union[str, Any] = 'unr', start: Optional[str] = None, end: Optional[str] = None, cache_dir: Optional[str] = None, max_workers: int = 8, max_step_m: Optional[float] = 10.0) -> xr.Dataset

Load GNSS series as a (station, time) Dataset.

Parameters:

Name	Type	Description	Default
`target`	`xr.Dataset, (west, south, east, north) or list of str`	A DISP dataset (its stations, with `x`, `y` in its CRS), a box in degrees, or station ids.	required
`source`	`str or geepers source`	See :func:`get_source`.	`'unr'`
`start`	`str`	Dates (ISO) to keep.	`None`
`end`	`str`	Dates (ISO) to keep.	`None`
`cache_dir`	`str`	geepers cache directory.	`None`
`max_workers`	`int`	Stations downloaded in parallel by geepers.	`8`
`max_step_m`	`float or None`	Largest 3-D position change between consecutive solutions that a station can make and keep its record. A larger one is not ground motion (near-field coseismic steps reach about 5 m) but a monument moved under the same id: the record is split there and only the segment with the most solutions is kept (:func:`drop_relocations`). None keeps every record.	`10.0`

Returns:

Type	Description
`Dataset`	`east`, `north`, `up` (meters, as geepers returns them) and `sigma_`, `corr_`; coordinates `station`, `time`, `lon`, `lat`, `alt` (and `x`, `y` for a dataset target), and `relocated` (True where a segment was dropped).

Examples:

>>> gps = dxr.gnss.load(ds, start="2017-01-01")
>>> gps = dxr.gnss.load(["P539", "P565"])

stations ¶

stations(target: Union[Dataset, DataArray, Bbox], *, source: Union[str, Any] = 'unr', cache_dir: Optional[str] = None)

Return the GNSS stations over a frame or a lon/lat box.

Parameters:

Name	Type	Description	Default
`target`	`(Dataset, DataArray or (west, south, east, north))`	A DISP dataset (stations inside its grid, with `x`, `y` in its CRS) or a box in degrees.	required
`source`	`str or geepers source`	See :func:`get_source`.	`'unr'`
`cache_dir`	`str`	geepers cache directory.	`None`

Returns:

Type	Description
`GeoDataFrame`	geepers station table (`id`, `lon`, `lat`, `alt`, geometry).

Examples:

>>> table = dxr.gnss.stations(ds)
>>> table = dxr.gnss.stations((-121.5, 35.5, -120.0, 37.0))

get_source ¶

get_source(name: Union[str, Any] = 'unr', cache_dir: Optional[str] = None) -> Any

Return a geepers GNSS source.

Parameters:

Name	Type	Description	Default
`name`	`('unr', 'unr_grid', 'sideshow')`	A source name, or an instance of a geepers `BaseGpsSource` subclass (returned as is).	`"unr"`
`cache_dir`	`str`	Where geepers caches downloads (its default otherwise).	`None`

Examples:

>>> source = get_source("unr", cache_dir="gnss-cache")

to_los ¶

to_los(gnss: Dataset, static: Dataset) -> xr.Dataset

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

The LOS unit vector is read at each station's pixel (nearest) in the static layers; the LOS sigma propagates the E/N/U sigmas and correlations (geepers get_sigma_los).

Parameters:

Name	Type	Description	Default
`gnss`	`Dataset`	Output of :func:`load` for a dataset target (with `x`, `y`).	required
`static`	`Dataset`	`los_east`, `los_north`, `los_up` on the frame grid (:func:`disp_xr.open_static`).	required

Returns:

Type	Description
`Dataset`	`gnss` with `los` and `sigma_los` (positive toward the satellite) and the LOS vector components per station.

Examples:

>>> los = dxr.gnss.to_los(gps, dxr.open_static(ds))
>>> los["los"].sel(station="P539").plot()

sample ¶

sample(ds: Dataset, gnss: Union[Dataset, Any], *, buffer_m: float = 0.0, layer: str = 'displacement') -> xr.DataArray

Read DISP at the stations: the median in a window around each one.

Uses geepers' XarrayReader.read_window (as its sample_insar does), with the station's pixel found by flooring (geepers rounds, which reads the neighboring pixel for half the stations). With buffer_m, geepers releases that have buffer_meters take a circular window of that radius; older ones (0.3) a square window of the same half-width in pixels.

Parameters:

Name	Type	Description	Default
`ds`	`Dataset`	DISP data with a CRS.	required
`gnss`	`Dataset or DataFrame`	Output of :func:`load` (`lon`, `lat` station coordinates), or a table with `lon`, `lat` (and `id`) columns.	required
`buffer_m`	`float`	Window radius in meters (0 reads the station's pixel).	`0.0`
`layer`	`str`	Layer to read.	`'displacement'`

Returns:

Type	Description
`DataArray`	`(station, time)` DISP at each station.

Examples:

>>> insar = dxr.gnss.sample(ds, gps, buffer_m=90)

align ¶

align(gnss: DataArray, times: Union[DataArray, DatetimeIndex, ndarray], *, reference: Optional[Any] = None, method: Literal['nearest', 'linear'] = 'nearest', tolerance: str = '1D') -> xr.DataArray

Put a GNSS series on the DISP epochs and the DISP reference date.

Parameters:

Name	Type	Description	Default
`gnss`	`DataArray`	`(station, time)` GNSS series, for example `to_los(...)["los"]`.	required
`times`	`array of datetimes`	DISP epochs (`ds.time`).	required
`reference`	`date - like`	DISP reference date; default the first epoch (a chained stack is relative to it).	`None`
`method`	`('nearest', 'linear')`	Nearest GNSS solution within `tolerance`, or linear interpolation. Daily solutions count as stamped at noon (:func:`centred`), so an epoch takes the solution of its own UTC day.	`"nearest"`
`tolerance`	`str`	Largest distance in time for `"nearest"` (pandas offset).	`'1D'`

Returns:

Type	Description
`DataArray`	GNSS at `times` minus its value at `reference`; NaN where the GNSS series does not reach an epoch. A station with no solution at `reference` (installed later, or a gap that day) is referenced to its first epoch instead, so its residual against DISP carries a constant offset, which rates and mean-removed statistics ignore. The `gnss_reference` coordinate gives each station's reference epoch.

Examples:

>>> on_epochs = dxr.gnss.align(los["los"], ds.time)

DISP against GNSS¶

dxr.residuals.compare reads DISP at the stations (or the downsampled cell holding each one), puts GNSS on the DISP epochs and reference, and lets geepers compute the metrics.

result = dxr.residuals.compare(ds, gps, static, buffer_m=100,
                               requirement=lambda d: (3 + 0.5 * d**0.5) / 1000)
result.stations          # velocities, normalized residuals, chi2, outliers
result.profile           # pairwise RMSE against distance, with the requirement
result.plot_velocities()

compare ¶

compare(data: Any, gnss: Dataset, static: Dataset, *, buffer_m: float = 0.0, layer: str = 'displacement', velocity_method: str = 'lsq', tolerance: str = '1D', min_common_dates: int = 3, min_stations: int = 3, min_epochs: int = 10, min_span_years: float = 2.0, n_bins: int = 10, requirement: Optional[Any] = DISP_S1_DISPLACEMENT, velocity_requirement: Optional[Requirement] = DISP_S1_VELOCITY, outlier_sigma: float = 3.0, epoch_flag: float = 3.0, disp_sigma: Optional[float] = None) -> Comparison

Compare DISP with GNSS at the stations (or the cells holding them).

geepers computes the metrics (velocities from one estimator applied to both series, pairwise differential RMSE against distance, its binned profile and semivariogram, epoch RMSE, agreement statistics); disp-xr prepares the series and adds normalized velocity residuals, reduced χ² per station and flags.

Parameters:

Name	Type	Description	Default
`data`	`xr.Dataset, point table or Quadtree`	DISP at 30 m (read at each station with :func:`disp_xr.gnss.sample`), or downsampled cells: a point table from :func:`disp_xr.to_points` or :func:`disp_xr.downsample_hex` (the nearest cell center) or a :class:`disp_xr.Quadtree` (the leaf holding the station).	required
`gnss`	`Dataset`	:func:`disp_xr.gnss.load` output for the same frame.	required
`static`	`Dataset`	LOS layers for :func:`disp_xr.gnss.to_los`.	required
`buffer_m`	`float`	Window radius for sampling 30 m data.	`0.0`
`layer`	`str`	DISP layer.	`'displacement'`
`velocity_method`	`('lsq', 'midas', 'trend')`	Estimator for both series (geepers `fit_velocities`); MIDAS needs about a year of data.	`"lsq"`
`tolerance`	`str`	Largest gap between a DISP epoch and a GNSS day.	`'1D'`
`min_common_dates`	`int`	Passed to geepers (pairs, epochs, profile).	`3`
`min_stations`	`int`	Passed to geepers (pairs, epochs, profile).	`3`
`n_bins`	`int`	Passed to geepers (pairs, epochs, profile).	`3`
`min_epochs`	`int`	Least number of epochs with both series, and least time between the first and last of them, for a station's velocities to count. A rate from a few epochs says nothing about the product, yet it would enter the velocity residuals, the outlier flag and the GNSS tie. Below either, the velocities and their sigmas are NaN and `velocity_usable` is False; the station's displacement series still enters the pairs and the epoch RMSE.	`10`
`min_span_years`	`int`	Least number of epochs with both series, and least time between the first and last of them, for a station's velocities to count. A rate from a few epochs says nothing about the product, yet it would enter the velocity residuals, the outlier flag and the GNSS tie. Below either, the velocities and their sigmas are NaN and `velocity_usable` is False; the station's displacement series still enters the pairs and the epoch RMSE.	`10`
`requirement`	`Requirement or callable`	Threshold of the pairwise relative displacement RMSE against distance, drawn on the profile with geepers' `fraction_passing`. Default :data:`DISP_S1_DISPLACEMENT`. A :class:`Requirement` follows the units of `layer`; a plain `requirement(distance_km)` must return the units of the data. None for no requirement.	`DISP_S1_DISPLACEMENT`
`velocity_requirement`	`Requirement`	Threshold of the relative velocity: every station pair's difference of velocity residuals is scored against it (`velocity_pairs`, `summary["velocity_fraction_passing"]` over the pairs inside its distance range). Default :data:`DISP_S1_VELOCITY`; None to skip.	`DISP_S1_VELOCITY`
`outlier_sigma`	`float`	Stations with `\|normalized_residual\|` or `sqrt(chi2)` above this are flagged.	`3.0`
`disp_sigma`	`float`	Noise of one DISP epoch at a station, in the units of `layer`, added in quadrature to the GNSS sigma in `chi2`. Default: estimated from all stations' residuals (1.4826 MAD, less the median GNSS variance), so one bad station does not raise it.	`None`
`epoch_flag`	`float`	Epochs whose RMSE exceeds the median plus this many MADs are flagged.	`3.0`

Returns:

Type	Description
`Comparison`

Examples:

>>> gps = dxr.gnss.load(ds)
>>> result = compare(ds, gps, dxr.open_static(ds))
>>> result.stations[["velocity_residual", "outlier"]]
>>> cells = compare(dxr.quadtree(ds, n_points=5000), gps, static)

Comparison `dataclass` ¶

DISP against GNSS at the stations; the metrics come from geepers.

Attributes are filled by :func:compare.

Examples:

>>> result = compare(ds, gps, static)
>>> result.summary["velocity_fraction_passing"]
>>> result.plot_velocities()

series `instance-attribute` ¶

series: Dataset

(station, time): insar, gnss (on the DISP epochs and reference), residual (insar - gnss) and sigma_gnss.

stations `instance-attribute` ¶

stations: Any

DataFrame per station: velocities of both series from the same estimator (geepers fit_velocities) with sigmas, velocity_residual, relative_residual (less the datum offset), velocity_sigma (combined), normalized_residual (relative residual over its sigma), chi2 (reduced, of the displacement residual against the GNSS and DISP sigmas in quadrature), n_dates, outlier, velocity_usable (enough shared epochs and time span for the velocities, which are NaN otherwise) and reference_epoch (the first epoch with both series, where both are zero: the DISP reference date unless the station has no GNSS that day).

pairs `instance-attribute` ¶

pairs: Any

geepers pairwise_differential_rmse: RMSE of the relative misfit of every station pair against distance.

profile `instance-attribute` ¶

profile: Any

geepers binned_rmse_profile of pairs (with requirement and fraction_passing when a requirement is given).

epochs `instance-attribute` ¶

epochs: Any

geepers epoch_rmse per DISP epoch, with flag for epochs above the median plus epoch_flag x MAD of the RMSE.

semivariogram `instance-attribute` ¶

semivariogram: Any

geepers misfit_semivariogram of pairs.

summary `instance-attribute` ¶

summary: dict

geepers comparison_stats of the station velocities (bias, MAD, RMSE, r², n), datum_offset (median velocity residual: DISP is relative to its reference, GNSS absolute), disp_sigma (the DISP noise used in chi2), with velocity_fraction_passing and the requirement names when requirements are given.

velocity_pairs `class-attribute` `instance-attribute` ¶

velocity_pairs: Any = None

Per station pair: distance, difference of the velocity residuals, the velocity requirement there, applies (inside its distance range) and passing; None without a velocity requirement.

plot_velocities ¶

plot_velocities(**kwargs: Any) -> Any

Scatter GNSS against DISP velocities (geepers plot_velocity_scatter).

to_dataset ¶

to_dataset() -> xr.Dataset

Return the station table and the summary as a station Dataset.

Every column of stations becomes a variable, with the station coordinates (lon, lat, x, y) from series; the summary is stored as JSON in the summary attribute.

Examples:

>>> result.to_dataset().velocity_residual.plot()

to_zarr ¶

to_zarr(store: Any, group: str = 'gnss', **kwargs: Any) -> None

Write :meth:to_dataset as a group of store (next to a cube).

disp_xr serve reads this group for its gnss route.

Examples:

>>> result.to_zarr("frame.zarr")

plot_semivariogram ¶

plot_semivariogram(**kwargs: Any) -> Any

Plot the misfit semivariogram (geepers plot_semivariogram).

residual_spectra ¶

residual_spectra(source: Any, *, frequencies: Optional[ndarray] = None, detrend: bool = True) -> xr.DataArray

Return the Lomb-Scargle power spectrum of each station's residual series.

Lomb-Scargle handles the irregular DISP epochs (gaps, 12 and 24 day revisits) without resampling. Peaks at 1 and 2 cycles per year are seasonal signal one series has and the other lacks; a red spectrum (power rising to low frequencies) is correlated error.

Parameters:

Name	Type	Description	Default
`source`	`Comparison or DataArray`	:func:`compare` output, or a `(station, time)` residual.	required
`frequencies`	`array`	Frequencies in cycles per year; default from one cycle over the record to the Nyquist frequency of the median epoch spacing.	`None`
`detrend`	`bool`	Remove each series' linear trend first.	`True`

Returns:

Type	Description
`DataArray`	`(station, frequency)` power (normalized, 0 to 1).

Examples:

>>> power = residual_spectra(result)

common_mode ¶

common_mode(source: Any, n_components: int = 1, min_coverage: float = 0.5) -> xr.Dataset

Return the common-mode residual across stations (geepers estimate_cme).

A signal shared by all stations' residuals (a reference or datum shift, a frame-wide atmospheric or orbit error) is a product error that the station-by-station metrics miss.

Parameters:

Name	Type	Description	Default
`source`	`Comparison or DataArray`	:func:`compare` output, or a `(station, time)` residual.	required
`n_components`	`int`	Principal components to keep.	`1`
`min_coverage`	`float`	Stations observed in fewer than this share of epochs are left out.	`0.5`

Returns:

Type	Description
`Dataset`	`common_mode` `(time, component)`, `spatial_response` `(station, component)`, `cleaned` `(station, time)` (residual minus the common mode) and `explained_variance` `(component)`.

Examples:

>>> cme = common_mode(result, n_components=1)

An EOF analysis of trajectory residuals with xeofs (eof extra) is in examples/eof_residuals.py; it runs in its own CI job.

GNSS¶

load ¶

stations ¶

get_source ¶

to_los ¶

sample ¶

align ¶

DISP against GNSS¶

compare ¶

Comparison dataclass ¶

series instance-attribute ¶

stations instance-attribute ¶

pairs instance-attribute ¶

profile instance-attribute ¶

epochs instance-attribute ¶

semivariogram instance-attribute ¶

summary instance-attribute ¶

velocity_pairs class-attribute instance-attribute ¶

plot_velocities ¶

to_dataset ¶

to_zarr ¶

plot_semivariogram ¶

residual_spectra ¶

common_mode ¶

Comparison `dataclass` ¶

series `instance-attribute` ¶

stations `instance-attribute` ¶

pairs `instance-attribute` ¶

profile `instance-attribute` ¶

epochs `instance-attribute` ¶

semivariogram `instance-attribute` ¶

summary `instance-attribute` ¶

velocity_pairs `class-attribute` `instance-attribute` ¶