Modeling movement probabilities within heterogeneous spatial fields

Recent efforts have focused on modeling the internal structure of space-time prisms to estimate the unequal movement opportunities within. This paper further develops this area of research by formulating a model for field-based time geography that can be used to probabilistically model movement opportunities conditioned on underlying heterogeneous spatial fields. The development of field-based time geography draws heavily on well-established methods for cost-distance analysis, common to most GIS software packages. The field-based time geographic model is compared with two alternative approaches that are commonly employed to estimate probabilistic space-time prisms— (truncated) Brownian bridges and time geographic kernel density estimation. Using simulated scenarios it is demonstrated that only field-based time geography captures underlying heterogeneity in output movement probabilities. Field-based time geography has significant potential in the field of wildlife tracking (an example is provided), where Brownian bridge models are preferred, but fail to adequately capture underlying barriers to movement.


Introduction
In recent years there has been considerable development and application of the core ideas stemming from Hägerstrand's seminal framework of time geography [14].Time geographic analysis has been conducted in a wide variety of problem contexts including, for example, gender and opportunities [27], public service provision [39], optimizing facility opening times [6], and mapping wildlife ranges [31].Through these applications there has been a number of developments to the core framework as originally posed by Hägerstrand.

LONG
Specifically, extensions to the classical space-time prism have been substantial so as to account for the structure of transportation networks [35], adjust for uncertainty in both space and time [23,25], take into consideration kinematic effects [26], and effectively visualize space-time paths in the space-time cube [24].The body of literature using time geographic methods continues to grow substantially owing to the broadly applicable nature of the time-geographic framework [32].
One of the main limitations of time geography, and specifically the use of space-time prisms, is that by definition space-time prisms delineate only the outer boundary of movement opportunity.Thus, new developments have attempted to quantify the unequal movement probabilities within the space-time prism.The most well developed ideas were originally proposed by Winter and Yin [55,56] who framed the problem in the context of random walks (termed probabilistic time geography).Song and Miller [46] extended the ideas from probabilistic time geography to incorporate a more robust statistical frameworknamely combining time geography with popularized Brownian bridge models.This work has been further extended to account for the discrete structure of networks and develop robust probabilistic network time prisms based on random walk models [47].Probabilistic time geography has also been approached through the use of kernel density estimation, both in planar 2D space [10] and for spatial networks [9].Long et al. [34] proposed a model for probabilistic time geography that simultaneously considered object kinematics, moving away from a random walk-based model.However, each of these models model movement as occurring in a homogeneous environment and fail to consider the context within which movement occurs.
Prior to these developments Miller and Bridwell [37] proposed the conceptual framework for a field-based time geography, which represents a more pragmatic approach to modeling the unequal movement possibilities within space-time prisms.Field-based time geography considers the potential movement limitations as a spatial field which serves as the basis for calculating internal movement possibilities.Miller and Bridwell [37] lay down the conceptual and theoretical framework for building field-based space-time prisms, and discuss the potential of this approach both in the context of movement along a spatial network and across a spatial lattice (i.e., a cost surface).In their example a transportation network is used to demonstrate the non-regular shapes resulting from consideration of variable movement speeds in the construction of network space-time prisms, and this approach is commonly used in the generation of isochrones (lines of equal travel time) in transportation research (e.g, [40]).A simple example of a lattice application is provided on synthetic data, but not further developed.More recently, context aware random walks have been proposed incorporating local decisions based on underlying covariates into probabilistic movement models and drawing on ideas from popular movement models in wildlife movement ecology [1].Context aware random walks represent a novel way of simulating movement patterns associated with important variables relating to the local environment (e.g., barriers, preferred features).However, context aware random walks are not explicitly bounded by the properties of the space-time prism.
The objective of this paper is three-fold: 1) to formalize the definition of field-based time geography for applications on heterogeneous spatial lattices, 2) to develop models for estimating probabilities of movement across a heterogeneous spatial lattice, and 3) provide the derivation of an algorithm for field-based time geography, along with an implementation in a free and open source software environment.Synthetic data are used to highlight the practical challenges of implementing field-based time geography and demonstrate the www.josis.orginfluence of different parameter choices.A case study involving the analysis of wildlife tracking data is provided to further demonstrate the approach.Discussion will centre on the development of the model, rather than on practical inferences.

Time geography
With time geography Hägerstrand [14] produced a conceptual model for understanding the constraints and limitations on movement of individuals within the space-time cube; a three dimensional space-time continuum consisting of geographic coordinates x and y and time (t).The main conceptual building block of time geography is the space-time prism, which delineates the potential locations (in both space and time) accessible to an individual given known start and end points (termed anchors).The space-time prism is constructed by intersecting the forward cone and past cone from the first and second anchor point, respectively.Taking a slice of the space-time prism at a single time point facilitates the mapping of what is termed the potential path space, a spatial representation of limits of movement at a particular time.Conceptually, the potential path space represents the intersection of two isochrones which are defined as mapped lines of equal travel time.Projecting the spacetime prism onto the spatial plane allows one to map the potential path area; which for classic space-time prisms is by definition in the shape of an ellipse.The key parameter in constructing space-time prisms is the individual's maximum travelling velocity which impacts the extent of space-time prisms in space and time.A set of formal quantitative definitions for time geographic analysis is provided by [36].

