PlanX DataCube Lab

Theoretical Background

The spatiotemporal data cube is a conceptual and computational structure that organises geospatial information along three orthogonal axes: two spatial dimensions ($X$, $Y$) and a temporal axis ($T$), yielding a hypercube of observations with an optional variable dimension ($V$). The approach formalises what Hägerstrand (1970) termed the space-time prism—a bounded region in space and time within which phenomena can be observed. Adopting a regular tessellation of space into a fishnet (lattice) grid, the data cube imposes a common spatial reference that makes observations from heterogeneous sources commensurable.

In GIScience, the data cube paradigm addresses a fundamental challenge: vector layers acquired at different epochs rarely share a common geometry, making direct temporal comparison impossible without spatial harmonisation. Rasterising each time step onto an identical grid resolves this mismatch. The resulting hypercube supports array-oriented analytics (slicing, dicing, aggregation along arbitrary axes) and interoperates with scientific data formats such as netCDF4 and Zarr, which are designed for chunked, compressed, parallel I/O on multidimensional arrays.

The algorithm accepts either vector layers (rasterised by field value) or pre-existing raster stacks. When multiple input layers are supplied, each is assigned to a discrete time index, producing a cube of dimensionality $|X| \times |Y| \times |T| \times |V|$ where $|X|$ and $|Y|$ denote the number of columns and rows of the fishnet grid, $|T|$ the number of time steps, and $|V|$ the number of variables. The design follows the Earth System Data Cube specification (Mahecha et al., 2020) and the Open Geospatial Consortium (OGC) Coverage Implementation Schema.

Mathematical Formulation

Grid cell assignment. Given a point $(x, y)$ and a regular grid with origin $(x_{\min}, y_{\min})$ and cell resolution $\Delta_g$, the column and row indices are:

Fishnet polygon geometry. For cell $(c, r)$, the axis-aligned rectangle is:

