GeoMotion: Decomposition, Averaging, Denoising and Mapping of EGMS Ground-Motion Data

Keywords: InSAR; EGMS; Sentinel-1; persistent scatterers; line-of-sight decomposition; error propagation; kernel smoothing; subsidence; ground motion

1Introduction

Interferometric synthetic aperture radar (InSAR) measures the change in the radar path length between repeated satellite acquisitions and therefore one component of the three-dimensional ground displacement (Hanssen, 2001). Multi-temporal techniques restrict the analysis to scatterers whose phase remains stable over the whole observation period, which yields displacement time series with millimetre precision at the location of buildings, infrastructure and exposed rock (Ferretti et al., 2001). The European Ground Motion Service (EGMS) of the Copernicus Land Monitoring Service applies such processing to all Sentinel-1 acquisitions over Europe and publishes the results in harmonised products (Crosetto et al., 2020; Costantini et al., 2021).

GeoMotion is a web application that queries the public EGMS endpoints through a caching proxy and performs all analyses in the browser. It offers three tools: a single-point probe, which selects the nearest ascending and descending measurement points and decomposes their LOS series into east–west and vertical motion; a multi-point selection, which averages several points per orbit before the same decomposition; and a displacement map, which animates the vertical and east–west Level 3 components inside a user-drawn polygon. The purpose of this note is to state the mathematical model behind each of these computations in the form in which it is implemented, including constants, thresholds and fallbacks. Where the implementation simplifies, we say so, and we distinguish throughout between quantities that EGMS measures, quantities that GeoMotion derives from them by exact transformations, and quantities that rest on a model assumption.

The note is organised as follows. Section 2 describes the EGMS products and their sampling. Section 3 derives the LOS observation model and analyses the observability of the north component. Section 4 treats the decomposition at a single probe point, including time alignment, error propagation and the derived velocities. Section 5 covers multi-point averaging and its uncertainty, Section 6 the accuracy bands of the time-series chart and Section 7 temporal denoising. Section 8 describes the displacement map, including the gradient model for horizontal motion. Section 9 compares the decomposition with EGMS Level 3, Section 10 discusses limitations and Section 11 concludes. Appendix A lists all parameters and Appendix B maps the equations to the source modules.

1.1Conventions and notation

We express ground displacement in a local east–north–up frame as $\vD = (E, N, U)\T$, in millimetres, with $E$ positive eastward, $N$ positive northward and $U$ positive upward, so that subsidence corresponds to $U \lt 0$. Subscripts $\mathrm a$ and $\mathrm d$ denote the ascending and the descending geometry. Time is measured in days, and one year is defined as $365.25$ days wherever a rate is expressed per year. A displacement time series is a sequence of pairs $(t_j, x_j)$; missing values are excluded from every sum, mean and fit. The symbols used throughout are collected in Table 1.

2Data

2.1EGMS products

GeoMotion uses two of the EGMS product levels (Copernicus Land Monitoring Service, n.d.-b; Crosetto et al., 2020). The calibrated product (Level 2B) contains, for each Sentinel-1 track, the measurement points identified by persistent- and distributed-scatterer processing. Each point carries a LOS displacement time series in millimetres and metadata, of which the application uses the LOS unit vector (los_east, los_north, los_up), the mean LOS velocity (mean_velocity) and its standard deviation (mean_velocity_std), the root-mean-square deviation of the time series from the EGMS displacement model (rmse_ts) and the temporal coherence. Level 2B differs from the basic Level 2A product by a calibration against a GNSS-based velocity model, which ties the otherwise relative InSAR measurements to a common reference (Copernicus Land Monitoring Service, n.d.-b). A track is identified by its orbit direction and relative orbit number; for example, A21-088 is the ascending track with relative orbit 88.

The ortho product (Level 3) provides the vertical and the east–west displacement components on a regular 100 m grid in the ETRS89 Lambert azimuthal equal-area projection (EPSG:3035). EGMS derives it from the Level 2B data of all overlapping ascending and descending tracks; the north component is not resolved (Copernicus Land Monitoring Service, n.d.-a). Each grid cell carries a time series per component and the corresponding mean velocity.

2.2Release and temporal sampling

EGMS publishes a new release approximately once per year. The catalogue may therefore list several releases of the same track, and the application retains, for every track (identified by orbit and relative orbit number) and for every Level 3 component, the dataset whose period ends latest. At the time of writing, all 85 datasets belong to the 2020–2024 release, which covers January 2020 to December 2024, so every analysis in this note refers to that period.