Cost surfaces and least cost paths
In the classic time geographic model, movement possibilities are limited only by a singular upper bound on movement (maximum travelling speed).However, practically movement is influenced by the characteristics of the environment through which the individual moves.The conceptual development of field-based time geography incorporated this idea into time geographic analysis.For field-based time geography on a regular lattice, a requirement is the calculation of a cost (resistance) surface, which provides a quantified measurement of impedance associated with every location in space (typically two-dimensional space).The inverse of impedance is conductance, which can be used to alternatively represent ease of navigation.Cost surfaces can be isotropic where local impedances are unidirectional (e.g., walking through thick forest) or anisotropic where local impedances are direction specific (e.g., walking up or down a hill [52]).Commonly, impedances are representative of the product of various factors (e.g, [41]).In many cases, the cost surface is used to represent a true cost (e.g., in monetary figures or environmental damage), but cost can also be represented as an impedance in terms of travel time, and thus the cost surface represents the time cost, and it is this definition that is of interest in the context of field-based time geography.
Measures of cost distance are commonly quantified using a node-link representation.With raster datasets nodes represent cells in a two dimensional regular lattice and links represents the connection between a cell and its neighbor.Using various data sources a model of travel time cost can be derived, for example based on the attributes of the envi-ronment (e.g., slope [52]).The generation of travel time costs is highly application specific and will be dependant on the type of individual, their mode of transport, and available datasets used to characterize the environment.Typically, the result of this analysis will be a lattice (commonly a raster) where each cell (i.e., pixel) is attributed a travel cost c i of crossing that particular pixel.
From a network with defined travel costs, path analysis can be developed using well established methods.For a given cell i with neighbor j let c ij be travel cost (in units of time) to go from cell i to neighbor cell j.Typically in this type of analysis movement between cells is limited to the 4 or 8 cardinal directions (i.e., rooks or queens case definition of neighbors) but other neighbor definitions are possible, but may require increased computational resources [38].The least cost path can be computed between any two points using a path optimization algorithm (commonly Dijkstra's algorithm is chosen).

Field-based time geography
The construction of field-based time geography follows classic time geography by considering the intersection of space-time cones.Consider any intermediate time point t between two anchors A{x a , y a , t a } and B{x b , y b , t b }, where t a < t < t b .For a location (typically a pixel) two different accumulated costs are defined: T ai , which is the cost (in units of time) from location A to location i based on the network N (and similarly compute T ib ).
To model the probability of travelling from A through location i at time t to B, a model for what is expected is useful.In field-based time geography, the expectation is that the object will follow the trajectory associated with the shortest-time path between A and B for which the time is computed, termed T * ab .Movement probabilities are then estimated from the deviations, measured as time, from the shortest time path.For any location i and time t the deviation from the shortest-time path, termed ∆T i,t , is defined as: where T * ab is time duration associated with the shortest-time path from A to B and δ t = t−ta t b −ta .This formulation of field-based time geography assumes the object will move proportionally along the shortest-time path, similar to how linear interpolation assumes the object moves proportionally along the beeline [29].The model draws on the theoretical idea that movement will typically follow the path of least resistance [15] which is based on the principle of least effort [59].The location i associated with the trajectory of the shortesttime path at time t will have ∆T i,t = 0.As a location deviates further from movement along the shortest-time path it will have increasing ∆T i,t values.With field-based time geography of interest is an estimate of the probability an object was at a location at a given time -P i,t .The P i,t can be used to study the internal movement probabilities within field-based space-time prisms, and thus based on classical definitions from time geography, the P i,t can be considered in the context of accessibility.If location i is accessible at time t (i.e., T ai ≤ t − t a ; and T ib ≤ t b − t) then location i is within the potential path space (PPS) at time t (i ∈ P P S t ).Any locations outside of the PPS are given P i,t = 0. Several types of further analysis allow the P i,t to be analyzed more practically.First, a map of the P i,t for any www.josis.orggiven t can be used to quantify movement potential at a specific time.Both [55] and [46] use incremental maps to demonstrate how P i,t change through time within a space-time prism.Such a mapping is useful to visualize and analyze the potential movement probabilities at a particular time.
A function is required to transform the time deviations (∆T i,t ) from equation ( 1) into probability values.There are, however, many potential mathematical functions that could be used to define P i,t (see Table 3.1 which is developed after [15,50]).The most straightforward way to model movement probabilities in the field-based space-time prism is to estimate the probability the individual visited location i at time t as proportional to the inverse of the time deviation.However, [17] discusses the growing trend to use inversesquared functions, typically in spatial interaction models.Alternatively, negative exponential functions have the firmest theoretical foundation for modeling the decreasing activities as a function of distance, cost or time [16,17,54].In Table 3.1, the scaling parameter c 1,t is used to standardize the P i,t so that P i,t = 1 at any time t.Scaling the P i,t is necessary to account for variations in the size and structure of the P P S t [46,56].
The tuning parameter (c 2 > 0) is a decay parameter controlling the strength of the decay function (from Table 3.1).Lower values (c 2 ≈ 0) are used to model weaker decay and model locations deviating from the shortest-time path with higher probabilities.Higher values (c 2 0) are associated with stronger decay and model locations deviating further from the shortest-time path with much lower probabilities.In nearly all applied scenarios c 2 will be unknown, but can be empirically estimated from the data (e.g., GPS tracking data) using, for example, a leave-one-out numerical estimation procedure (similar to that proposed by [19] for Brownian bridges).
Calculating the cumulative visit probability (P i ) for any location i over the entire time interval between t a and t b can be done by integrating the P i,t over time.

LONG
In practice, the integral in equation ( 3) is not easy to calculate, but can be approximated numerically by taking a set of n k equally spaced times between t a and t b (i.e., t a < t k < t b ) and performing numerical integration using the trapezoid rule.
where A is what is termed the anchor raster which is defined as 1 in the cells associated with the two anchors and 0 elsewhere.For any space-time prism the sum of the P i is equal to the time budget of the prism, that is This definition of P i is powerful because it facilitates easy interpretation of modeled probabilities relative to the overall time budget and can be interpreted as the expected value of time spent at each location or can be considered as relative visit probabilities (as in a Brownian bridge, [19]).The map of the P i for the entire space-time prism represents the probabilistic version of the potential path area -the projection of the space-time prism onto the spatial plane.

Factoring in anchor location uncertainty
Anchor points are subject to location uncertainty, typically due to data collection technology, which adds a further complication to time geographic analysis.Here, anchor uncertainty is incorporated by assuming that location error can be quantified by a bivariate Gaussian distribution.A bivariate Gaussian filter, with variance σ 2 , is applied to the resulting P i,t values.The parameter σ is used to quantify the level of uncertainty associated with he anchor point locations in a similar fashion to what is done with Brownian bridge models [19].
where δ t is defined as in (1).The function in equation ( 5) models σ 2 t as time varying between anchor points.As time moves away from the anchors σ 2 t decreases to a minimum at the mid-point between the two anchors.Modeling location uncertainty dependent on time provides the desirable effect that as the object moves away from anchor points the location uncertainty has a lesser effect on the modeled movement probabilities, because location uncertainty is most important where actual location data exists (i.e., at the anchors).The modeling of location uncertainty in this way in effect smooths the output P i,t values, with the degree of smoothing dependent on σ, the effect of this parameter is further examined in the following demonstrations.

Implementation in R
The implementation of the algorithm described above is reproduced using equations (1-5).The P i,t are computed recursively for an arbitrary number n k of time points between t a and t b .The P i are then computed from equation (3) via numerical integration and the trapezoid approximation.Increasing n k leads to more precise estimates of the integral in equation ( 3) (see Supplementary Material).A function is provided that facilitates the leave-one-out numerical estimation of the c 2 parameter.The computational speed of the algorithm is discussed at the end of the results.The framework outlined above has been implemented in the statistical software R [51] as a set of functions as part of the wildlifeTG (Time geographic analysis of wildlife movement) package available via GitHub [30].The algorithm draws heavily on the previously developed packages for spatial analysis and cost surface analysis in R, most notably the gdistance package [53].

Three example scenarios
Three synthetic example scenarios are used to demonstrate the new field-based time geography model for estimating movement probabilities within the space-time prism.The first scenario represents the case where a corridor separates the two anchor points.The second scenario represents the case where a circular barrier is situated between the anchors A and B. The third scenario represents the more realistic case of movement through a heterogeneous environment.The scenarios are represented on a 100x200 grid, where location A is at position (40,50) and location B is at position (160,50).The time budget for movement between A and B is 200 time units and n k = 200.The R code for generating the three scenarios is provided as supplementary material.
As a spatial comparison, P i surfaces are compared between field-based time geography (c 2 = 1, σ = 2), the Brownian bridge, and time-geographic kernel density estimation via raster differencing to assess differences between the methods in terms of output probabilities.For the Brownian bridge and time geographic KDE models, the resulting probability surfaces were scaled so that P i = 1 (i.e., comparing probability densities).The distance between any two surfaces is calculated as the Bhattacharyya distance (D B ) [3] defined as: where r 1 and r 2 are respective P i surfaces to be compared.The D B = 0 when two surfaces are identical and increases unbounded with increasing distance between the two surfaces.Output P i surfaces can also be examined with respect to the most probable location of movement at any given time (i.e., the location of the maximum value of P i,t -or the peak of the surface -termed l t ).For both Brownian bridge and time geographic kernel density estimation the l t locations are equivalent to the linear interpolation between the anchors or more specifically, equally spaced points that follow the straight line from A to B [29].However, for field-based time geography, the l t locations will follow closely along the shortest-time path from A to B, and will not necessarily be equally spaced.Thus, as LONG a spatial and temporal measure of similarity between the two surfaces the Euclidean distance between l t locations from field-based time geography is compared to the linear path (associated with thel t locations from both Brownian bridges and time geographic kernel density estimation) for each of the three scenarios.

Comparing the different parameter combinations
The selection of different input parameters will inevitably influence model results in fieldbased time geography.The cost surface defined by scenario 3 (Figure 1c) is used to test differences between different parameter combinations for computing the output P i surfaces.First, the functions for calculating the P i (i.e., those in Table 3.1) are compared visually in terms of their output surfaces, and quantitatively using the Bhattacharya distance and distance from the shortest time path.Second, differences between different values of the tuning parameter c 2 = [0.01,0.1, 1] and σ = [0, 2, 4] (in pixel units) are evaluated.In evaluating c 2 and σ the exponential model is presented (as similar results were observed for the other models) and analysis is restricted to visually comparing differences between P i surfaces from different parameter combinations.Finally, the choice of n k and the pixel resolution of the underlying surface are evaluated in order to assess the impact of the spatial and temporal granularity of analyisis.

Computational speed of the algorithm
The computational speed of the algorithm is examined by testing various parameter choices for each of the four key parameters: the number of intermediate time slices (n k ), the size of the raster grid (i.e., number of pixels) used to generate the cost surface, the magnitude of σ, and the choice of time weighting function.The n k are varied from 10 to 10,000, the size of the raster grid is varied from 4,000 to 4,000,000 pixels, the magnitude of σ is varied from 0 to 10 (in pixel units), and all six of the time-weighting functions (from Table 3.1) are evaluated (n = 10 times).In each test, only one parameter is varied and the other three are held constant.

Three example scenarios
Simply visualizing the output probability surfaces (i.e., the maps of the P i ) provides an initial assessment of field-based time geography in comparison to the other approaches (Figure 2).In particular, field-based time geography model estimates a unique P i surface for each scenario, whereas the Brownian bridge and time geographic density estimation models result in the same P i surfaces regardless of the scenario (Figure 2).Time geographic kernel density estimation results in the lowest maximum probability value (and lower variation in output surface values) compared to field-based time geography and Brownian bridge outputs.Whereas the range of values is comparable between field-based time geography and Brownian bridge methods (based on the parameters used here).
The deviations between field-based time geography and the Brownian bridge and time geographic kernel density estimation are clearly demonstrated using image subtraction.The Bhattacharyya distance was greatest with scenario 1 for both the Brownian bridge www.josis.orgThe distance (in pixel units) between the most probable location (l t ) of field-based time geography and a linear interpolation demonstrates some notable differences (Figure 4).First, in scenario 1 (Figure 4a) it would be expected that there is little differences here as both follow the direct route between the two anchors, however, due to resistance during the middle portion of the corridor, the deviations from zero are a result of changes in movement speed.With the second scenario (Figure 4b) it can be seen that the field-based time geography model adequately captures the barrier, and the l t locations are associated with the southern route around the barrier due to the initial conditions of the scenario).The maximum deviations between the field-based time geography model surface and linear interpolation occurs at the midpoint in time between the two anchors for scenario 2. Finally, in the third scenario (Figure 4c) we see that the maximum difference in l t occurs at time t = 154, and the deviation from linear interpolation varies considerably over time.

