Timescape: A Simple Spatiotemporal Interpolation Tool

We developed a novel approach in the field of spatiotemporal modelling, based on the spatialisation of time: the Timescape algorithm. It is especially aimed at sparsely distributed datasets in ecological research, whose spatial and temporal variability is strongly entangled. The algorithm is based on the definition of a spatiotemporal distance that incorporates a causality constraint and that is capable of accommodating the seasonal behaviour of the modelled variable as well. The actual modelling is conducted exploiting any established spatial interpolation technique, substituting the ordinary spatial distance with our Timescape distance, thus sorting, from the same input set of observations, those causally related to each estimated value at a given site and time. The notion of causality is expressed topologically and it has to be tuned for each particular case. The Timescape algorithm originates from the field of stable isotopes spatial modelling (isoscapes), but in principle it can be used to model any real scalar random field distribution.


Introduction
One of the major issues in ecological modelling is the associated occurrence of spatial and temporal variability of the observations underlying the models [1]. Often one has to merge old-fashioned and modern measurements in a coherent mixture, especially when following the evolution of some phenomenon over time. Concerning both the purely spatial and purely temporal dependence of the variables, we have a rich modelling toolbox at hand, ranging from time series analysis [2,3] to spatial statistics [4][5][6]. However, when one cannot decide whether to favour spatial or temporal variability, the toolbox offers fewer options, commanding greater statistical intricacy [7][8][9][10], see for example the moving average processes [11] or the covariance-based methods [12] implemented in R gstat package [13,14]. In some cases we could have a dynamical model at hand [15] so the modelling is generally conducted via partial differential equations integration, as it is often the case with population dynamics [16]. Other trending approaches include neural networks and deep leaning techniques [17,18]. Furthermore, one is often confronted with the problem of model scaling [19], though the computational limitations are less and less important.
We chose to treat our observations as a random scalar field through standard geostatistical techniques, without introducing an explicit dynamic model or any kind of machine learning. Spatial statistics demands that the distribution of our observations obeys more or less strict constraints, depending on the kind of modelling we would like to perform [4,20]. On the other hand, the distribution of observations in ecology if often sparse, both in time and space, especially when merging new data, collected according to a robust spatial pattern, with old observations, typically collected just at spot sites or along transects: data non-stationarity and a spatially skewed distribution of the observations is often encountered [21], so the modelling should be conducted, as a matter of principle, also with sub-optimal datasets, maybe at the expense of accuracy.
In order to deal with entangled spatial and temporal variability, we propose a novel spatiotemporal interpolation technique -the Timescape algorithm -that estimates the spatiotemporal distribution of a given variable (technically, a real random field) from a set of observations. The key idea is borrowed from the machinery of relativistic physics [22,23], which incorporates a notion of causality rooted into its very topological structure, the Minkowski spacetime [24]. Time is treated as new spatial dimension [25], via the multiplication by a constant with the dimensions of a velocity, and we also forbid the instantaneous propagation in space of the field values by filtering the possible sources when estimating the random field.
In a nutshell, the Timescape algorithm works at a topological level, introducing a spatiotemporal distance that mimics the Minkowskian structure, then exploiting any ordinary spatial interpolation tool for actually estimating the field values. Our spatiotemporal distance can be tuned by a couple of parameters, the aforementioned velocity and a causal strictness parameter, plus a form factor related to seasonality. The Timescape distance is exposed in section 2.1 and 2.4 and the actual implementation is discussed at length in [26].
On practical grounds, we sought for a modelling tool that is: 1) not too choosy about the distribution of the input data, 2) capable to accomodate seasonal variability, 3) suitable for both both projected and geographical coordinates, 4) capable to run on affordable hardware in decent times, 5) distributable as open source. The resulting software implementations are a Python module [26] and two Java applications, one for projected [27] and one for geographical coordinates [28]. The software is exposed in section 2.5. A few case studies are discussed in section 3.

Methods
A Timescape model is a discrete representation of a scalar field on a space-time continuum with a temporal dimension and any number of spatial dimensions (figure 1), i.e. a lattice of spatiotemporal cells (here and henceforth the target) whose values are interpolated from a sparse set of observed field values (the source). A Timescape can be regarded as a stack of sheets, each of which is an ordinary discrete spatial model, indexed by the time, as well as a sheaf of time series, identified by their spatial location.
To define a spatio-temporal distance, we mimic the structure of the Minkowskian relativistic metric [22,23,29]. According to the conventions of relativistic physics, we call event any space-time element x = (x o , x 1 , x 2 . . . x n ), where x 0 is the time and x 1 , x 2 . . . x n are the spatial coordinates.
In the following, we will represent the time and only a single spatial component, regardless the actual number of dimensions (generally two, in both geographical and whatever projected coordinates) -see figure 2. The Minkowskian space-time incorporates a notion of causality due to its topology [24]. In order to implement this geometric realization of causality we must impose an equivalent causal structure on our space-time.