The temporal sampling of Level 2B is irregular. Until December 2021 both Sentinel-1A and Sentinel-1B were operational and each track was imaged every 6 days; after the loss of Sentinel-1B the revisit interval became 12 days (Torres et al., 2012). A typical point of track A21-088 has 209 epochs, with intervals of 6 and 12 days and occasional 24-day gaps. Level 3, in contrast, is sampled on a regular 6-day lattice: in the current release every cell ends on 18 December 2024 and starts on 8, 14 or 20 January 2020, which gives 302, 301 or 300 epochs. This difference matters for the time alignment of Section 4.2 and for the time cube of Section 8.2.

2.3Tile decoding and coordinates

Point positions are obtained from the vector tiles that the official EGMS viewer uses. The tiles store Cartesian coordinates on a sphere of radius 128 in arbitrary units, as the sum $\vp = \mathbf c + \mathbf v$ of a tile centre $\mathbf c$ and a per-point offset $\mathbf v$. With $\vp = (x, y, z)\T$, where $y$ points to the north pole and $z$ to the intersection of the equator with the prime meridian, longitude $\lambda$ and latitude $\varphi$ follow from

We treat $\lambda$ and $\varphi$ as WGS 84 coordinates. For display, the application converts them to Web Mercator (EPSG:3857) with

and it measures distances between points with the haversine formula on a sphere of radius $R = 6\,371\,008.8$ m. Since the Mercator projection is conformal with scale factor $\sec\varphi$, map lengths at the latitude of Brussels ($\varphi \approx 51^\circ$) exceed ground lengths by a factor of about 1.59; Section 8 accounts for this where it matters.

3Line-of-sight geometry

3.1The observation equation

A radar interferometer is sensitive only to the change in distance between the satellite and the scatterer. For displacements that are small compared with the satellite range, this change equals the projection of the ground displacement onto the line of sight. Let $\los = (e, n, u)\T$ denote the unit vector pointing from the ground towards the satellite. The LOS displacement is then

so that $d \gt 0$ indicates motion towards the satellite, that is, a shortening of the range. EGMS publishes the three components of $\los$ per measurement point, rounded to three decimals; the application renormalises them as $\los \leftarrow \los/\lVert\los\rVert$ before use.

3.2The LOS vector from incidence angle and heading

Sentinel-1 observes to the right of its flight direction. If $\alpha$ denotes the heading, measured clockwise from north, and $\theta$ the incidence angle at the target, the look direction has azimuth $\alpha + 90^\circ$, and the horizontal projection of the ground-to-satellite vector points towards azimuth $\alpha - 90^\circ$. A horizontal unit vector at azimuth $\beta$ has east and north components $(\sin\beta, \cos\beta)$; substituting $\beta = \alpha - 90^\circ$ and scaling by $\sin\theta$ gives

GeoMotion does not evaluate (4), because the vectors are part of the metadata, but the relation reproduces the published vectors to within the rounding of 0.001 for all tracks we have examined. Table 2 lists the geometry of the four tracks that cover the centre of Brussels, and Figure 1 illustrates the configuration.

3.3Why the north–south component is not observed

The sensitivity of a LOS measurement to horizontal motion depends on the direction of that motion. For a horizontal unit displacement at azimuth $\beta$, $\hat{\vh}(\beta) = (\sin\beta, \cos\beta, 0)\T$, equations (3) and (4) give

The sensitivity vanishes for $\beta = \alpha$ and $\beta = \alpha + 180^\circ$: a radar acquisition in zero-Doppler geometry is blind to motion parallel to the flight direction (Hanssen, 2001). Sentinel-1 flies a near-polar, sun-synchronous orbit with an inclination of 98.18° (Torres et al., 2012), and at the latitude of Belgium its heading deviates from the meridian by only about 10°. For the four tracks of Table 2, the blind directions lie between 11.4° west and 10.7° east of north (Figure 2). Combining an ascending and a descending track does not remove the blindness, because their blind directions differ by only about 20°. The same fact explains the small north components of the LOS vectors: by (4), $|n| = \sin\theta\,|\sin\alpha| \le 0.11$ for all four tracks.

We therefore adopt the following assumption, which is also made in the EGMS Level 3 product (Copernicus Land Monitoring Service, n.d.-a).

If the assumption is violated, the estimates absorb the neglected term. Writing the decomposition of Section 4.3 as $\vA\vx = \vb$ with $\vx = (E, U)\T$, a true north displacement $N$ adds $N(n\sa, n\sd)\T$ to $\vb$, and the estimate becomes