Rasterisation. For vector input, each polygon is scan-converted onto the grid. The value assigned to cell $(c, r)$ at time $t$ is the area-weighted mean of all polygons intersecting the cell (or, at the user's choice, the modal value):

where $P_t$ is the set of polygons at time $t$, $a_{p \cap (c,r)}$ is the area of intersection between polygon $p$ and grid cell $(c, r)$, and $f_p$ is the value field of polygon $p$.

Cube dimensionality. The resulting data cube has the shape:

where $|X| = \lceil (x_{\max} - x_{\min}) / \Delta_g \rceil$ and $|Y| = \lceil (y_{\max} - y_{\min}) / \Delta_g \rceil$.

netCDF4 / Zarr serialisation. The cube array $\mathbf{C} \in \mathbb{R}^{|X| \times |Y| \times |T|}$ is written with CF-1.8 conventions, including coordinate variables (lon, lat, time) and grid-mapping metadata.

Parameters

Output Description

Interpretation Guide

The grid cell size is the primary determinant of information density. Choose a cell size that is at least twice the typical feature size of the input data to avoid aliasing, and no larger than the coarsest acceptable analytical unit. A useful heuristic is to set $\Delta_g$ to the mean nearest-neighbour distance of input features divided by two.

netCDF4 output is recommended for single-file portability and compatibility with Panoply, QGIS, and Python's xarray. Zarr output is recommended for cloud-native workflows where parallel, chunked reads are needed. Both formats preserve the full precision of the rasterised values; no rounding or quantisation is performed unless the input type requires it.

The optional grid output layer is useful for inspecting which cells are populated and for linking subsequent panel-format analyses back to the original grid. The cell_id column is the join key for all downstream spatiotemporal analyses.

References

1. Hägerstrand, T. (1970). What about people in regional science? Papers of the Regional Science Association, 24(1), 7–21. DOI: 10.1007/BF01936872
2. Andrienko, G., Andrienko, N., & Gatalsky, P. (2003). Exploratory analysis of spatial and temporal data: A systematic approach. Springer Science & Business Media.
3. Mahecha, M. D., Gans, F., Brandt, G., et al. (2020). Earth system data cubes unravel global multivariate dynamics. Earth System Dynamics, 11(1), 201–234. DOI: 10.5194/esd-11-201-2020
4. Pebesma, E. (2012). spacetime: Spatio-temporal data in R. Journal of Statistical Software, 51(7), 1–30. DOI: 10.18637/jss.v051.i07
5. Baumann, P., Hirschorn, E., & Maso, J. (2017). OGC Coverage Implementation Schema. Open Geospatial Consortium, Document 09-146r8.
6. ESRI (2016). ArcGIS Pro: How Create Space Time Cube Works. Environmental Systems Research Institute, Redlands, CA.
7. Eaton, B., Gregory, J., Drach, B., et al. (2022). NetCDF Climate and Forecast (CF) Metadata Conventions, Version 1.10. cfconventions.org

Theoretical Background

Space-time aggregation transforms a collection of discrete point events into a structured panel dataset by simultaneously binning observations into spatial cells and temporal intervals. This operation generalises the classical quadrat-count method of spatial point-pattern analysis (Diggle, 2013) into the temporal domain, producing a balanced spatiotemporal lattice that serves as the input for every subsequent analytical algorithm in the plugin.

The aggregation step addresses two core needs. First, it converts irregularly sampled point data into a regular structure amenable to array-based computation. Second, it provides a natural framework for summarising event intensity within each space-time bin through user-selectable statistics (count, sum, mean, maximum, minimum). The choice of aggregation mode is consequential: count captures frequency of occurrence, sum captures cumulative magnitude (when a value field is supplied), and the central-tendency measures capture typical magnitude per bin.

The output is a long-format panel layer where each row corresponds to a single cell-time combination, identified by a composite key $(c, r, t)$. This structure is directly compatible with the fixed-effects panel regression tradition in econometrics (Wooldridge, 2010) and the lattice-based spatiotemporal models in climatology (Cressie & Wikle, 2011). The cell_id column, computed as $\text{row} \times n_{\text{cols}} + \text{col}$, provides a flat integer key for efficient joins.

Mathematical Formulation

Temporal bin assignment. For an event with timestamp $t$ and a time-step size $\Delta_t$, the temporal bin index is:

Spatial binning. The spatial bin is determined by the same grid-cell assignment as Algorithm 1 (Equations 1.1–1.2). Each event is assigned to the unique space-time tuple $(\text{col}, \text{row}, t_{\text{bin}})$.

Aggregation statistics. For each space-time bin $b = (c, r, \tau)$, the accumulated statistics are computed incrementally from the set of events $E_b$:

Panel cell identifier. The flat integer key for each spatial cell is:

where $n_{\text{cols}} = |X|$ is the number of grid columns.

Drop-empty-cell filter. When enabled, any row for which $n_b = 0$ is excluded from the output. This reduces memory footprint and avoids propagating structural zeros into downstream statistical tests where they could bias variance estimates.

Long-format panel. The output matrix $\mathbf{P}$ has $N$ rows and columns $[\text{cell\_id}, \text{col}, \text{row}, \text{time\_step}, \text{time\_value}, \text{value}, \text{count}]$, where $N \leq |X| \cdot |Y| \cdot |T|$.

Parameters

Output Description

Interpretation Guide

The time-step size should be large enough that typical bins contain at least 5–10 events for stable statistics, yet small enough to resolve the temporal dynamics of interest. A useful diagnostic is to inspect the histogram of bin counts after aggregation: a heavy left tail (many bins with 0–1 events) indicates that the time step is too fine or that events are spatially sparse.

The Modifiable Areal Unit Problem (MAUP; Openshaw, 1983) applies to both the spatial grid and the temporal bin. Results should be interpreted with awareness that alternative grid resolutions or temporal partitions may yield different patterns. Sensitivity analysis across a reasonable range of cell sizes and time steps is recommended before drawing substantive conclusions.

When a value_field is supplied, the count field still records the number of contributing events. Ratios (e.g., sum/count = mean) are arithmetic identities, so reporting only the aggregated value without the count is discouraged—the count serves as a measure of reliability.

References

1. Peuquet, D. J. (1994). It's about time: A conceptual framework for the representation of temporal dynamics in geographic information systems. Annals of the Association of American Geographers, 84(3), 441–461. DOI: 10.1111/j.1467-8306.1994.tb01873.x
2. Langran, G. (1992). Time in Geographic Information Systems. Taylor & Francis, London.
3. Openshaw, S. (1983). The modifiable areal unit problem. Concepts and Techniques in Modern Geography, 38. Geo Books, Norwich.
4. Diggle, P. J. (2013). Statistical Analysis of Spatial and Spatio-Temporal Point Patterns, 3rd ed. Chapman & Hall/CRC.
5. Cressie, N. & Wikle, C. K. (2011). Statistics for Spatio-Temporal Data. John Wiley & Sons.
6. Wooldridge, J. M. (2010). Econometric Analysis of Cross Section and Panel Data, 2nd ed. MIT Press.
7. Worboys, M. F. (1994). A unified model for spatial and temporal information. The Computer Journal, 37(1), 26–34. DOI: 10.1093/comjnl/37.1.26

Theoretical Background

Emerging Hot Spot Analysis (EHSA) is a two-stage procedure that combines the Getis-Ord $G_i^*$ local spatial autocorrelation statistic with the nonparametric Mann-Kendall trend test to classify each spatial unit into one of seventeen spatiotemporal pattern categories. The method was popularised by ESRI's Space Time Pattern Mining toolbox and has been applied extensively in criminology, epidemiology, and urban analytics to detect locations where event intensity is not merely high or low but evolving in statistically meaningful ways over time.

The first stage computes $G_i^*$ independently at each time step. The $G_i^*$ statistic (Getis & Ord, 1992; Ord & Getis, 1995) identifies local clusters of high values (hot spots) and low values (cold spots) by comparing the weighted local mean to the global mean, standardised by the sample standard deviation. The spatial weights matrix $\mathbf{W}$ may be defined by a fixed distance band (binary adjacency) or by a $K$-nearest-neighbour scheme. The second stage applies the Mann-Kendall trend test (Mann, 1945; Kendall, 1975) to the time series of $G_i^*$ $z$-scores at each location, detecting whether the intensity of spatial clustering is increasing, decreasing, or stable.

The seventeen output classes are derived from a decision tree that considers the significance pattern of the $G_i^*$ $z$-scores across time, the direction and significance of the Mann-Kendall trend, the presence of consecutive runs of significant hot or cold spots, and oscillation between hot and cold states. This taxonomy captures nuanced dynamics such as intensifying, diminishing, sporadic, and oscillating hot and cold spots, providing a richer characterisation than a single cross-sectional hot-spot map.

Mathematical Formulation

Getis-Ord $G_i^*$ statistic. For each time step $t$ and location $i$, let $x_{i,t}$ be the value and let $w_{ij}$ be the spatial weight between locations $i$ and $j$. The statistic is:

where $\bar{X}_t = \frac{1}{n}\sum_{j=1}^{n} x_{j,t}$ and $S_t = \sqrt{\frac{1}{n}\sum_{j=1}^{n} x_{j,t}^2 - \bar{X}_t^2}$.

Kendall's $S$ statistic. For the time series of $G_i^*$ $z$-scores $\mathbf{z}_i = (z_{i,1}, \dots, z_{i,T})$, the Mann-Kendall test statistic is:

where $\operatorname{sgn}(\theta) = 1$ if $\theta > 0$, $-1$ if $\theta < 0$, and $0$ otherwise.

Variance of $S$. Under the null hypothesis of no trend, the variance incorporates a correction for tied values:

where $m$ is the number of tied groups and $t_g$ is the size of the $g$-th tied group.

Theil-Sen slope estimator. The trend magnitude is the median of all pairwise slopes:

17-class decision tree. Each location is classified by traversing the following rule hierarchy: (i) determine whether the location is ever a significant hot spot or cold spot; (ii) apply the Mann-Kendall test to the $z$-score series; (iii) count consecutive significant steps at the end of the series; (iv) detect oscillation if the location alternates between hot and cold significance. The resulting 17 EHSA codes (0–16) map to labels: No pattern detected, New hot spot, Consecutive hot spot, Intensifying hot spot, Persistent hot spot, Diminishing hot spot, Sporadic hot spot, Oscillating hot spot, Historical hot spot, and their cold-spot analogues, plus the New (most recent step only) variants. The exact decision tree is documented in ESRI (2016).

Parameters

Output Description

Interpretation Guide

The EHSA classification is inherently multi-dimensional, and no single class is universally "concerning" or "desirable." A persistent hot spot indicates sustained high intensity over the entire study period, which may represent chronic conditions (e.g., endemic disease areas, persistent crime corridors). An intensifying hot spot indicates that intensity is not only high but statistically increasing, warranting closer monitoring or intervention.

The choice of spatial relationship is critical. A distance band assumes isotropy and uniform interaction range, suitable for phenomena with a known spatial range (e.g., pollutant dispersion). KNN adapts to variable density but can mask true spatial structure in sparse regions. Conduct sensitivity analysis across multiple distance thresholds or $K$ values to assess stability.

The Mann-Kendall test assumes serially independent observations. Temporal autocorrelation in the $G_i^*$ $z$-scores can inflate the Type I error rate (Yue et al., 2002). When $T$ is small ($T < 10$), the test has low power; when $T$ is large, consider pre-whitening or block bootstrapping to control for autocorrelation.

References

1. Getis, A. & Ord, J. K. (1992). The analysis of spatial association by use of distance statistics. Geographical Analysis, 24(3), 189–206. DOI: 10.1111/j.1538-4632.1992.tb00261.x
2. Ord, J. K. & Getis, A. (1995). Local spatial autocorrelation statistics: Distributional issues and an application. Geographical Analysis, 27(4), 286–306. DOI: 10.1111/j.1538-4632.1995.tb00912.x
3. Anselin, L. (1995). Local indicators of spatial association—LISA. Geographical Analysis, 27(2), 93–115. DOI: 10.1111/j.1538-4632.1995.tb00338.x
4. Mann, H. B. (1945). Nonparametric tests against trend. Econometrica, 13(3), 245–259. DOI: 10.2307/1907187
5. Kendall, M. G. (1975). Rank Correlation Methods, 4th ed. Charles Griffin, London.
6. Sen, P. K. (1968). Estimates of the regression coefficient based on Kendall's tau. Journal of the American Statistical Association, 63(324), 1379–1389. DOI: 10.1080/01621459.1968.10480934
7. ESRI (2016). ArcGIS Pro: How Emerging Hot Spot Analysis Works. Environmental Systems Research Institute, Redlands, CA.
8. Yue, S., Pilon, P., Phinney, B., & Cavadias, G. (2002). The influence of autocorrelation on the ability to detect trend in hydrological series. Hydrological Processes, 16(9), 1807–1829. DOI: 10.1002/hyp.1095

Theoretical Background

Local Moran's $I$, introduced by Anselin (1995) as a Local Indicator of Spatial Association (LISA), decomposes the global Moran's $I$ statistic into location-specific components that identify spatial clusters and spatial outliers. Where global Moran's $I$ provides a single summary of spatial autocorrelation, the local variant produces a map of statistically significant clusters of high values (HH), clusters of low values (LL), high-value outliers surrounded by low values (HL), and low-value outliers surrounded by high values (LH). This decomposition is essential for understanding the spatial heterogeneity of a phenomenon.

In the spatiotemporal context of DataCube Lab, the algorithm first aggregates the panel data to a single value per location using a user-selected temporal reduction (mean, last observation, sum, or maximum). This produces a cross-sectional dataset to which the local Moran's $I$ is applied. The temporal aggregation step is not a compromise but a deliberate analytical choice: it collapses the temporal dimension to reveal the average spatial structure, which can then be compared with the trend results from Algorithm 5 and the anomaly results from Algorithm 6 to construct a multi-faceted characterisation of each location.

Significance is assessed via a conditional permutation test (Anselin, 1995). For each location, the observed value $x_i$ is held fixed while the values of the neighbours are randomly drawn (without replacement) from the remaining $n-1$ observations. A pseudo $p$-value is computed as the proportion of permutations in which the absolute local $I$ statistic equals or exceeds the observed value. This nonparametric approach avoids distributional assumptions and is robust to non-normality.

Mathematical Formulation

Local Moran's $I$. For location $i$ with deviation from the mean $d_i = x_i - \bar{X}$, the local statistic is:

Spatial lag. The spatially lagged deviation at location $i$ is:

The spatial lag is positive when neighbours tend to be above the global mean and negative when they tend to be below it.

Conditional permutation test. For $N_{\text{perm}}$ permutations at location $i$, let $I_i^{(k)}$ be the statistic computed using a random draw of the neighbour set from $\{x_j : j \neq i\}$. The pseudo $p$-value is:

The unit increment in both numerator and denominator ensures the $p$-value is never zero (a conservative adjustment; see North et al., 2002).

$z$-score standardisation. The observed $I_i$ is converted to a $z$-score using the permutation mean $\mu_{\text{perm}}$ and standard deviation $\sigma_{\text{perm}}$:

Quadrant classification. Each location is assigned to one of five categories based on the sign of the standardised deviation $z_i = \operatorname{sign}(d_i) \cdot |z_i|$ and the spatial lag:

Parameters

Output Description

Interpretation Guide

HH clusters indicate spatial agglomeration of high values—locations where both the unit and its neighbours exhibit above-average intensity. These are conventionally interpreted as "hot spots" in the LISA framework (distinct from the $G_i^*$ hot spots of Algorithm 3, which use a different null hypothesis). LL clusters indicate low-value agglomerations.

HL outliers (high surrounded by low) are especially informative: they identify "islands" of high intensity in an otherwise low-intensity region, suggesting localised anomalies that may be driven by site-specific conditions rather than regional processes. LH outliers (low surrounded by high) indicate under-performing locations in a generally high-intensity area. Both outlier types warrant case-level investigation.

The permutation count should be at least 999 for publication-quality $p$-values (the minimum resolvable $p$-value is $1/(N_{\text{perm}}+1)$). Set permutations to 0 to disable the test and obtain only the point estimates. The temporal aggregation method affects the interpretation: mean smooths transient fluctuations, last captures the most recent state, and sum reflects cumulative magnitude over the study period.

References

1. Anselin, L. (1995). Local indicators of spatial association—LISA. Geographical Analysis, 27(2), 93–115. DOI: 10.1111/j.1538-4632.1995.tb00338.x
2. Anselin, L. (2019). A local indicator of multivariate spatial association: Extending Geary's $c$. Geographical Analysis, 51(2), 133–150.
3. Cliff, A. D. & Ord, J. K. (1981). Spatial Processes: Models and Applications. Pion, London.
4. North, B. V., Curtis, D., & Sham, P. C. (2002). A note on the calculation of empirical $p$-values from Monte Carlo procedures. American Journal of Human Genetics, 71(2), 439–441. DOI: 10.1086/341527
5. Sokal, R. R., Oden, N. L., & Thomson, B. A. (1998). Local spatial autocorrelation in a biological model. Geographical Analysis, 30(4), 331–354.
6. Tiefelsdorf, M. (2002). The saddlepoint approximation of Moran's $I$'s and local Moran's $I_i$'s reference distributions and their numerical evaluation. Geographical Analysis, 34(3), 187–206.

Theoretical Background

Temporal trend analysis addresses the simplest and most frequently asked question in spatiotemporal data analysis: for each location, is the observed time series increasing, decreasing, or stable? The algorithm employs the Mann-Kendall nonparametric trend test (Mann, 1945; Kendall, 1975) paired with the Theil-Sen slope estimator (Theil, 1950; Sen, 1968). This combination is the de facto standard for environmental monitoring, recommended by the World Meteorological Organization for trend detection in hydro-climatic time series.

The Mann-Kendall test is a rank-based procedure that tests the null hypothesis $H_0$ of no monotonic trend against the two-sided alternative $H_a$ of an increasing or decreasing monotonic trend. Its nonparametric nature confers two advantages: it does not require normally distributed data, and it is robust to outliers, which are common in geospatial time series (e.g., extreme weather events, sudden land-use changes). The test statistic $S$ counts the number of concordant minus discordant pairs across all possible pairwise comparisons.

The Theil-Sen slope provides a robust estimate of the trend magnitude. Unlike ordinary least squares, which can be severely distorted by a single outlier, the Theil-Sen estimator is the median of all pairwise slopes and has a breakdown point of approximately 29.3%, meaning that nearly 30% of the data can be contaminated before the estimator becomes unreliable (Rousseeuw & Leroy, 1987). The Sen intercept (median of residuals) completes the trend line, making it suitable for generating predicted values and visualisations.

Mathematical Formulation

Kendall's $\tau$ correlation coefficient. For a time series of length $n$:

where $\tau \in [-1, 1]$, with $+1$ indicating perfect increasing monotonicity and $-1$ perfect decreasing monotonicity.

Standardised Mann-Kendall $Z$ statistic.

The continuity correction ($\pm 1$) improves the normal approximation for small $n$.

Theil-Sen slope. The robust trend magnitude per time step:

Sen intercept. The intercept of the trend line, estimated as:

$p$-value via normal approximation. The two-sided $p$-value is obtained from the standard normal cumulative distribution function $\Phi$:

The implementation uses the Abramowitz & Stegun (1964) rational approximation for $\Phi$, accurate to $7.5 \times 10^{-8}$.

Parameters

Output Description

Interpretation Guide

Kendall's $\tau$ measures the strength of monotonic association, not the rate of change. A $\tau$ of 0.6 with a Sen slope of 0.01 indicates a highly consistent but very slow increase; a $\tau$ of 0.3 with a Sen slope of 5.0 indicates a weaker consistency but a rapid magnitude of change. Report both statistics. Use the Sen slope for substantive interpretation (e.g., "the value increases by $\beta$ units per time step").

The Mann-Kendall test has low power for $n < 8$—trends may be visually apparent but fail to reach significance. Conversely, for $n > 30$, even trivially small trends may achieve statistical significance. Contextualise the results with the Sen slope to distinguish statistical significance from practical importance.

Temporal autocorrelation inflates the Type I error rate (Yue et al., 2002). If the time series exhibit strong lag-1 autocorrelation (common in monthly data), consider applying the trend-free pre-whitening (TFPW) procedure before analysis. The algorithm reports the raw (unadjusted) $p$-value; users working with serially correlated data should note this caveat.

References

1. Mann, H. B. (1945). Nonparametric tests against trend. Econometrica, 13(3), 245–259. DOI: 10.2307/1907187
2. Kendall, M. G. (1975). Rank Correlation Methods, 4th ed. Charles Griffin, London.
3. Sen, P. K. (1968). Estimates of the regression coefficient based on Kendall's tau. Journal of the American Statistical Association, 63(324), 1379–1389. DOI: 10.1080/01621459.1968.10480934
4. Theil, H. (1950). A rank-invariant method of linear and polynomial regression analysis, I, II, III. Proceedings of the Koninklijke Nederlandse Akademie van Wetenschappen, 53, 386–392, 521–525, 1397–1412.
5. Yue, S., Pilon, P., Phinney, B., & Cavadias, G. (2002). The influence of autocorrelation on the ability to detect trend in hydrological series. Hydrological Processes, 16(9), 1807–1829. DOI: 10.1002/hyp.1095
6. Abramowitz, M. & Stegun, I. A. (1964). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. National Bureau of Standards, Washington, DC.
7. Rousseeuw, P. J. & Leroy, A. M. (1987). Robust Regression and Outlier Detection. John Wiley & Sons, New York.
8. Helsel, D. R. & Hirsch, R. M. (2002). Statistical Methods in Water Resources. Techniques of Water-Resources Investigations, Book 4, Chapter A3. U.S. Geological Survey.

Theoretical Background

Anomaly detection in spatiotemporal data identifies individual time points at each location that deviate markedly from the location's typical behaviour. Unlike trend analysis, which characterises the overall direction of change, anomaly detection pinpoints when unusual values occurred, making it valuable for event detection (e.g., disease outbreaks, traffic incidents, ecological disturbances) and data-quality screening (identifying sensor malfunctions or recording errors).

The algorithm offers two complementary approaches. The robust modified $z$-score, based on the median and Median Absolute Deviation (MAD), is recommended for most geospatial applications because the median and MAD have breakdown points of 50%, rendering them immune to the influence of up to half the data being anomalous (Iglewicz & Hoaglin, 1993). The scaling factor $0.6745$ calibrates the MAD so that the robust $z$-scores are comparable to standard $z$-scores for normally distributed data. The standard $z$-score ($z = (x - \mu)/\sigma$) is also available for users who prefer a parametric approach, but it should be used with caution because both the mean and standard deviation are themselves distorted by the anomalies they aim to detect (Leys et al., 2013; Hampel, 1974).

The output is a per-location summary that includes the number and rate of anomalies, the maximum absolute $z$-score, the $z$-score of the last observation (useful for real-time monitoring), and the timestamps of the first and last detected anomalies. The last_anomaly flag (0/1) answers the operational question: "Is the most recent observation anomalous?"

Mathematical Formulation

Robust modified $z$-score. For location $i$ with observations $\mathbf{x}_i = (x_{i,1}, \dots, x_{i,n})$:

The constant $0.6745 = \Phi^{-1}(0.75)$ is the 75th percentile of the standard normal distribution, chosen so that $\operatorname{E}[0.6745 \cdot \operatorname{MAD}] = \sigma$ when $X \sim \mathcal{N}(\mu, \sigma^2)$.

Standard $z$-score. The parametric alternative:

Anomaly flag. An observation is flagged as anomalous if its absolute $z$-score exceeds a user-defined threshold $\theta$:

Recommended thresholds are $\theta = 3.5$ for the robust $z$-score (Iglewicz & Hoaglin, 1993) and $\theta = 2.0$ or $2.5$ for the standard $z$-score.

Anomaly rate. The proportion of anomalous time points at location $i$:

First and last anomaly time. If anomalies exist at location $i$, then $t_{\text{first}} = \min\{t : \text{anomaly}_{i,t} = 1\}$ and $t_{\text{last}} = \max\{t : \text{anomaly}_{i,t} = 1\}$.

Parameters

Output Description

Interpretation Guide

A high anomaly rate (say, $r > 0.15$) does not necessarily indicate a problematic location; it may reflect a genuinely volatile process or a time series that is poorly modelled by a constant central tendency. Conversely, a location with only one anomaly but a very high max_abs_z deserves scrutiny—it may represent a one-off extreme event or a data-entry error.

The last_anomaly flag is designed for operational monitoring workflows: run the algorithm on updated data, filter for last_anomaly = 1, and dispatch field verification or mitigation measures to those locations. The last_z field provides the signed magnitude (positive for anomalously high, negative for anomalously low).

When using the robust method, be aware that the MAD can be exactly zero if more than half the observations share the same value (common in count data with many zeros). In this case, the algorithm falls back to the standard $z$-score with a warning in the log. The standard method is appropriate only when the underlying distribution is approximately symmetric and uncontaminated; if in doubt, use the robust method.

References

1. Iglewicz, B. & Hoaglin, D. C. (1993). How to Detect and Handle Outliers. ASQC Quality Press, Milwaukee, WI.
2. Hampel, F. R. (1974). The influence curve and its role in robust estimation. Journal of the American Statistical Association, 69(346), 383–393. DOI: 10.1080/01621459.1974.10482962
3. Leys, C., Ley, C., Klein, O., Bernard, P., & Licata, L. (2013). Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median. Journal of Experimental Social Psychology, 49(4), 764–766. DOI: 10.1016/j.jesp.2013.03.013
4. Rousseeuw, P. J. & Croux, C. (1993). Alternatives to the median absolute deviation. Journal of the American Statistical Association, 88(424), 1273–1283. DOI: 10.1080/01621459.1993.10476408
5. Hawkins, D. M. (1980). Identification of Outliers. Chapman & Hall, London.
6. Barnett, V. & Lewis, T. (1994). Outliers in Statistical Data, 3rd ed. John Wiley & Sons.
7. Chandola, V., Banerjee, A., & Kumar, V. (2009). Anomaly detection: A survey. ACM Computing Surveys, 41(3), Article 15. DOI: 10.1145/1541880.1541882

Theoretical Background

Forecast backtesting (also known as time-series cross-validation with a fixed-origin hold-out) is the gold-standard procedure for selecting the most appropriate forecasting model for a given time series. Rather than committing to a single model class a priori, the algorithm fits three structurally distinct models—Exponential Smoothing (ETS), ARIMA, and Random Forest (RF)—on a training portion of each location's time series and evaluates their predictive accuracy on a held-out test set. This empirical approach operationalises the principle that "the data choose the model," which is especially valuable in spatial contexts where the optimal model class may vary across locations due to heterogeneous data-generating processes.

The three model classes represent different structural assumptions. Exponential Smoothing (Hyndman & Athanasopoulos, 2018) captures locally adaptive level via an exponentially weighted moving average; it is parsimonious, interpretable, and often performs well for short, noisy series. ARIMA (Box & Jenkins, 1970) models the autocorrelation structure through autoregressive and moving-average terms after optional differencing, making it suitable for series with persistent memory. Random Forest (Breiman, 2001) constructs an ensemble of regression trees on lagged features, capturing nonlinear interactions and threshold effects that linear models miss, but at the cost of greater data hunger and reduced interpretability.

Model selection is based on Root Mean Squared Error (RMSE) on the test set, with Mean Absolute Error (MAE) and Mean Absolute Percentage Error (MAPE) reported as supplementary diagnostics. RMSE penalises large errors quadratically, making it sensitive to outlier forecasts, while MAE provides a scale-dependent measure in the original units. MAPE expresses error as a percentage of the observed value, facilitating comparison across locations with different scales, but is undefined when any observed value is zero. The hold-out length is user-specified, typically 10–30% of the series length, balancing the need for a reliable training fit against a meaningful out-of-sample evaluation.

Mathematical Formulation

Exponential Smoothing (simple). The one-step-ahead forecast at time $t+1$ given observations up to time $t$ is the exponentially weighted moving average:

The smoothing parameter $\alpha \in (0, 1)$ is estimated by minimising the sum of squared one-step-ahead forecast errors on the training set. The forecast for all $h$ steps ahead is held constant at the last fitted level: $\hat{x}_{T+h} = F_{T+1}$.

ARIMA model. The ARIMA$(p, d, q)$ model is:

where $B$ is the backshift operator ($B x_t = x_{t-1}$), $\phi(B) = 1 - \phi_1 B - \cdots - \phi_p B^p$ is the autoregressive polynomial, $\theta(B) = 1 + \theta_1 B + \cdots + \theta_q B^q$ is the moving-average polynomial, $d$ is the order of differencing, and $\varepsilon_t \sim \text{WN}(0, \sigma^2)$ is white noise. For backtesting, a fixed ARIMA(1,0,1) order is used for computational tractability across potentially thousands of locations.

Random Forest recursive forecast. A feature matrix is constructed from lagged values of the training series using a sliding window of length $w$. For each window position $p$ from $w$ to $T_{\text{train}}$:

An ensemble of $B$ regression trees is trained on $\{(\mathbf{X}_p, y_p)\}$. Forecasts are generated recursively: $\hat{x}_{t+1} = f_{\text{RF}}([x_{t-w+1}, \dots, x_t])$, and this predicted value is fed back into the window for subsequent steps.

Forecast accuracy metrics. Let $h$ be the hold-out length and $y_t, \hat{y}_t$ the observed and predicted values for $t = 1, \dots, h$:

MAPE is computed only over $t$ for which $y_t \neq 0$; if all $y_t = 0$, MAPE is reported as null.

Best method selection. For location $i$, the method with the lowest RMSE is selected:

Parameters

Output Description

Interpretation Guide

The best_method field indicates which model class the data support for each location. A heterogeneous map (different methods dominating in different regions) suggests spatially varying data-generating processes and justifies the backtesting approach over a one-size-fits-all model. If one method dominates uniformly, that method may be preferred for forecasting (Algorithm 8) with stronger theoretical grounding.

RMSE values should be interpreted relative to the scale of the data. A useful normalisation is $\text{RMSE} / \sigma_{\text{train}}$, where values below 0.5 indicate good predictive performance and values above 1.0 indicate that the model is worse than simply predicting the training mean. MAPE below 10% is generally considered excellent, 10–20% good, and above 50% poor (Lewis, 1982).

Locations with very few observations ($n_{\text{obs}} < 10$) should be interpreted cautiously. ARIMA and RF both require a minimum training length; if this condition is not met, the method is skipped and its error columns are populated with null. The best method is then selected from the remaining candidates.

References

1. Box, G. E. P. & Jenkins, G. M. (1970). Time Series Analysis: Forecasting and Control. Holden-Day, San Francisco.
2. Breiman, L. (2001). Random forests. Machine Learning, 45(1), 5–32. DOI: 10.1023/A:1010933404324
3. Hyndman, R. J. & Athanasopoulos, G. (2018). Forecasting: Principles and Practice, 2nd ed. OTexts, Melbourne. otexts.com/fpp2
4. Makridakis, S., Andersen, A., Carbone, R., et al. (1982). The accuracy of extrapolation (time series) methods: Results of a forecasting competition. Journal of Forecasting, 1(2), 111–153. DOI: 10.1002/for.3980010202
5. Lewis, C. D. (1982). Industrial and Business Forecasting Methods. Butterworth Scientific, London.
6. Bergmeir, C. & Benitez, J. M. (2012). On the use of cross-validation for time series predictor evaluation. Information Sciences, 191, 192–213. DOI: 10.1016/j.ins.2011.12.028
7. Tashman, L. J. (2000). Out-of-sample tests of forecasting accuracy: An analysis and review. International Journal of Forecasting, 16(4), 437–450. DOI: 10.1016/S0169-2070(00)00065-0
8. Armstrong, J. S. (2001). Principles of Forecasting: A Handbook for Researchers and Practitioners. Springer, New York.

Theoretical Background

Time-series forecasting extends the analytical pipeline from description to prediction by producing future values for each spatial location. This algorithm serves as the operational counterpart to Algorithm 7 (Forecast Backtest): where the backtest evaluates which modelling approach best fits the historical data, the forecast algorithm uses that knowledge (or an automatic selection heuristic) to produce actual out-of-sample predictions. Together, they provide a complete model-selection-and-deployment workflow that is rarely available in GIS software outside of specialised statistical environments.

The auto method implements a pragmatic model hierarchy tailored to short and noisy geospatial time series: ARIMA is preferred when the series is sufficiently long ($n \geq 8$), Random Forest when the series is at least $\text{lookback} + 5$ points long, and Exponential Smoothing as a fallback for the shortest series. This hierarchy reflects the principle that more complex models should be used only when the data volume supports reliable estimation. The lookback window for Random Forest is set to $\min(5, \lfloor n/3 \rfloor)$, balancing feature richness against training-data availability.

All three methods generate multi-step-ahead forecasts. Exponential Smoothing produces a constant forecast equal to the last fitted level. ARIMA forecasts recursively through the differencing equation and the estimated ARMA coefficients. Random Forest generates recursive forecasts by sliding a fixed-length window, where each predicted value becomes an input for the next step. This recursive strategy, while propagating errors forward, is the standard approach when exogenous predictors are unavailable.

Mathematical Formulation

ARIMA(p,d,q) forecast. Given the estimated model on the full series, the $h$-step-ahead forecast is obtained by iterating the ARIMA equation forward, setting future errors $\varepsilon_{T+k}$ to their expectation (zero) for $k > 0$:

where $\hat{x}_{T+k} = x_{T+k}$ for $k \leq 0$ (observed values) and $\hat{\varepsilon}_{T+k} = 0$ for $k > 0$, while $\hat{\varepsilon}_{T+k}$ for $k \leq 0$ are the model residuals. Differencing is inverted after ARMA forecasting to return to the original scale.

Random Forest recursive forecast. With a window size $w$ and a trained Random Forest model $f_{\text{RF}}$, the $h$-step-ahead forecast is computed recursively:

Exponential Smoothing forecast. The forecast is constant at the last smoothed level:

Auto method selection hierarchy. For a location with $n$ observations:

where $w = \min(5, \lfloor n/3 \rfloor)$ is the lookback window size for Random Forest.

Forecast columns. The output includes $h$ forecast columns. For step $k \in \{1, \dots, h\}$, the value is $\hat{x}_{T+k}$ computed by the selected method. The method used is recorded in the forecast_method field.

Parameters

Output Description

Interpretation Guide

Forecast values should always be reported with an indication of uncertainty. Since this algorithm produces point forecasts only, users are strongly encouraged to consult the backtest results (Algorithm 7) for the corresponding RMSE at each location as a proxy for forecast uncertainty. A location with a low backtest RMSE and ARIMA as the best method is a strong candidate for reliable forecasting; a location with a high RMSE and Exponential Smoothing as the best method should be interpreted as having inherently unpredictable dynamics.

Multi-step forecasts beyond 5–10 steps should be treated with increasing caution. The recursive forecasting strategy feeds prediction errors forward, so the forecast variance grows with the horizon. As a rule of thumb, forecast up to $h \leq n/3$ steps ahead, where $n$ is the length of the historical series. Beyond this, the forecast is dominated by model drift rather than data signal.

When the auto method is used, inspect the forecast_method field to understand which model was selected at each location. If the backtest (Algorithm 7) identified a different best method for some locations, consider running this algorithm with the explicit method indicated by the backtest results rather than relying on the heuristic hierarchy.

References

1. Box, G. E. P. & Jenkins, G. M. (1970). Time Series Analysis: Forecasting and Control. Holden-Day, San Francisco.
2. Breiman, L. (2001). Random forests. Machine Learning, 45(1), 5–32. DOI: 10.1023/A:1010933404324
3. Hyndman, R. J. & Athanasopoulos, G. (2018). Forecasting: Principles and Practice, 2nd ed. OTexts, Melbourne. otexts.com/fpp2
4. Hyndman, R. J. & Khandakar, Y. (2008). Automatic time series forecasting: The forecast package for R. Journal of Statistical Software, 27(3), 1–22. DOI: 10.18637/jss.v027.i03
5. Brockwell, P. J. & Davis, R. A. (2002). Introduction to Time Series and Forecasting, 2nd ed. Springer, New York.
6. Makridakis, S., Spiliotis, E., & Assimakopoulos, V. (2018). Statistical and machine learning forecasting methods: Concerns and ways forward. PLoS ONE, 13(3), e0194889. DOI: 10.1371/journal.pone.0194889
7. Chatfield, C. (2000). Time-Series Forecasting. Chapman & Hall/CRC, Boca Raton.
8. Shumway, R. H. & Stoffer, D. S. (2017). Time Series Analysis and Its Applications: With R Examples, 4th ed. Springer. DOI: 10.1007/978-3-319-52452-8