Space-time distances
Each event x can be split in its time component x 0 and space components x = (x 1 , x 2 . . . x n ), i.e. x = (x 0 , x). Now define a spatial distance D s (x, y) between the events x and y using any given ordinary metric function d : D × D −→ R + , where D is the spatial domain of interest (a subset of R n or any suitable manifold, e.g. the sphere S 2 of geographical coordinates). A metric d(·) must satisfy the the following conditions: identity: d( x, y ) = 0 ⇐⇒ x ≡ y, non-negativity: d( x, y ) ≥ 0, symmetry: d( x, y ) = d( y, x ) and subadditivity: d( x, y ) ≤ d( x, z ) + d( z, y ) ∀ z ∈ D, also known as triangle inequality. The latter, in particular, gives the metric its distinctive character of nearness measurement and it is often the key property for the convergence of the interpolation algorithms. We define the spatial distance simply as The metrics of practical interest in ecological modelling are the ordinary Euclidean distance in R n for projected coordinates (2) and the geodesic distance in S 2 for geographic coordinates D s (x, y) = R arccos sin x 2 sin y 2 + cos x 2 cos y 2 cos where R is the radius of the Earth, x 1 , y 1 the longitudes and x 2 , y 2 the latitudes of the events x and y [30,31].
Since the temporal coordinates live in (a subset of) R, we could naïvely take |x 0 − y 0 | as the temporal distance between x and y but the dimensions of D s (x, y) (length) and |x 0 − y 0 | (duration) are not coherent. In order to integrate space and time we introduce the first of the Timescape parameters: a time-to-space conversion factor c, i.e. a velocity (∈ R + ) associated with each Timescape model. 1 Thus we define the temporal distance Now we combine the spatial and temporal distances via the Pythagorean theorem, thus assuring the fulfillment of identity, non-negativity, symmetry and subadditivity, taking D 2 t (x, y) + D 2 s (x, y). This function is a well-behaved metric but lacks a causal structure: the point is that the symmetry property commands no distinction between x and y, whatever x 0 and y 0 , so we renounce the symmetry constraint, resorting to what is technically a distance in a pseudo-metric space [32,33]. The distinctive "measure of closeness" character of a metric is given by its subadditivity so, respecting the latter constraint, we impose a directionality to our distance, restricting to the subset x 0 ≤ y 0 , thus introducing the asymmetric function D(x, y) from x to y: In a context of spatial interpolation, what (5) is saying is that something happening at x might be a cause of y but y cannot influence x (i.e. no time traveling in the past). We can be more strict requiring that x 0 < y 0 is not enough for a source in x to propagate up to y, i.e. we forbid instantaneous actions-at-a-distance, introducing the second Timescape parameter, the causal parameter k, positing where k ≥ 0. Said E = T × D the space of all the events, where T ⊆ R contains all the events' times, the ∆(x, y) in (6) coincides with D(x, y) in (5) only in a conical subset whose tip is x ( figure 3). We call this subset C + x , the future causal cone of x. It is the locus of the possible outcomes of something happening at x. Formally, Figure 3. A geometrical interpretation of ∆: ∆(x, y) is the hypotenuse of the triangle formed by D s (x, y) and D t (x, y) if y ∈ C + x . For all the other events / ∈ C + x the distance is infinte. For example: Exploiting the Heaviside theta function 3 the formal definition of the Timescape distance from x to y can be rephrased as 4 note that (7) is parametrized by c and k: each pair of parameters gives a distinct distance, so a different Timescape model. In section 2.4 we discuss in detail how to choose c and k, maximizing the accuracy of the interpolation.

Causal Structure: Topology of the events space
Analogously of C + x , we introduce the past causal cone of x, defined as C − x = y ∈ E ∆(y, x) < ∞ , something occurring at an event in C − x is a potential cause of x. The value of k tunes the causal cones width. The bigger k, the larger the cones, the looser the causal constraint. While the smaller k, the thinner the cones, the stricter the constraint (figure 4). As limiting cases, as k → 0, the cones shrink to vertical lines, while as k → ∞, C + x becomes, E + x the whole future of x and C − x becomes E − x ; k → ∞ corresponds to an infinitely fast propagation of information in E. Figure 4. Sheaves of future and past causal cones. x sits at the common tip of the cones. The cones widen as k increases. For k = 0 the cones shrink to segments and for k → ∞ the cones coincide to the whole past (E − x ) and future(E + x ) of x.
The Timescape algorithm, in a nutshell, consists in filtering topologically the possible causes for each target event from a set of given source events, picking only the source events that fall into the past causal cone of the target event.
Further mathematical details can be found in the manual accompanying the Python Timescape module [26].