For the pair A20-161/D18-110, $\vA^{-1}(n\sa, n\sd)\T = (0.009, -0.137)\T$: one millimetre of northward motion appears as 0.14 mm of apparent subsidence and 0.01 mm of eastward motion. Because $n\sa \approx n\sd$, the neglected term acts as a common-mode LOS bias and leaks almost entirely into the vertical component.

4Decomposition at a single probe point

The probe tool converts one ascending and one descending LOS time series into east–west and vertical time series and velocities. The procedure consists of point selection (Section 4.1), alignment of the two series in time (Section 4.2), the solution of a $2\times2$ system per epoch (Section 4.3), error propagation (Section 4.4), the estimation of velocities and totals (Section 4.5) and the presentation of the result (Section 4.6).

4.1Point selection

Given a probe location $(\lambda_0, \varphi_0)$, the application first identifies candidate tracks whose coverage spheres contain the location; this test is coarse and typically returns more than ten tracks per orbit, most of them without points at the location. For each orbit separately and for each search radius $r_{\max} \in \{100, 250, 500\}$ m in turn, it loads all points of all candidate tracks inside the bounding box of half-width $r_{\max}$ and ranks them with the local equirectangular metric

The best-ranked point of each track is accepted if its haversine distance does not exceed $r_{\max}$; the nearest accepted point over all tracks of the orbit is selected, and the search stops at the first radius that yields a point. The ascending and the descending point are thus selected independently. They are, in general, different scatterers, possibly on different structures, up to 500 m apart; the reported distances allow the user to judge the following assumption.

4.2Time alignment

The two series are not acquired on the same days. At Brussels, for instance, the descending track D18-110 is imaged three days before the ascending track A20-161 (Figure 3). Instead of pairing nearest dates, the application interpolates both series onto a common regular grid over their overlapping period,

where $t^{\mathrm a}_{1}$ and $t^{\mathrm a}_{\mathrm{last}}$ are the first and last epochs of the ascending series with a finite value, and likewise for the descending series. Each series is interpolated linearly between the bracketing acquisitions $t_j \le t_k \le t_{j+1}$,

without extrapolation. Each series is then referred to a robust estimate of its level at the first grid epoch rather than to the single, noisy acquisition at that epoch,

where $\hat x_{\mathrm{LL}}(t;\sigma)$ is the Gaussian-weighted local-linear estimate of (32) (Section 7.1), evaluated on the original acquisitions of the series with a kernel scale $\sigma_{\mathrm{ref}}$. It is the intercept at $t_0$ of a straight line fitted to the acquisitions within $3\sigma_{\mathrm{ref}} = 180$ days, each weighted by its Gaussian distance from $t_0$.

Because all acquisitions of a track fall on a 6-day lattice, the grid (8) coincides with the acquisitions of the later-starting series as long as the 6-day revisit lasts, whereas the other series is always interpolated (at mid-interval when the offset is three days). After December 2021 every second grid epoch is interpolated for both series. Interpolation at mid-interval averages two acquisitions, so for independent errors of variance $\sigma^2$ the interpolated value has variance $\sigma^2/2$, and consecutive grid values become correlated, which the uncertainty of the velocities in Section 4.4 takes into account.

4.3The two-by-two system and Cramer's rule

Under Assumption 1, equation (3) written for both geometries at grid epoch $t_k$ gives the linear system

with determinant

Cramer's rule yields the east–west and vertical displacement at every grid epoch at which both referenced series are defined,

The application refuses the decomposition if the overlapping period contains fewer than two grid epochs or if the geometry is too poorly conditioned, as quantified below in (16). Because both series are referred to robust levels rather than to their values at $t_0$, $E_0$ and $U_0$ are small but in general not exactly zero.

The determinant has a simple geometric meaning. Let $\vp\sa = (e\sa, u\sa)\T$ and $\vp\sd = (e\sd, u\sd)\T$ be the projections of the two LOS vectors onto the east–up plane, and let $\gamma$ be the angle between them. Then

the signed area of the parallelogram spanned by $\vp\sa$ and $\vp\sd$ (negative because $\vp\sd$ lies clockwise of $\vp\sa$). Each row of (11) confines $(E_k, U_k)$ to a line perpendicular to the corresponding $\vp$; the solution is the intersection of the two lines, and $|\det\vA|$ measures how transversally they intersect (Figure 4). Substituting (4) with $\cos\alpha\sa \gt 0$ and $\cos\alpha\sd \lt 0$ gives