Comparing the different parameter combinations
The different functions (Table 3.1) for computing P i resulted in different P i surfaces (Figure 5).The exponential and normal models (see Figure 5b-c) showed the steepest distribution and resulted in the largest maximum values.The inverse, root exponential and log-normal models all resulted in a less peaked distribution surface, with a lower maximum value of ≈ 0.3 (Figure 5a,e,f).The surfaces showed varying levels of similarity when evaluated using D B , and pairwise comparison of each of these surfaces showed that D B ranged from 0.013 to 0.463 (Table 4.2).The exponential and normal models were the most similar, while the inverse and normal models were the most dissimilar (for the parameter combinations chosen here; Table 4.2).
None of the methods for computing P i resulted in l t locations that deviated from the least-cost path by more than about 4 pixels (Figure 6).The log-normal model deviated furthest from the least-cost path of all six methods, with maximum deviations of 3 to 4 pixels (Figure 6).The other five methods had more or less similar results with l t locations deviating from the least-cost path by up to 2 pixels.
The effects of the tuning parameter c 2 and locational uncertainty parameter σ were also tested.The c 2 parameter substantially influences the resulting P i surface (see Figure 7).As www.josis.org  the c 2 value is increased the output surface is more peaked and the probability density is confined closer to the shortest time path.Very low values of c 2 spread the distribution more evenly across locations contained within the potential path space.The location uncertainty parameter results in a similar pattern, where higher values of σ result in a smoother output P i surface.
The granularity of the analysis in terms of the choice of the n k parameter and the pixel resolution of the underlying cost surface can significantly influence the resulting P i surface (see Supplementary Material).As would be expected, lower n k values result in a less continuous output P i surface and blocky artefacts are present at small n k values.The output surface appears comparatively smooth at much higher values for n k .Similarly, a coarser pixel size results in a much coarser P i surface while a finer pixel size results in a much smoother P i surface (see Supplementary Material).