Accommodating seasonal variability
In ecological modelling, it is of paramount importance to be able to cope with periodic phenomena. We propose a refinement of the Timescape topological structure that includes seasonality by reshaping the cones by means of a multiplicative positive form factor ψ(t) that shrinks and inflates the cone width as ψ(t) < 1 or ψ(t) > 1, thus modifying the time distance (4) as 3 The Heaviside function is Θ(t) = 1 if t ≥ 0 or 0 otherwise, we also posit finite 0 = ∞. 4 It is important to distinguish the starting event (from x) from the ending one (to y). We cannot simply say "distance between x and y" since ∆(x, y) lacks symmetry; the remnant shadow of symmetry lies in the relation y ∈ C − x iff x ∈ C + y .
Earth 2022, 1 and the Timescape distance (7) accordingly: Figure 5 shows an example of periodic shrinking with period T, using the form factor ψ α , tuned by the parameter α, defined as 5 with 0 ≤ α ≤ 1 (the smaller α, the larger the correction). The resulting "Xmas tree" past cone C − x shows clearly on-and off-season peaks. The trick operated by ψ α (t) in (10) consists in shrinking C − x (k), the regular causal cone of aperture k, down to C − x (αk), i.e. a smaller cone of aperture αk, during the off-season times, thus reducing the number of possible causes of x becomes a union of saucers connected at single events on the cone axis. The n-dimensional shrunk causal cone C − x is obtained rotating the t ψ(t) function's profile about the horizontal axis (figure 5). In bi-dimensional models it is an ordinary solid of revolution; the rotational symmetry is required by the statistical hypothesis of spatial isotropy, while the widening along time is due to the (finite speed) propagation of the modelled field.

The Algorithm
A Timescape model consists in a set of estimated values of a scalar field Φ, generally arranged as cells of a regular lattice (the target T). The field is interpolated from a set of known values (the source S). Each model is characterized by the choice of a spatial distance (the parameters c and k, and the form factor ψ) and the spatial interpolation method. In general, there are no constraints on the distribution of the source events -see figure 6. Some interpolation methods (e.g. Kriging) require the analysis of the variogram and are also able to evaluate the accuracy of the estimate, some others (e.g. Inverse Distance Weighting) are less demanding but do not provide an accuracy figure. The Timescape algorithm operates a topological filtering of source events before the actual interpolation, so in principle one can choose any tool from the spatial statistical catalogue.
The evaluation of a Timescape variogram deserves a note of caution. First of all, each target event x possesses its own input set S x , the causal neighborhood of x (figure 7), defined as so, as a matter of principle, each x ∈ T possesses its own variogram. However, a global variogram can be evaluated by the pairwise comparison of the field values within the source events S. Figure 7. Two target events x and y with their input sets S x and S y .
To evaluate the variogram we introduce the set S δ (h) of the causally connected pairs of S with a distance within the interval of width δ about h: the corresponding variogram value is where S δ (h) is the number of elements of S δ (h). The bin width δ has to be chosen in order to have no overlapping bins, typically δ = max{∆} N , where N is the optimal number of bins [34] and h is taken at regular intervals h n = n− 1 2 δ counting n = 1 . . . N. There is not the usual variogram factor 1 2 in (13) since the asymmetry of ∆ in (12) prevents double-counting. 6 Generally speaking, for each target event x, the value of the field at x is a weighted average of the field values of the source: the actual functional form of the weight W depends on the spatial interpolator of choice, distance being the only common input to all the interpolators. Note that if S x = ∅ there is no way to estimate Φ(x), if this is the case, we attribute a null value to Φ(x). Other than the field value Φ(x), one has to calculate the accuracy of the estimate σ(x); some interpolators provide such a figure, but some others do not, in the latter case, it is advisable to find the accuracy by jackknifing techniques upon the completion of the interpolation, see section 2.4 for details.