where the approximation holds because $|\cos\alpha| \approx 0.98$ for both orbits. For A20-161/D18-110, $\theta\sa + \theta\sd = 75.5^\circ$ and $\det\vA = -0.953$. Over Belgium the incidence angles range from about 30° to 45°, so $|\det\vA|$ lies between about 0.85 and 1.0 for every opposite-orbit pair.

A threshold on $|\det\vA|$ alone does not separate good from poor geometries, because it also depends on the lengths of $\vp\sa$ and $\vp\sd$: two ascending tracks with different incidence angles (A21-088 and A20-161) give $|\det\vA| = 0.16$ although their projected vectors are only 9.1° apart. The application therefore tests the angle itself and requires

where $\kappa(\tilde\vA)$ is the condition number of $\vA$ with its rows normalised to unit length. Opposite-orbit pairs over Belgium have $\gamma$ between about 60° and 90° and pass; any pair of the same orbit has $\gamma \lt 20^\circ$ and is rejected.

4.4Error propagation

We propagate the published uncertainties to first order under the following assumption.

The solution $\vx_k = \vA^{-1}\vb_k$ is linear in the observations, so its covariance matrix is

Evaluating the diagonal of (17) gives the standard deviations used by the application,

and the off-diagonal element, which the application does not use, is

For A20-161/D18-110 with $\sigma\sa = \sigma\sd = \sigma$, equation (18) gives $\sigma_E = 1.17\,\sigma$ and $\sigma_U = 0.90\,\sigma$, and (19) gives a correlation coefficient of only 0.02, because $e\sa$ and $e\sd$ have opposite signs. The vertical component is better determined than the east–west component: both geometries are steep ($u \approx 0.76$–$0.82$), so the vertical signal enters both observations with the same sign, whereas the east–west signal enters with opposite signs and smaller weights.

The application applies (18) twice: with the per-epoch noise rmse_ts, which yields the bands $\pm\sigma_E$ and $\pm\sigma_U$ of the decomposed series (Section 6), and with the velocity standard deviations mean_velocity_std, which yields the propagated velocity uncertainties $\sigma^{\mathrm{prop}}_{v_E}$ and $\sigma^{\mathrm{prop}}_{v_U}$. The latter enter the velocity uncertainty of Section 4.5 together with an estimate from the decomposed series itself.

4.5Velocities and totals

Velocities are the ordinary least-squares slopes of the decomposed series against time in years. With $\tau_k = (t_k - t_0)/\yr$ and the means $\bar\tau$ and $\bar x$ taken over the epochs with a finite value,

applied to $x_k = E_k$ and $x_k = U_k$, which gives $v_E$ and $v_U$.

The propagated uncertainty $\sigma^{\mathrm{prop}}_v$ of Section 4.4 is optimistic in three respects: the published mean_velocity_std refers to the EGMS fit of the full LOS series rather than to the overlapping period, it is typically only 0.1–0.2 mm/yr, and it does not reflect the serial correlation that interpolation (Section 4.2) and residual seasonal motion introduce in the decomposed series. The application therefore also computes the standard error of (20) from the residuals $\varepsilon_k$ of the fit, corrected for lag-one autocorrelation following Santer et al. (2000),

where $n$ is the number of grid epochs with a finite value and the lag-one sum runs over consecutive finite epochs. The reported uncertainty is the larger of the two estimates,

The second term accounts for the actual scatter of the decomposed series, including seasonal motion, and for its serial correlation; the first protects against an optimistic fit of a short or smooth series. At the Brussels probe point of Figure 1 (A20-161 and D18-110), for example, the propagated values are 0.12 mm/yr east–west and 0.09 mm/yr vertically, whereas (21) gives 0.22 and 0.19 mm/yr, which are therefore reported. Errors that are common to a whole region, such as residual calibration or reference-frame errors, remain outside both estimates.

The totals shown in the interpretation are the displacements implied by the fitted trend over the length $T = (t_{K-1} - t_0)/\yr$ of the grid,

We prefer (23) to the difference between the last and the first decomposed epoch because that difference carries the full noise of two single epochs: by (18), with $\sigma = 3$ mm its standard deviation is about $\sqrt2 \times 1.17 \times 3 \approx 5$ mm in the east–west component, whereas the trend uses all epochs. The price is that non-linear motion, such as acceleration or a seasonal cycle, is summarised by its best linear approximation over the period. Where a model-free total is needed, namely in the fallback of the significance test below, the application uses the difference of robust levels, $\Delta^{\mathrm{L}} x = L_x(t_{K-1}) - L_x(t_0)$, with $L_x$ the local-linear level of the decomposed series as in (10).