Computational speed of the algorithm
The field-based time geography algorithm as implemented is linearly dependent on both the user-defined number of time slices n k (Figure 8a) and on the size (number of pixels) of the raster grid upon which the calculation is conducted (Figure 8b).The locational uncertainty parameter (σ) showed a substantial increase in computation time in changing from σ = 0 to σ > 0, however further, after σ becomes non-zero, no further computational costs was observed.The choice of function used to derive the movement probabilities has negligible effect on the computation time.From these results, in most applications the main considerations (in terms of computational time) will be trade-offs between the size of the raster upon which calculations are made and the choice of n k .www.josis.org Figure 6: Distance (number of pixel units) of the location of maximum probability (l t ) from the shortest time path for each of the six potential functions (see Table 3.1) using the fieldbased time geography model applied to Scenario 3.

Deriving the conductance surface
To demonstrate the method using a real world dataset the movement of a single caribou in northern British Columbia, Canada was analysed.The data captures the location of a single caribou every 4 hr during the summer (July -August) of the year 2000.Two landscape attributes were used to define the resistance of the landscape to caribou movement: i) a digital elevation model (DEM; http://geogratis.gc.ca/) was used to derive slope estimates and ii) a land cover dataset (EOSD; [58]) was used to define different land cover barriers.These datasets were resampled from their original spatial resolutions (17 m and 25 m respectively) to a 100 m spatial resolution and clipped to the area associated with the caribou movement (Figure 9a-b).From these two datasets a conductance surface was generated which describes the speed at which a caribou can cross a given pixel.The slope variable was modified to reflect caribou movement speed using a modified version of Tobler's hiking function [52] adjusted to caribou mobility based on empirical findings [11].The land cover types from the EOSD data were used to scale the movement speeds associated with different slopes following empirical findings associating caribou mobility with land cover attributes [20].The output conductance surface then reflects the ease at which caribou can navigate the landscape taking into consideration slope and land cover (Figure 9c).The R code (with links to the data) used in this analysis is provided in the supplementary material.

Calculating the utilization distribution
The field-based time geography model requires the definition of four key parameters prior to implementation: 1) the mathematical time-decay model (i.e., from Table 3.1), 2) the number of time slices to use in the algorithm (n k ), 3) the location uncertainty parameter (associated with the telemetry data), and 4) the c 2 time-decay tuning parameter.The exponential model was chosen for this analysis based on it's widespread adoption in the geographical analysis literature and its performance in the initial tests described above.It is important to note that it is the choice of c 2 that is the most important.The number of time slices was chosen to be n k = 100 based on the earlier findings and in consideration of computational speed.The location uncertainty parameter associated with the telemetry data was chosen conservatively to be 100 m, which relates directly to the spatial resolution of the conductance surface.Finally, an empirical estimate of the c 2 parameter was derived from the telemetry data using a similar leave-one-out optimization procedure to that described in [19].The procedure generates a log-likelihood estimate for a user-defined range of values and from which a value for c 2 = 0.002 was selected as optimal.The output P i surface was then computed using the parameter set described above.To compare, the Brownian bridge utilization distribution [19] was also computed using the kernelbb function in the 'adehabitat' package [4].The movement variance parameter (σ 2 ) was estimated using the liker function following the methods described in [19].The location uncertainty was set to 100 m.All surfaces were normalized to sum to 1 for comparison and presented as the √ UD to better show the distribution of values.The 99% volume contour of the UD raster was used to delineate a home range polygon for both the field-based time geography and Brownian bridge methods for further comparison.