Timescape Algorithm
OUTPUT START END Figure 8. Timescape flowchart. Users must choose a target set T, a spatial interpolator, the distance parameters c and k and form factor ψ (thick boxes). The accuracy can be evaluated by the spatial interpolator itself or a posteriori, depending on the chosen spatial interpolator.
Summing up, the Timescape workflow proceeds as follows (a block flowchart is given in figure 8): firstly, one must create an empty target 7 T, choose a distance function (equation 7 or 9) and a spatial interpolator, then: The loop above can be parallelised in a number of ways, since all the events in T are all considered independent of each other. A useful performance figure can be calculated during the interpolation, it is the number η(x) of elements in S x : the larger η(x), the more robust are Φ(x) and σ(x). Upon the completion of the loop one has to estimate the model's accuracy, if it has not already be done by the interpolator itself. If it is not the case, the accuracy can be estimated by jackknifing the source set [35], removing a few random elements at a time, and repeating such process a number of times N [36], a consistent standard deviation estimate is thus [37]: is the model's field estimate at x and Φ k (x) is its estimate on the subset S k . Note that this accuracy estimate is very time-consuming. On practical grounds, (16) can be performed only on a few S k sub-sources, since the running time of each subset is comparable with that of the full Timescape model. Note also that one cannot remove too many events from S since it could lead to too many empty S k sets, i.e. too many null Φ k (x). Further details on the procedure can be found in Python Timescape module manual [26]. The overall accuracy can be assessed by two figures: a global variance σ 2 model and the ratio of null target events over the total target size η model where T * = x ∈ T Φ(x) = null is the set of not-null elements of T. The η model index, in particular, measures the predictive power of the model. σ model is expressed in the same units of Φ while η is adimensional and 0 ≤ η ≤ 1.

Model Tuning
Given the same input set S and output T, the very character of a Timescape model is given by a number of factors: the spatial interpolator, the metric parameters c, k and the form factor ψ. We give some suggestions in order to choose wisely an optimal distance.

Spatial Interpolator
The vast family of Gaussian process regression techniques, (in short, Kriging) provides a wide choice if interpolators, furthermore, such techniques provide an accuracy estimate as a variance figure, so that there is no need to perform the post hoc jackknifing. On the other hand, some old-style deterministic interpolation techniques, like the Inverse Distance Weighting, though lacking the accuracy, are extremely lightweight and fast to perform, so they should be considered, especially when geographic coordinates are involoved, with a spatial distance like (3). 8 As a rule of thumb, it is advisable to pick the same spatial interpolator one would choose if all the events were located at the same time.

Metric Parameters
c and k affect the metric of E: tuning these values is of paramount importance for a successful interpolation. Other than a trial-and error strategy, one should pick the best performing (c, k) pair. What best performing means is a compromise between fast running and accuracy [38]. We recommend a procedure based on residuals minimization by source jackknifing [35]. The idea is based on the cross-estimation of the field values of S with S itself, similar to (16) with S in place of T, for a whole ensemble of (c, k) values, picking the best performing pair (c,k).
Let's define Φ ck (u) the Timescape estimate of Φ(u) -which is known-with parameters (c, k). The estimate is performed on the reduced source set S\{u}. Mimicking (16) we introduce two indicators, the cumulative residuals R 2 (c, k) and the non-null events ratio N(c, k): the best performing pair (c,k) can be found by minimizing R 2 (c, k) and maximizing N(c, k) at the same time. Figure 9 shows an example of R 2 and N surfaces as functions of c and k on a regular lattice. Note that, for a given c, N(c, k) is a decreasing function of k: this is the consequence of the fact that a larger k implies a looser causality (see figure 4). As a rule of thumb, any ensemble of (c, k) values should initially include k 0 = 1, i.e. treating space and time on equal grounds. As of c, if there is a characteristic velocity in the model (transport, diffusion, etc.) one can play around such a value, otherwise a naïve first guess could be, given u, v ∈ S: ensuring a fair coverage of the causal cones. 9

Form Factor
The topology of the event space-time E is controlled by the form factor ψ. The causal cones are actual cones if ψ ≡ 1 or any positive constant, which by the way can be absorbed into k. The form factor is used to modify the shape of the cones and so the topology of E, promoting or depressing the contribution of the events according to their temporal distance. We propose two cases: Threshold: we can impose a threshold T as a maximum time distance requiring c y 0 − z 0 ≤ T, posing ψ(t) = Θ(T − t) whose effect is to close the cone C + y . See figure 10 for a comparison with a regular C + x , with ψ(t) = 1. Periodicity: a periodic ψ(t) ≤ 1 accommodates seasonality into a Timescape model shrinking the causal cone (figure 5), thus filtering out the off-season events. Given a period T, a straightforward form factor is ψ(t) = cos 2 ωt, where ω = π/T. If one wants to be more or less selective some variations on the theme are at hand (figure 11). Higher even powers of the cosine have sharper maxima and flat minima, so the cone is more selective. On the other hand, roots of cosine give smoother maxima and sharp minima, corresponding to a looser causal constraint. The cross section of the C ± cones is given by the product function t ψ(t).
As a rule of thumb, as a first guess take ψ(t) = 1 for non-seasonal field values and ψ(t) = cos 2 ( π T t) for seasonal fields modelling, then α-tempering the form with ψ α in (10), as discussed therein. The larger the area below ψ, the looser the seasonal selectivity.