4.6Significance and presentation

A decomposed component is flagged as not significant if

which corresponds to a two-sided test at approximately the 5 % level for Gaussian errors. If $\sigma_v$ is unavailable, the test is applied to the robust total $\Delta^{\mathrm L}x$ with the per-epoch $\sigma_E$ or $\sigma_U$. The same rule is used in the map callout, in the side panel and in the chart. The interpretation is drawn on the map as a plan-view arrow for the east–west component and a symbol for the vertical component. The arrow is drawn only if the east–west velocity is significant by (24) and $|v_E| \ge 0.05$ mm/yr, the displayed precision. It is then horizontal, points east if $v_E \gt 0$ and west otherwise, and has a screen length, independent of the zoom level, of

It is anchored at the probe location or, for a selection, at the mean longitude and latitude of the selected points. The vertical component is shown as a disc at the anchor, green with a central dot for uplift ($v_U \ge 0$, motion out of the drawing plane) and red with a cross for subsidence, following the engineering convention for vectors normal to the plane. A callout states $v_U \pm \sigma_{v_U}$ and $v_E \pm \sigma_{v_E}$, the totals (23), the period and the fact that the north–south component is not resolved. Otherwise no arrow is drawn: a neutral dashed ring around the vertical symbol and the callout note that no significant east–west motion was found, so that the minimum length of (25) never exaggerates a velocity that is indistinguishable from zero. The temporal coherence of each point is displayed with the heuristic classes good ($\ge 0.8$), fair (0.6–0.8) and poor ($\lt 0.6$).

5Averaging of multi-point selections

A single scatterer may move for reasons unrelated to the ground, for example structural or thermal deformation of a building. The multi-point tool therefore averages several points of each orbit before the decomposition. This section describes the selection, the averaging and the uncertainty assigned to the average.

5.1Selection

A click adds the nearest ascending and the nearest descending point, found as in Section 4.1 with radii of 30 m and then 100 m. A box drawn on the map adds all points of both orbits inside it from every candidate track; if an orbit contributes more than 100 points, the application keeps 100 of them by uniform striding, that is, the points with indices $\lfloor i\,m_{\mathrm{box}}/100 \rfloor$, $i = 0, \dots, 99$, where $m_{\mathrm{box}}$ is the number of points of that orbit in the box. The points of one orbit may therefore belong to different tracks, for instance A21-088 and A20-161.

5.2Averaging on a common grid

Let the selection contain $m$ points of a given orbit. The application constructs the grid (8) over the period covered by all $m$ series, or over the union of their periods if they do not overlap. Each series is interpolated by (9) and referred to its own robust level at its first grid epoch by (10). At each epoch the average is taken over the set $\mathcal M_k$ of members with a finite value,

and the LOS vector of the average is the arithmetic mean of the member unit vectors,

which is deliberately not renormalised. If all members experience the same displacement $\vD$, the average of their LOS displacements is $\tfrac1m\sum_i\los_i\cdot\vD = \boldsymbol\ell_o\cdot\vD$ exactly, so that (3) holds for the averaged series with $\boldsymbol\ell_o$ in place of $\los$. Normalising $\boldsymbol\ell_o$ would scale the decomposed displacement by $\lVert\boldsymbol\ell_o\rVert \le 1$, by 0.3 % for an equal mix of the two ascending tracks of Table 2, whose LOS vectors differ by 9.1°. The mean velocity and the coherence of the average are the arithmetic means of the member values. The averaged ascending and descending series then enter Section 4.3 unchanged.

5.3Uncertainty of an average

The application assigns to the average the root mean square of the member uncertainties,

both for the per-epoch noise and for the velocity uncertainty, and uses $\sigma_o$ as $\sigma\sa$ or $\sigma\sd$ in (18). The familiar reduction by $\sqrt m$ is deliberately not applied. The following proposition shows that (28) is conservative.

Equality holds when all members have the same $\sigma$ and fully correlated errors. For independent errors, (29) reduces to $\sigma_{\mathrm{mean}} = \sigma_o/\sqrt m$. The choice of the bound is motivated by the error structure of InSAR. The dominant error sources between nearby scatterers, namely the atmospheric phase screen, which is correlated over kilometres, and errors of the reference and the calibration, are common to all points of a small selection (Hanssen, 2001). For $m$ members with equal $\sigma$ and a common correlation $\rho$,