Case study results
The utilization distribution (P i surface) produced by the field-based time geography is very similar to that of the Brownian bridge (D b = 0.028; Figure 10a-c) and only small deviations were observed between them in this case.The home ranges were similar in shape (based on the parameter combinations shown here) with the field-based time geography 99% volume contour home range being 46% larger.The most striking differences in the resulting UD and home range polygons are specifically associated with the lake areas encountered during a longer movement period between two more stable areas of movement (e.g., the area in the middle of Figure 10d).In particular, the Brownian bridge assumes movement directly across a large lake, while the field-based time geography model is able to capture this barrier directly into the resulting UD and home range estimates.The ability of the field-based time geography model to capture both hard (such as lakes) and soft (such as variations in landcover/slope) barriers to movement is evident in the results of the caribou case study.

Discussion
The new field-based time geography method as presented improves upon the original idea described in Miller and Bridwell [37] in three fundamental ways.First, Miller and Bridwell only provided a small example of how field-based time geography could be applied to two-dimensional spatial fields and focused almost exclusively on transportation networks.Thus, the developments here provide a substantial contribution in terms of the application of field-based time geography to movement scenarios across continuous spatial lattices.Second, Miller and Bridwell focus exclusively on the temporal component of the output results, that is the mapping of variation in the accessibility timings within space-time prisms.Here, focus is placed on emphasising movement probabilities using a suite of functions commonly employed in accessibility modeling (that is those in Table 3.1).The calculation of movement probabilities from the accumulated cost times represents an extension of current models for quantifying movement probabilities within space-time prisms (e.g., [46,55]), and thus contributes to this expanding area of research.Finally, the field-based time geography model as presented, is fully implemented within free and open source software (R; [51]) making it accessible for other researchers.The tools are available as part of the wildlifeTG package [30] which draws heavily on existing packages within R and also allows users to customize how underlying cost surfaces are generated.Field based time geography captures the underlying complexity of the environment explicitly and uses this information in the estimation of movement probabilities between two anchor locations.Both the Brownian bridge and time geographic kernel density estimation assume that a linear path between the two anchors is the most probable, whereas fieldbased time geography assumes the most probable path is associated with the shortest-time LONG path, based on an underlying cost surface.This development is significant; probabilistic time geography was developed with the idea that within space-time prisms movement opportunities were unequal, and thus attempted to model these probabilistically using random walks.However, random walks are poorly suited for most applications, thus fieldbased time geography offers a more realistic alternative for scenarios where an appropriate cost surface can be defined.Moreover, field-based time geography explicitly consider the time limitations on movement in the construction of probabilistic space-time prisms.