Timescape Implementations
Two algorithm implementations have been published: a Python and two Java versions. All the software packages are released under the GNU-GPL license. Other than going deep into the mathematical subtleties of sections 2.1 and 2.2, the manuals accompanying the software distributions treat in detail the practical implementation and describe the code.
In this paper, we only sketch the key features of the software. Both implementations are limited to two spatial dimensions and provide some statistical analysis tools to explore the source dataset and the finished models.

Python
The Timescape Python module [26] consists in a single module, exposing all the objects, methods and functions needed to perform source data analysis, parameters tuning and the actual model interpolation, plus a bunch of functions dedicated to data plotting and exporting. All the examples in section 3 have been computed with this software.
The usage is command line based, but the functions can be easily integrated into users' own code. The default spatial interpolator exploits the PyKrige module functions [39]; a finished model is a Python object that can be saved in binary format and can be included in a GIS workflow as a multi-layered geotiff image [40], one band per time sheet. The storage is RAM-based, so this is the limiting factor for the number of target cells.
The Python package is designed to be easy configurable, allowing users to code their interpolators. Predefined interpolation methods include Kriging (examples in sections 3.1 and 3.2) and Inverse Distance Weighting (sections 3.3 and 3.4). Other than the Euclidean (2) and Geodesic (3) spatial distances, square and diamond metrics are provided. The allowed form factors are ψ ≡ 1 and a periodic ψ = cos 2 πt T (eq.10, example 3.2). Defining other distances or form factors require a non-trivial intervention on the code.

Java
The Timescape Java implementation consists in two stand-alone applications, one for projected coordinates, the local version [27,41], and one for geographical coordinates, the global version [28,42]. These are distributed as self-consistent jar executables.
The storage of the Timescape models is based on a Hibernate midlayer [43] that manages the connection with any compatible relational database management system (MySql [44] scripts are included in the distribution). The user interaction is mediated by a Graphical User Interface based on standard swing widgets. The database storage is more flexible, in that each target event is just a record of a table, but there is the need of a dedicated database service.
The Java applications allow a full customization of distance and form factor via a javascript-like scripting language but, since it requires run-time interpretation, the performance is dramatically affected. It is also possible to include ancillary variables along with -or even instead of-the source field values.
The java applications come with a utility for the configuration of the database connection and other parameters. A set of tools is also provided to explore, subset and export the finished models.

Results
We present four case studies. They cover a variety of coordinate systems, time scales and peculiarities encountered in ecological geostatistics: projected and geographical coordinates, seasonal effects, uneven distribution and number of observations. Three examples come from the field of stable isotopes spatial and temporal distribution, also known as isoscapes [45], an increasingly widespread tool in the complex field of spatial ecological analysis [46][47][48]. The measured isotope ratios are usually expressed as differences with respect to a given reference standard, the so-called δ-notation: where R is the molar ratio of the heavy n X stable isotope vs the lighter, most abundant m X, the factor 1000 rescales the otherwise small numbers to easily readable figures; the deltas are adimensional; for example, δ 13 C represents the abundance of 13 C vs 12 C in a given sample. A positive δ 13 C indicates an enrichment in 13 C vs the common 12 C with respect to the standard, while a negative δ 13 C denotes a 13 C depletion. The interpolated field in such isoscapes is the isotopic abundance (δ 2 H, δ 15 N, δ 18 O). In one more case study, that demonstrates seasonality, the interpolated field is temperature.  In the following we briefly describe each case.

Fungi δ 15 N
This dataset comes from a study about host trees -mycorrhiza relationship conducted via geostatistical methods exploiting carbon and nitrogen stable isotopes abundances [51]. The model's scalar field represents the hypogeous fungis' ascomata δ 15 N (Tuber aestivum Vittad. -edible black truffles). The complex tree-fungi relations include the nitrogen transport from fungi to the host tree, in particular, the δ 15 N difference, or fractionation, is a clue of tree-fungi symbiosis [52][53][54], with a significant difference in terms of δ 15 N, resulting in a fungal 15 N enrichment with respect to the host [55,56]. From a spatio-temporal point of view, the collection took place within 64 days of continuous fungi-tree interactions, so it was impossible to treat spatial and temporal variability on different grounds. This study led to the development of the Timescape algorithm.
The source consists of 62 events, repeated measures at the same sites, skew-distributed both in space and time. The target is a stack of 64 time sheets counting 128×128 spatial cells each (about one million cells). Figure 12 illustrates the typical Timescape behaviour, time increasing from bottom to top: as a new observation comes in, there is a sort of perturbation in the field that pushes its value towards that of the following ones. This perturbation will eventually consolidate or fade according to the closest subsequent source events. The round shapes are due to the choice of Euclidean metric (2) as the spatial distance.
When treating cases like this one, it is important to check the accuracy of the interpolation, so it is advisable to choose a statistical interpolator like Kriginig, which is capable to calculate a target event's variance estimate. Such variance could be relatively high but it is noteworthy that a time-unaware geostatistical interpolation would be impossible to perform altogether. As it appears from figure 12, accuracy stabilizes about the value 1 as time goes on; the upper right corner is meaningless since, spatially, it lies out of the convex hull of the source events, it would be a non-informative area with any ordinary geostatistical tool.