With $m = 15$ and $\rho = 0.8$, (30) gives $0.90\,\sigma$, whereas the independent model would give $\sigma/\sqrt{15} = 0.26\,\sigma$ and overstate the precision by a factor of 3.5. The RMS rule slightly overestimates the uncertainty of the uncorrelated part of the error, which we consider the appropriate side on which to err.

6Accuracy bands in the time-series chart

The chart shows the LOS series of the selected points, referred to their robust levels at the first epoch as in (10), the averaged series of Section 5 and the decomposed series, each with a shaded band intended to convey the scatter of individual epochs. When denoising is off, the band is centred on a running mean over a 30-day window,

and spans $\mu_k \pm \sigma$, where $\sigma$ is rmse_ts for a single point, $\sigma_o$ of (28) for an average and $\sigma_E$ or $\sigma_U$ of (18) for the decomposed series. Because the window is defined in time, it contains five acquisitions at 6-day sampling and three at 12-day sampling. When denoising is on, the band is centred on the denoised series $\hat x_k$ of Section 7 and spans $\hat x_k \pm \sigma_{\mathrm{res}}$, with $\sigma_{\mathrm{res}}$ from (40).

The band is a one-sigma band for individual epochs about a local mean. It is neither a confidence band for the local mean, which would be narrower by up to $\sqrt5$ for independent errors, nor a band for the displacement relative to the start, which would also contain the error of the reference level (Remark 2). The numerical readouts beside the chart show the EGMS mean LOS velocity with mean_velocity_std for each point and, after interpretation, $v_U \pm \sigma_{v_U}$ and $v_E \pm \sigma_{v_E}$ with the significance flag of (24).

7Temporal denoising

InSAR time series contain per-epoch noise of a few millimetres from atmospheric delay, processing and residual unwrapping errors, as well as seasonal motion of the scatterers themselves. To make slow deformation easier to see, the application offers an optional denoiser, shared by the chart and the displacement map, with two independent steps: removal of the annual cycle (Section 7.3) followed by Gaussian-weighted local-linear smoothing (Section 7.1). In the chart, denoising affects only the display; velocities and the decomposition are always computed from the raw series. In the displacement map, the denoised series replace the raw ones (Section 8.4). The same local-linear estimator, with a fixed scale of 60 days, provides the reference level of (10).

7.1Local-linear Gaussian smoothing on irregular sampling

The smoothed value at epoch $t_i$ is the intercept of a straight line fitted by weighted least squares to the available samples within three kernel scales,

with $J_i = \{\, j : |u_{ij}| \le 3,\ x_j \text{ finite} \,\}$ and $\sigma_t = 30$, 60 or 90 days for the levels Light, Medium and Strong; the default is Medium. This is local-linear kernel regression with a truncated Gaussian kernel (Fan and Gijbels, 1996). With the weighted moments $S_k = \sum_{j\in J_i} w_{ij}u_{ij}^k$ and $T_k = \sum_{j\in J_i} w_{ij}u_{ij}^k x_j$, the normal equations of (32) give the closed form

where $\ell_{ij}$ is the equivalent kernel. It satisfies $\sum_j \ell_{ij} = 1$ and $\sum_j \ell_{ij}u_{ij} = 0$ for any sampling, so constants and straight lines are reproduced exactly, including at the ends of a series and next to gaps. Where the samples are placed symmetrically about $t_i$, as in the interior of a regularly sampled series, $S_1 = 0$ and $\ell_{ij} = w_{ij}/S_0$: the estimator then coincides with the Gaussian kernel mean, the Nadaraya–Watson estimator (Nadaraya, 1964; Watson, 1964). If the local design is degenerate, that is, if $J_i$ contains fewer than three samples or the weighted spread of the sample times is below $0.2\,\sigma_t$ ($S_0S_2 - S_1^2 \lt 0.04\,S_0^2$), the application falls back to the kernel mean $T_0/S_0$, which avoids extrapolating a poorly determined slope. Because all weights are computed from the actual acquisition dates, the estimator handles the transition from 6-day to 12-day sampling without modification; it is also evaluated at epochs whose own value is missing, which fills such epochs whenever samples lie within $3\sigma_t$.

7.2Frequency response