www.josis.org
Six different functions have been proposed for estimating the movement probabilities (i.e., Table 3.1) from the accumulated costs (or times) associated with a movement segment.Each of the functions represent a decay function used to convert the deviation in time from the shortest cost path (∆T i,t ) into movement probabilities, taking a similar approach to that employed with time geographic kernel density estimation [10].Based on the results from the three scenarios, the normal model and exponential model were most similar, whilst the normal model and the inverse model were most dissimilar (based on the chosen parameters).This is not surprising as these different models have been universally employed in many classical spatial interaction models [12] and spatial interpolation [5].However, further examination of the tuning parameter c 2 suggests that the choice of c 2 is what controls the shape of the output probability surface (see Figure 7), rather than the choice model function.A similar effect if found in traditional kernel density estimation, which is commonly used to estimate the underlying probability surface from a sample of spatial points, where the kernel shape is far less important than the choice of the kernel bandwidth [45].Estimating the c 2 parameter (from some available trajectory data) is straightforward using a leave-one-out estimation procedure (similar to that for the Brownian bridge model; [19]) and a tool for this is provided as part of the accompanying R package.
The granularity of analysis, both spatial and temporal, is an important consideration when implementing field-based time geography models.The temporal granularity is associated with the choice of the number of intermediate time slices (n k ) upon which to base the calculation.This parameter can significantly influence the output results, specifically, when n k is too small, it can result in unrealistic discontinuities in the output P i surfaces (Supplementary Material).Thus, it is important that n k be sufficiently large to avoid this effect.However, the choice of n k linearly increases computational time (Figure 5a).Thus, the selection of n k should be made by carefully weighing the computational costs whilst avoiding values which are too low.The choice of spatial granularity (i.e., pixel size) will typically be a function of the underlying data used in the analysis (i.e., variables used to describe the underlying resistance to movement).Coarser spatial granularities will be unable to capture fine-scale movement processes, whereas finer granularities may be unnecessarily costly in terms of computational speed (Figure 5b).Issues with spatial granularity are well known in least cost path analysis, specifically that when data are aggregated into coarser spatial resolutions the least cost path may be altered significantly [13,44].Thus, the effect of spatial resolution on output results should be considered carefully in relation to the scale of the movement process under study.
Location uncertainty in anchor fixes has been incorporated into the model in the form of a Gaussian filter, where the shape of the Guassian filter is dependent on parameter σ.Here, σ has been modeled as a time-varying function following the model employed with popular Brownian bridge models [19].Other formulations for location uncertainty could be used, and might take on a different form.Location uncertainty can be considered absent by setting σ = 0.The result of removing location uncertainty is a much more peaked P i www.josis.orgsurface due to the fact that it is certain the object is located at the anchor points.A similar effect is noticed in Brownian bridge models when the location uncertainty parameter is set to ≈ 0. In practice, the location uncertainty parameter will be based on the data and application, but in wildlife applications employing Brownian bridge models it is common to set the location uncertainty parameter to the standard deviation of the positional error associated with data collecting technology (e.g., GPS error; [19]).
Existing movement models for estimating movement probabilities based on random walks -notably Brownian bridges [19,46] -are attractive because they treat the underlying movement process as stochastic.However, it is important to note that the estimates of movement probabilities that are derived from such models are deterministic.That is, given a set of parameters, for any pair of anchors with the same time budget and the same distance between them the output probability surfaces (i.e., P i 's) will be identical.This is not the case with the field-based time geography model, as it is entirely dependent on the unique constraints imposed by the underlying environment -the context within which the movement occurs.Recent efforts have developed alternative models that consider context within a random walk framework [1], and are capable of similarly considering impedances to movement.Similarly, random walks can be readily incorporated into cost-path analysis in order to model so-called randomised shortest paths [42].The additional development of randomised shortest paths into the field-based time geography framework can be readily implemented using existing functionality within the R package gdistance [53], however, randomised shortest paths are similarly limited as with Brownian bridges in that they do not explicitly consider the temporal bounds of movement, and thus the development of truncated randomised shortest paths (similar to the truncated Brownian bridge of [46]) would be valuable.
The most significant potential application of field-based time geography is to the analysis of wildlife tracking data.Use of field-based time geography offers new potential to improve and refine probabilistic measures of space use -commonly referred to as the utilization distribution (UD; [57]) where the UD is defined by combining each of the P i surfaces for a larger telemetry dataset (the Case Study provides an example).Brownian bridges (and to a lesser extent time geographic kernel density estimation) represent the current state-ofthe-art for estimating wildlife UDs.However, both the Brownian bridge and time geographic kernel density estimation methods fail to capture the heterogeneous nature of the environments within which animals move.In applications involving terrestrial animals, field-based time geography might use slope, land cover, and anthropogenic barriers as variables in defining an appropriate cost surface [28].Similarly defining an appropriate cost surface for avian applications might depend on accurate wind speed measurements [43].Importantly, following Long and Nelson [31,33] field-based time geography can be used to highlight areas of commission error (areas that could not have been visited) in existing UD models, for example due to the presence of barriers [2] or the assumption of unrealistic movement speeds [33].Quantifying barriers in individual-level wildlife movement is still a largely under-developed component in present wildlife movement models.
The use of field-based time geography should be attractive to wildlife ecologists due to the prevalence of cost-path analysis in addressing other types of research problems.To date, cost path analysis is often employed to study connectivity across large spatial extents, for example when identifying corridors connecting important habitat patches [28].The use of cost path analysis is implicit within field-based time geography, but is used at a much finer spatial and temporal scale; to model individual movement probabilities within LONG a space-time prism (for example with GPS tracking data).Given the prevalence at which individual-based tracking is now being used in the study of wildlife ecology [22], it will be important to identify ways field-based time geography can be used in wildlife space-use studies.Future work will develop field-based time geography for application specifically to the calculation of wildlife UDs.
Many studies are interested in studying human movement in non-network spaces where field-based time geography could be employed as an improved model of movement opportunities.For example, Kay [21] employed cost path analysis in their study of human route choices at an orienteering event in the United Kingdom, drawing on similar ideas first proposed by Douglas [8].Similarly, there is increasing research on eco-tourism and fieldbased time geography might be used to model the potential movements of recreationalists such as hunters [48] and hikers [49] in a similar fashion to how terrestrial wildlife might be analyzed.Similarly, field-based time geography could be used to improve probabilistic models of missing person locations in wilderness search and rescue activities (e.g., [7]).Movement data associated with marine vessels could also benefit from the field-based time geography approach where barriers associated with land/shallow water, currents, and speed restricted zones could all be incorporated into the generation of a suitable cost surface.Further, the methods as described here could be applied to existing field-based time geography model along spatial networks (as presented by [37]) to derive movement probabilities based on resistance along the network (e.g., associated with traffic).

Conclusion
A new model for field-based time geography is presented where movement probabilities within the space-time prism are estimated from cost-path analysis, as is commonly employed in two-dimensional lattices.When compared with two existing approaches -Brownian bridges and time geographic kernel density estimation -field-based time geography offers the distinct advantage of characterizing underlying factors that may restrict or limit movement opportunities.Specifically, the application of field-based time geography offers significant potential for studying wildlife movement and estimating utilization distributions when wildlife move within heterogeneous environments where barriers are presentwhich is the case in nearly all potential applications.Similarly, there is significant potential to develop field-based time geography as an analysis tool to study human mobility outside of traditional transportation networks for example associated with outdoor recreational activity.The new field-based time geography method is implemented as part of an existing package in the free and open source statistical computing software R. Three key datasets are included as part of the supplementary material found at https://github.com/jedalong/wildlifeTG.
• m3-the telemetry data for the caribou over the summer season, • dem-a digital elevation model clipped to the study area, • eosd-a land cover dataset (see [58]) clipped to the study area.
First we will use the DevTools package to load in the wildlifeTG package from it's repository on GitHub.The three datasets are stored within the package.library(devtools) install_github("jedalong/wildlifeTG") library(wildlifeTG) First, let's take a look at the caribou telemetry data (Figure 11).The telemetry data was recorded using satellite collars, where a position fix was attempted every 4 hours.These data are for the individual "M3" and were recorded during the year 2000.Here we are only going to focus on the summer season which is defined as the months of July and August.In total there are n = 327 fixes in the dataset.The data is stored as an ltraj object from the package adehabitatLT [4].

library(adehabitatLT) data(m3) plot(m3)
Next, we will look at the DEM and land cover data (Figure 12).The DEM and landcover data have both been resampled to a resolution of 100m, which is adequate for this demonstration (and given the large study area seems appropriate).Both datasets have also been clipped to the area surrounding the caribou telemetry data, so as not to include extra area in the analysis and improve the speed of computation.The data are stored as a RasterLayer object from the package raster [18].

A.2 Generating the cost surface
We are going to use Tobler's hiking function [52] combined with empirical findings on caribou mobility [11] to generate the impact of slope (in degrees) on movement by caribou.We need to generate a conductance surface where the value in each pixel reflects how quickly that pixel can be traversed by a caribou.We start by focusing solely on slope, where Tobler's hiking function is adjusted so that a maximum of 4.5 km/h is used (based on the findings by [11]).We need to use the gdistance package [53] here.Next, the conductance surface (named speed) will be modified according to the land cover attributes (i.e., we will reclassify the EOSD land cover categories onto a scale of 0-1 where 0 = barrier, 1=easy movement).These conversions are based (approximately) on the findings of [20].We can now take a look at the conductance surface in order to see how it delineates barriers on the landscape (Figure 13).The values in the cells represent the speed at which a caribou can cross a pixel in m/s (i.e., this is what is meant by a conductance surface).We can clearly see the location of some major lakes in the region which would be barriers (the root-exponential normal log-normal We will use n k =100 as a fairly good default for the number of time-slices.The location uncertainty for this data (satellite tags) is assumed to be about 100m so we will set σ=100.We will use the "exponential" model as the time-decay function.We will use a provided function (named estc2) to compute c 2 using a "leave-one-out" estimation procedure.The estc2 function produces a curve identifying an optimized (maximum likelihood) c 2 parameter value based on the estimation procedure (Figure 15).Here based on previous runs of the function, the maximum value of 0.01 is used (which means we are estimating c 2 between 0 and 0.01).

www.josis.org
n_k <-100 sigma <-100 timefun <-'exp' #Note: this function is slow (˜2 minutes).c_2 <-estc2(m3,speed,timefun=timefun,max=0.01) Here we can see for this example that the optimized c 2 value is 0.0018367.Once we have estimated c 2 from the data we are ready to implement the field-based time geography model to estimate the animal's UD.The fbtgUD function (for the moment) is extremely slow in part because we are woking with a large raster (˜95,000 pixels), so it is best to let this run overnight.The sum of the output UD pixel values is equal to the entire time budget of the animal?stracking period (i.e., t n − t 1 ).cellStats(m3ud,sum)/(60 * 60 * 24) #Convert seconds to days We will take a look at the output UD in a plot, and we can also plot the square root of the UD to get a better look at the values (Figure 16).The output is a raster which will sum to the time budget of the entire trajectory (i.e., the difference in time between the first and last fix).par(mfcol=c(1,2),mar=c (3,3,3,6)) plot(m3ud) plot(m3udˆ0.5) We can then estimate the polygon home range from this UD using, for examle, the 95% volume contour (Figure 17).A function volras is provided to do this estimation, however there are many other approaches that could be used here.

A.4 Final remarks
Field-based time geography is a new-take on estimating animal UD's that explicitly considers barriers through the creation of the underlying cost surface.This approach is flexible in that any number of environmental layers can be incorporated into the creation of cost surface.Here we used slope and land cover as these represent important factors limiting caribou movement, but other applications would have other variables of importance.As currently implemented the algorithm is too slow to be practical on large datasets covering large spatial extents (i.e., due to the large raster datasets involved).Any suggestions for improvements to the speed of the algorithm would be appreciated.

Figure 1 :
Figure 1: Three scenarios used to test field-based time geographic model; a) corridor, b) circular barrier, c) heterogeneous landscape.