Lowest temperatures
This case study shows the ability of the Timescape algorithm to cope with seasonal variability. The source dataset consists in the monthly lowest temperatures recording of Umbria region, in central Italy [57,58].
The temperature records span 20 years with monthly time resolution (1/12 yr), their spatial resolution increases over time as more and more stations have been added in the two decades considered. The interpolation has been conducted with a periodic causal cone (T = 1 yr) using universal Kriging. The great number of source events, more than 2000, added further computational load: in fact, the throughput is as low as 30 target events per second, see table 3 in section 3.5.
The source events are shown in figure 13. Note the temporal gaps in the measurements. Also, new stations started logging data over time, resulting in denser events from year 1990   [61]; in our case, the Timescape algorithm filtered out the off-season source values, thus reproducing the correct periodic behavior. The plots shown in figure 14 have different starting times, corresponding to the earliest meaningful time (i.e. the first non-null Φ(x)) for the the given site. The time resolution of the series is one month. The B series of figure 14 is initially influenced by outlying low values, due to the influence of a distant mountain station, which eventually becomes irrelevant, as closer measuring stations started logging data.

Extra Virgin Olive oil δ 18 O
This dataset comes from a study on the geographic traceability of Italian Extra Virgin Olive Oil through isoscape analysis [62] based on the spatial distribution of Carbon [63][64][65] and Oxygen [66][67][68] isotopic composition of precipitation and source water [69]. The δ 18 O isoscapes, in particular, are related to the elevation and to the distance from the sea and correlate with the year's precipitation [66,67]; the original results in [62] have been obtained by standard spatial statistical procedures. The Timescape interpolation of the same data gives results comparable with [62] for the first year only ( figure 15, left), while the subsequent years (figure 15, centre and right) do not match well with the source events. Note in particular how the extreme values are averaged out in 2010 and 2011: this example shows clearly that the Timescape algorithm is not a universal tool, in this case, each year's data are self-contained and independent of each other, so mixing three years of observations is not correct. The computation per se, however, is straightforward: the high voxels throughput (more than one thousand per second, see table 3) is due to the relatively scarce number of source events and, mostly, to the computational speed of the IDW. In this case the weight function (14) is simply W(u, x) = ∆(u, x) −1 . It is worth noting that the same model interpolation by Kriging gave rise to a significant number of errors, about 35% of the target events, in spite of the causal cones relatively large width (k = 5). This is due to the fact that each year's δ 18 O integrates that particular year's water intake of the olive tree, and it is unrelated to other years' values, resulting in different values with the same spatial coordinates, which is particularly disturbing for Kriging [4] -the high number of target events' computational errors is an alarm bell signalling a possible inconsistency of the source events' values or a badly chosen interpolation tool, as it is in this case.

Precipitation δ 2 H
This example is based on the renowned Global Network of Isotopes in Precipitation (GNIP) dataset [71,72] from the International Atomic Energy Agency [73]. The collection of isotopic abundance data started in the early 1960s and the stations network has grown dramatically over time. For historical and technical reasons, the spatial coverage is uneven ( figure 16), with the earlier, denser coverage in central Europe.
The interpolation algorithm chosen is a smooth version of the Inverse Distance Weighting method, with a weight function where m 2 is a positive inertial term, that prevents the overweighting of the closest source events. In our case m 2 = 2. Moreover, a maximum number of ten neighbors has been imposed. Nonetheless, the computation is CPU-intensive since the geodesic metric (3) involves slow trigonometric calculations. Limiting the causal neighborhood requires S x to Earth 2022, 1 be ordered by distance, then trimmed retaining only the closest events to x; this is based on the heuristic observation that only the closest elements contribute significantly to the weight functions. Both Python and Java Timescape implementations allow this kind of trimming. This example shows how the interpolation results are affected by the addition of new sampling sites. The dates covered are form year 1975 to 1984 (included). This is of course a very crude model, based only on the GNIP δ 2 H observations, further improvements could benefit of a regionalised regression based on climatic ancillary data [67,74]. An acknowledged source of Hydrogen isoscapes is the web-based application IsoMAP [75], that also provides insightful information about precipitation hydrogen and oxygen isoscapes in general [68].