In the interior of a regularly sampled series, (33) reduces to a convolution with the normalised Gaussian $g(\tau) = (2\pi\sigma_t^2)^{-1/2}\exp\bigl(-\tau^2/2\sigma_t^2\bigr)$. Its Fourier transform is again Gaussian,

so the smoother is a zero-phase low-pass filter. Setting $H(f_c) = 1/\sqrt2$, the half-power or −3 dB point, gives the cut-off frequency and period

Table 3 and Figure 5 give the values for the three levels. We verified numerically that the truncation at $\pm3\sigma_t$ and discrete sampling at 6 or 12 days change $|H|$ by less than 0.004 at all frequencies. An important consequence is that the smoother is not neutral with respect to seasonal signals: in the interior, Medium smoothing retains 59 % of the amplitude of an annual cycle and Strong smoothing 30 %. Removal of the annual cycle by Section 7.3 is the explicit alternative.

7.3Removal of the annual cycle

Thermal dilation of structures and moisture-driven soil motion produce an approximately annual cycle. The application models a series as a linear trend plus a sinusoid with a period of one year,

where $\bar t$ is the mean epoch of the samples with a finite value; centring the time variable makes the intercept and the slope nearly uncorrelated and improves the conditioning. With $\vg_i = \bigl(1, \tilde t_i, \sin 2\pi\tilde t_i, \cos 2\pi\tilde t_i\bigr)\T$, the least-squares estimate $\vbeta = (\beta_0, \dots, \beta_3)\T$ solves the normal equations

which the application solves by Gauss–Jordan elimination with partial pivoting. The fit requires at least eight samples and is abandoned, leaving the series unchanged, if a pivot falls below $10^{-12}$. Only the periodic part is removed,

and the amplitude $A$ is reported in the chart legend. Fitting the trend jointly with the sinusoid prevents the trend from leaking into the seasonal coefficients, which would otherwise happen whenever the series does not span an integer number of years. When both steps are active, (39) is applied before (32).

7.4Residual scatter

The band of a denoised series uses the scatter of the raw values about the denoised curve,

where $\mathcal K$ is the set of epochs at which both values are finite. This quantity measures what the filter removes. It underestimates the per-epoch noise, because the smoother partly follows the noise, and it includes any genuine high-frequency signal as well as the annual cycle when that has been removed.

8Displacement map

The displacement map visualises the vertical and east–west motion inside a user-drawn polygon as an animated height map with an arrow field. It is built from Level 3 rather than Level 2B data, because Level 3 already provides both components on a regular grid, whereas the use of Level 2B would require tens of thousands of time-series requests and a decomposition per cell. Figure 6 summarises the processing, which the following subsections describe in order.

8.1From polygon to grid

The user draws a polygon $\Pi$ in Web Mercator coordinates. Its bounding box may not exceed 25 km on either side, where widths are computed with 111.32 $\cos\varphi_c$ km per degree of longitude and 110.57 km per degree of latitude at the central latitude $\varphi_c$. Because the Mercator scale factor is $\sec\varphi$ in every direction, the ground area is estimated from the planar shoelace area $A_M$ of the polygon in map units as

and the number of Level 3 cells as $100\,A/\mathrm{km}^2$. The application then loads the vertical and east–west Level 3 cells inside the bounding box at full tile depth and keeps those inside $\Pi$ by the even–odd ray-casting test.

If a component has $M \gt 4000$ cells inside $\Pi$, a regular subset is retained. The cell centres are projected to EPSG:3035 (Lambert azimuthal equal-area with origin 52° N, 10° E, false easting 4 321 000 m and false northing 3 210 000 m on GRS 80), and integer lattice indices are computed as

and likewise $j_c$ from the northing. The offset $\phi_x$ is the circular mean of the fractional cell position and makes the indices insensitive to whether the coordinates refer to cell centres or corners. Cells with $i_c \equiv j_c \equiv 0 \pmod s$ are kept, where the stride $s$ starts at $\max\bigl(2, \lceil\sqrt{M/4000}\,\rceil\bigr)$ and is increased until at most 4000 cells remain or the stride exceeds 64. The result is an exactly regular subgrid with a spacing of $100\,s$ m.

For rendering, each retained cell is mapped to one pixel of a raster in Web Mercator. The pixel size is $s$ times the median nearest-neighbour distance of up to 64 sampled cells, accepted if it lies between 0.6 and 1.6 times the nominal $100\ \mathrm m/\cos\varphi_c$ and replaced by the nominal value otherwise. Each pixel takes the cell nearest to its centre within 0.75 pixel; remaining empty pixels take the nearest cell within 1.5 pixels, and enclosed single-pixel holes with at least three of four neighbours filled take the nearest cell within 2.2 pixels. These tolerances compensate for the slight rotation of the EPSG:3035 lattice relative to Web Mercator, while genuine gaps in the data, such as heathland or fields without scatterers, remain empty. Pixels whose centres lie outside $\Pi$ are masked.