Figure 2 :
Figure 2: a) -c) Output probability surfaces (P i ) from the three scenarios using the fieldbased time geography model, d) the Brownian bridge model, and e) time geographic kernel density estimation.

Figure 3 :
Figure 3: Differences in the P i surfaces for each of the three scenarios.Tiles a) to c) show field-based time geography against the Brownian bridge model and tiles d) -f) show fieldbased time geography against the time geographic kernel density estimation model.The Bhattacharyya distance (D B ) is an overall measure of dissimilarity between the two surfaces ranging from 0 to 1 (where D B = 0 implies identical surfaces).

Figure 4 :
Figure 4: Distance of the location of highest probability (l t ) from the field-based time geography model to the linear interpolated path (the bee-line) for each of the three scenarios.Example P i,t surfaces show time-slices of the field-based space-time prism associated with specific time points during the segment.

Figure 5 :
Figure 5: Output probability surfaces (P i ) from the six potential functions (see Table 3.1) using the field-based time geography model applied to Scenario 3.

Figure 7 :
Figure 7: Output probability surfaces (P i ) for different values of the tuning parameter c 2 and locational uncertainty parameter σ using the exponential model applied to Scenario 3.

Figure 8 :
Figure 8: Computational speed of the field-based time geography algorithm (as currently developed) when varying a) the number of time slices parameter (n k ), b) the raster size (spatial resolution, c) the σ value (in number of pixels), and d) the choice of the P i function.