Performance Analysis
The interpolation throughput has been estimated in terms of evaluated target events per second, typical values range from a few tens to more than one thousand, depending on three factors: the number of source events |S|, the spatial interpolator (see the method column in table 2, Kriging is slower than IDW), the distance function (geodesic is much slower than Euclidean) and the causal cone form factor (periodic is slower than straight). The cone shape and width also influences the number of null target events.  Table 3. Models' performances, measured in target events per second. η model (17) is expressed in percent, running times in minutes' seconds".
The lowest throughput is associated with the temperature model (section 3.2). In fact, this model associates a large source set, a periodic form factor and a Kriging spatial interpolator, each factor adding its computational load.
A fair estimate of the python Timescape throughput can be set roughly about a few hundreds events per second with a standard hardware configuration. A conservative estimate of 100 evt/s is a good starting point for a new model.
We give no figures for the Java applications since their time performance depends essentially on the database connection speed.

Discussion
The outcome of a Timescape model is based on three elements: the source set distribution, the spatiotemporal geometry (i.e. the shape of the causal cone and the parameters c and k) and the chosen spatial interplation method. Picking a spatial interpolator is a matter of personal taste and computational performance; anyway it is worth stressing that the Timescape distance is independent of (logically, it precedes the) spatial interpolation itself, so two models based on the same source, same topology, but different interpolators shouldn't be much different than two ordinary spatial models are.
As for the scale of modelling, using geographical coordinates adds further complexity, due to the geodesic distance (3). If it is possible, 11 it is advisable to divide the whole target set into patches which are tractable with projected coordinates, then merge the results, as it could be done with any ordinary spatial interpolator. 12 We can explore further topological diversions modifying the structure of distance; section 4.1 introduces a few suggestions. Altering the definition of the causal cones has a deep impact on causality, so the results can be unexpected even if the mathematical procedure itself is sound. A possible scenario could be inferring the past from the future (properly, retrodiction), exchanging the role of past end future cones C − x and C + x . Accuracy assessment is the most delicate phase of Timescape modelling. On the one hand, statistical interpolators provide an accuracy figure σ(x) along with the estimated field Φ(x) at running time, but it is also possible to apply jackknifing (16) on the finished target for every interpolator. The Python Timescape module [26] offers a set of functions (see figure 9) for an initial guess of model's accuracy.

Timescape Topology Modification
The distance (7) can be modified in a more general way than than just applying a form factor that is only a function of time (9). As a matter of principle, the past causal cone C − x of x can have any shape, provided that it is a subset of the past of x: figure 2. As a limiting case, C − x could be the whole x's past E − x . Any such cone shape would retain the property ∆(x, y) < ∞ =⇒ ∆(y, x) = ∞ but not necessarily y ∈ C − x ⇐⇒ x ∈ C + y .
x y z y 0 Figure 17. Symmetric double causal cones. Note that ∆(x, y) < ∞ and ∆(y, z) < ∞, nonetheless Another modification of (7) consists in including retrodiction; just glue C − x and C + x by the common tip x as a double cone C x . This restores symmetry in the metric function (x ∈ C y ⇐⇒ y ∈ C x ), but at the price of subadditivity, that only holds if, given x, y, z ∈ E, the three events are mutually causally connected, i.e. x ∈ C y , y ∈ C z and z ∈ C x (the converse is guaranteed by symmetry). Otherwise symmetry is violated, as clarified in figure 17. Furthermore, allowing retrodiction requires a factor one-half in front of variogram (13) to prevent double-counting. Giving up causality altogether, the topological filtering can still be performed, with any x-neighbourhood in lieu of C x double cone, but hardly it would be of any use in ecological modelling. It is worth mentioning that taking the whole E as C x for all x, though respecting both symmetry ad subadditivity, just means performing a standard spatial interpolation in one more dimension, with no topological filtering at all: this is not a Timescape model any more.
Other than geometrically, the distance (7) can be modified according to a set of external factors, i.e. exploiting the knowledge of one or more ancillary variables. For example, the form factor-modified time distance D t (8) could be a function of some weather-related variables other than the time D t (x, y) = c x 0 − y 0 ψ x 0 , y 0 Ψ a 1 (x) . . . a n (x), a 1 (y) . . . a n (y) (23) where the a k are the observations of ancillary variables at the x and y events. Note that ψ x 0 , y 0 is more general than ψ |x 0 − y 0 | in (8).
The Java applications [27,28] allow the mentioned modifications of distance function and the inclusion of ancillary variables, via a scripting language. The Python module [26] can only use symmetrical cones in the event's past, allowing at most C − x = E − x as a limiting case, letting k → ∞. Practically, a large enough k suffices, say k = max D s (u, x)/ min D t (u, x) where u ∈ S, x ∈ T belong to the source and target sets, respectively.