8.2The time cube

The time series of the retained cells are requested in batches of 150 cells, the vertical component first. Each component is stored as a time cube $C \in \mathbb R^{F\times M}$ on a regular grid of frames,

where $t_s$ and $t_e$ are the earliest and the latest epoch among all series loaded so far, and $F = \lfloor (t_e - t_s)/\Delta t \rfloor + 1$, increased by one if $t_{F-1} \lt t_e$. When a batch contains an epoch outside $[t_s, t_e]$, the grid is extended to the union of the old grid and the new span, a new cube is allocated and all series loaded before are resampled onto it; the original series are kept for this purpose. No epoch of any loaded series is therefore dropped. Each series $(t_{c,j}, x_{c,j})_{j=1}^{N_c}$ is resampled onto the frames by linear interpolation (9) inside its period and by holding its end values outside it,

Since all Level 3 cells of the current release share one 6-day lattice (Section 2.2), the frames coincide with Level 3 epochs and (44) reduces to a copy within each series; holding applies only to the first frames of cells whose series begin after $t_s$, at most two frames in the current release. Cells for which no series is available keep the velocity model of (46).

The map displays the displacement at a playhead time $t$ relative to the start $t_0$ of the selected period. The playhead is converted to a fractional frame index, and the cube is interpolated linearly between the neighbouring frames (Figure 7),

The start is represented, as in (10), by the robust level $L_c(t_0)$ of the cell: the local-linear estimate (33) with $\sigma_{\mathrm{ref}} = 60$ days, evaluated on the frames of the displayed cube within $\pm 3\sigma_{\mathrm{ref}}$ of $t_0$. The displacement is then

where $v_c$ is the mean velocity of the cell from the tile. The second case provides an immediate preview before the series arrive and is replaced cell by cell as they are loaded. The levels $L_c(t_0)$ depend only on $t_0$ and on the cube contents; they are computed once when either changes and reused for every rendered frame, so the robust reference does not slow down playback. For independent noise and 6-day frames, the reference contributes $0.41\,\sigma$ of noise at the first frame and $0.17\,\sigma$ in the interior (Remark 2) instead of the full $\sigma$ of a single frame.

8.3Regional trend removal

By default the application removes, at every frame, the median displacement of all $M$ cells of the component,

where the median of an even number of values is the mean of the two central ones. The operation removes any signal common to the polygon, such as regional subsidence or uplift, a residual offset of the reference frame of the calibrated product, or a residual atmospheric contribution shared by the area, so that differential motion stands out. The median is used instead of the mean because its breakdown point is 50 %: a local subsidence bowl that occupies less than half of the cells does not bias the reference. Conversely, if a deformation feature covers more than half of the polygon, the median lies inside the feature and the remaining area appears to move in the opposite direction. The removed value is displayed so that the user can judge this. With the option enabled, all map quantities are relative to the polygon; with the option disabled, $D'_c \equiv D_c$ in what follows.

8.4Temporal denoising of the cube

When the denoiser is active, (39) and (33) are applied to every cell with a loaded series, along the frame axis, with the frame epochs (43) as sampling times and the same settings as in the chart. Because the frames are regular, the Gaussian weights depend only on the frame offset; they are tabulated once per run, and the per-cell computation allocates no memory and evaluates no exponential, which keeps the cost at about 40 ms per 1000 cells. The computation runs in a background worker and its result replaces $C$ in (45)–(47), in the colour range and in the gradient model. The cells are processed independently, without spatial smoothing; an optional $3\times3$ median filter of the rendered field is described in Section 8.5. Because the local-linear estimator reproduces linear trends at both ends of the cube (Remark 3), smoothing does not reduce the displacement between the first and the last frame.

8.5Height map

The vertical field $D'_U$ is coloured with a diverging red–white–blue palette, red for subsidence and blue for uplift, through a lookup table of 256 entries interpolated between nine colour stops. For a symmetric range $\pm R$, the table index is

Unless the user sets $R$ manually, it is the 98th percentile of the absolute displacement sampled at sixteen evenly spaced times of the selected period $[t_0, t_1]$,