Figure 9 :
Figure 9: Data used in caribou example: a) slope (degrees) derived from a digital elevation model (DEM) of the study area, b) land cover data for the study area, and c) the derived conductance surface.

Figure 10 :
Figure 10: Utilization distributions for a single caribou during the summer (July-August) of the year 2000 using the a) field-based time geography, and b) Brownian bridge methods, along with the 99% volume contour home ranges; c) the difference between the field-based time geography and Brownian bridge methods, and d) the home range polygons overlaid on the conductance surface of the study area with the caribou relocation points.

[A. 1
59] ZIPF, G. Human Behaviour and the Principle of Least Effort.Addison-Wesley Press, Cambridge, MS, 1949.A Example: Caribou telemetry data in Northern BC, Canada Loading and viewing the dataThis example will demonstrate the use of field-based time geography to estimate the utilization distribution (UD) of an animal which is a common spatial measure of movement and typically used to study the animal's home range.

Figure 14 :
Figure 14: Output of conductance surface with telemetry points.

Figure 15 :
Figure 15: Output of conductance surface with telemetry points.

Figure 16 :
Figure 16: Output UD and square root of UD.

Table 1 :
Potential functions used to derive probabilities from ∆T i (t) in field-based time geography.*c 1,t is a time-varying scaling parameter, c 2 is a time-decay tuning parameter and s = ∆T i,t .

Table 2 :
Bhattacharyya distance (D B ) comparing different functions for deriving probabilities from ∆T i (t) in field-based time geography as applied to Scenario 3 (D B = 0 for identical surfaces).