Timescape vs Minkowskian Geometry
Since the seminal work of Minkowski [76] the causal structure of relativistic space-time has been encoded in the so-called light cone. Causally connected events lie within each other's light cone, while no signal can reach two events outside the cones. Furthermore, the speed of light c acts as the obvious conversion factor between spatial and temporal coordinates. So the Minkowskian space-time is naturally equipped with a topological structure that incorporates causality straightforwardly.
The Timescape geometry purposely introduces a light cone like structure in the ordinary, non-relativistic, space-time. The role of the speed of light is shared between our c and k parameters, and the causal connection has to be introduced somehow forcefully, but the resulting topology is almost identical.
As far as topology is concerned, the Timescape and Minkowskian spaces act similarly, but the metric behaviour is completely different. Using our symbols, the Timescape D T and Minkowskian D M squared distances are 13 The minus sign of D M in (24) encodes the light cone structure. Furthermore, D M is not even a distance since D 2 M is not positive definite and, even when D 2 M ≥ 0, it does not behave as a distance; on the boundary of the light cone D t = D s =⇒ D M = 0, which describes perfectly the light in the relativistic arena, but it is far from what an ordinary distance is expected to do. On the other hand, the plus of D T in (24) ensures subadditivity, i.e. D T restricted to the causal cone is a proper distance, in fact, it is just a Pythagorean sum of squared sides (figure 3). Summarizing, our Timescape distance mimics the Minkowskian topological structure, not its metric. The trick is operated by the Θ(kD t − D s ) in the denominator of (7) that makes ∆ singular. 14

A universal tool?
The Timescape algorithm can be exploited, in principle, with any set of observed variables varying in space and time. The Timescape algorithm should not be intended as a universal tool, though. In figure 18 we sketch a decision tree to help decide whether it should be chosen or not, depending on the distribution of input observations. In particular, it has not to be chosen whenever there is a clear temporal or spatial pattern in the source data that allows the disentanglement of spatial and temporal variability by subsetting the input set -this is the case, for example, of the Olive oil δ 18 O case study, 13 The sign of D 2 M (x, y) defines, for each event x, three different space-time regions: the so-called time-like one, corresponding to the interior of our causal cones, the space-like region of causally unreachable events, and their common boundary, the locus of light-like events y whose D M (x, y) = 0, reserved for light (or other massless particles). This is not the kind of space one expects to find in ecological modelling: D T (x, y) = 0 means that x an y coincide in space and time, unlike D M (x, y) = 0 that includes all the light-like events' pairs. 14 The topological correspondence is: Minkowski space-like interval → Timescape ∆ = ∞, time-like → ∆ < ∞, light-like → ∆ < ∞ too. There is nothing special to the events laying on the boudary of the causal cones. section 3.3. On the other hand, when the input dataset can be clearly partitioned as a bunch of time series, the algorithm can be applied, but it is unnecessarily complicated and time-consuming; in this case, the result is comparable with ordinary interpolated one-dimensional time series.

Conclusions
The Timescape algorithm is a modelling tool aimed at the interpolation of spatiotemporal field distributions from recorded observations presenting an entangled variability in space and time. The spatial and temporal records of the modelled field are embedded into a spacetime inspired by the causal topology of the Minkowskian space, the well-known arena of relativistic physics. This is achieved applying the distance function (7) for filtering the causally connected events before submitting the actual field estimate to an ordinary spatial interpolator. The choice of the interpolator is up to the user, as well as a couple of parameters that tune the distance to fit the particular input data. Seasonality, often encountered in ecological modelling, can be dealt with modifying the distance function (7) by means of a form factor (9).
The Timescape algorithm has been already implemented as a Python module and as two Java stand-alone GUI-based applications. The software is distributed under the GNU-GPL open license. The software is meant to be used with georeferenced input data and can be included in a GIS workflow, eventually exporting geotiff files.
Our modelling tool can cope with seasonality by playing with the shape of the causal cone, without subsetting the input set season-by-season. Moreover, the algorithm can be applied with a variety of interpolation methods. Ideally, any spatial interpolator would do, the choice being driven by the input data and the user's taste.