4D Time Density of Trajectories: Discovering Spatiotemporal Patterns in Movement Data

Abstract

Modern positioning and sensor technology enable the acquisition of movement positions and attributes on an unprecedented scale. Therefore, a large amount of trajectory data can be used to analyze various movement phenomena. In cartography, a common way to visualize and explore trajectory data is to use the 3D cube (e.g., space-time cube), where trajectories are presented as a tilted 3D polyline. As larger movement datasets become available, this type of display can easily become confusing and illegible. In addition, movement datasets are often unprecedentedly massive, high-dimensional, and complex (e.g., implicit spatial and temporal relations and interactions), making it challenging to explore and analyze the spatiotemporal movement patterns in space. In this paper, we propose 4D time density as a visualization method for identifying and analyzing spatiotemporal movement patterns in large trajectory datasets. The movement range of the objects is regarded as a 3D geographical space, into which the fourth dimension, 4D time density, is incorporated. The 4D time density is derived by modeling the movement path and velocity separately. We present a time density algorithm, and demonstrate it on the simulated trajectory and a real dataset representing the movement data of aircrafts in the Hong Kong International and the Macau International Airports. Finally, we consider wider applications and further developments of time density.

1. Introduction

Recent ubiquity and widespread use of modern positioning and context-aware devices have enabled the acquisition of movement positions and attributes of almost any type of moving object, and thus, have produced large amounts of trajectories data [1]. These data are usually collected as a series of trajectories; that is, when an object moves in the basic 3D geographic space of our physical world, the movement of each object can be presented as a tilted 3D polyline in space [2]. Analysis and exploration of large movement datasets is of great significance to many fields, such as the exploration of the movement laws of a particle, the analysis of human behavior, and the search for ‘bottlenecks’ in transportation networks [3]. They study the movement of groups or individuals on different spatial scales, different time scales, and with varying degrees of complexity [1,4].

In wildlife ecology, the decreased size and widespread use of animal tracking labels allow animal ecologists to collect large amounts of trajectory data describing animal movements. These data are usually composed of trajectories in space and time [1,5], which are commonly analyzed and visualized using the methods for home range/utilization distribution estimation [6]. Home range is defined as an area within which animals usually confine their normal activities (foraging, mating and taking care of young) [7]. Another concept associated with the home range is the utilization distribution [8], which is a probability surface on the 2D region that represents the possibility of finding animals in a particular area [1,9,10]. The home range is usually defined by the probabilistic contour of a certain value of the utilization distribution surface. It usually uses a 0.95 probability, but the choice is subjective and can vary depending on the study [1]. However, these two concepts often focus on the spatial distribution of the measured positions only in 2D space and ignore the time series of the measurements.

One notion that has persisted in wildlife ecology is the space usage patterns of animals’ movement. Animals prefer to spend more time at particular locations, or visit a given place frequently, or move slowly/quickly in some areas, making their living environment uneven [11,12,13]. Ecologists are particularly interested in exploring the role of time in this heterogeneous behavior [1,14]. However, most of the methods for home range/utilization distribution estimation [15,16,17,18,19] are only associated with two spatial dimensions, and seldom consider the time dimension; this makes it difficult to visually discover the spatiotemporal patterns of movement.

Thus, another potential development we can see is to extend the calculation of home range/utilization distribution into a real 3D geographical space (i.e., using elevation as the third dimension). This is of great significance to animals moving in the air or in water, and animals often changing their vertical distribution relative to external environmental factors [1,20,21]. We expect that the movement range of birds or marine creatures can be viewed as a 3D geographic space [1,22], and time, specifically the function of time (3D analogy of utilization distribution), would now represent the fourth dimension. A 3D space, with three spatial axes, forms a three-dimensional geographic space, and is used to visualize the spatial aggregation of the collected movement trajectories. This is where the data cube is incorporated into the trajectory data visualization.

We propose 4D time density as an effective method for analyzing and identifying spatiotemporal movement patterns in large trajectory datasets. This approach was inspired by the space-time density of trajectories [1,2] and the home range/utilization distribution concept in wildlife ecology [10,17]. However, we changed the algorithm for density calculation, incorporating the fourth dimension-4D time density into 3D geographic space, instead of the 2D geographic space. The 4D time density is derived by dividing the normalized path length by aggregated velocity. In addition, time density and utilization distribution are conceptually slightly different. Time density is a measure of the frequency and intensity of a space use of a species, which is associated with the concept of utilization distribution, since the space use intensity is large for places with a large probability of finding a given species. We subsequently describe the trajectory division and cube cell construction methods to establish the computational range of time density, where the traveled distance and aggregated velocity are calculated. The method yields results of space use intensity that are highly correlated with the true probabilities of occurrence (i.e., utilization distribution), and successfully depicts temporal variations in density of occurrence [22]. In a real application case, we present an application of the 4D time density of trajectories on a real dataset that represents the movement data of the aircrafts at the Hong Kong International and the Macau International Airports. The results show that the proposed approach produces density volumes that successfully capture temporal variations in the density of occurrence, and visually identify the specific spatiotemporal patterns of the movement of aircrafts.

2. Related Works

In computer programming contexts, a data cube is a multidimensional array of values, generally used to describe multidimensional data, such as time series of image data [23], spatiotemporal sensor data [24], statistical data [25,26], geographic data [27], and simulation data [28]. Even though it is called a ‘cube’, it can be 1-, 2-, or 3-dimensional, or even higher dimensional, and the semantic of the dimension is not necessarily spatial or temporal (https://en.wikipedia.org/wiki/Data_cube). Each dimension of the data cube typically represents an attribute of the datasets, and the cells of the data cube represent the facts of interest [29,30]. As a form of organizing data, a data cube can thus be used to describe a large volume of trajectory data. However, on the one hand, as larger movement datasets become available, the display of a data cube quickly becomes unsatisfactory, especially when there are many trajectories [2]. On the other hand, it is difficult to illustrate the trajectory attributes in space and time, as they involve many aspects. It is therefore necessary to aggregate or simplify such data, and visually explore their attributes [2,31].

Kaya et al. [32] proposed an experimental study to verify whether the advantages and limitations of 2D and 3D representations are valid for spatiotemporal data visualization. They conducted an experiment that identified different scenarios of spatiotemporal data analysis, each of which is visualized with both a density map and density cube techniques. The experimental results revealed that 3D representations perform better than 2D representations when analyzing spatiotemporal data.

Demšar and Virrantaus [2] used the space-time density to solve the problem of overlapping and clutter of trajectories when there is a large amount of data displayed in the space-time cube. However, they provided information about traffic and routes but did not display other movement attributes.

Spretke et al. [33] applied color coding to multiple trajectory segments on a 2D map based on differences in attributes, in order to display different categories of segments. This helped to separate migratory birds’ daytime flight from the night flight and stops. However, over-plotting hindered the detection of spatial behavior along each trajectory.

Ware et al. [34] developed the GeoZui4D system to show the multiple attributes of 3D trajectories of whale underwater movement by using color, texture, and glyphs, but there were still deficiencies, such as over-rendering and view occlusion.

Kraak and Huisman [35] used a combination of a time graph, scatter chart, cartographic map, and space-time cube to identify interesting events. The colored trajectory segments were used to represent the selected attributes through a map and a data cube. However, it only considered a single trajectory, despite their collection.

Gao [36] presented a spatiotemporal analysis framework to explore human mobility patterns and internal urban dynamics. Based on large-scale detailed records of mobile phone calls in a city, the combination of space-time density estimation, spatiotemporal visualization and spatiotemporal autocorrelation analysis not only helped to visually represent spatiotemporal data, but also provided a quantitative analysis to determine the spatiotemporal patterns of human mobility. Therefore, the above methods often focused on a particular attribute, and do not provide sufficient visualization support to study the spatial, temporal, and attribute components of trajectories in real 3D geographical space, according to our research.

Static methods for home range and utilization distribution derivation calculate these two concepts from the sampled locations of animal movement. However, they do not take into account the temporal aspect of the data [1]. The concept of home range was first proposed by Burt [37] as follows: the normal area of daily activity for animals foraging, mating, and taking care of young. The definition, however, is vague about the application of “normality”, and was later revised by other scholars. Kernohan et al. [38] proposed a more comprehensive concept of home range, namely, an area where a certain animal appears with a certain probability (e.g., 95%) over a period of time. Three methods are commonly used to calculate the home range: polygonal methods, such as the minimum convex polygon (MCP); grid cell methods, such as the grid cell counting method; and probabilistic methods such as kernel density estimation (KDE) [39]. In 1949, Hayne [40] introduced the concept of the “activity center”, and emphasized that the biological connotation of home range should emphasize the utilization intensity of animals at different parts within the range of activity.

Since the 1980s, the concept of utilization distribution has been widely utilized in the study of home range [11]. Utilization distribution is a probability distribution surface constructed from the discrete activity sites of individuals acquired at different times [41]. Kernel density estimation is most commonly used to generate a probabilistic surface from animal location measurements [1], which usually uses a standard 2D point kernel density estimation to place a probability decay function at each point, and then incorporates them into the surface, thereby generating a probability estimation of the animal being in each position [42]. In most cases, the decay function is a symmetric Gaussian kernel, and other decay functions are sometimes used [1,17], such as linear, bisquare, Epanechnikov, and Brownian. As mentioned above, most of the methods for deriving home range/utilization distribution were generally based on two-dimensional kernel functions in geographic space, and seldom accounted for time. However, Demšar et al. [1] proposed an alternative geographic visualization approach that extended the concept of home range/utilization distribution into three-dimensional space-time (x-y + t), based on the concept of the space-time cube. However, it only regards time as the third dimension upon a 2D space (i.e., two-dimensional plane x-y) to estimate the home range/utilization distribution.

We thus utilize the ideas of calculating the home range/utilization distribution in 2D space (x-y) or 3D space (x-y + t) to derive our 4D time density (x-y-z + t) in 3D geographic space, which seeks to solve the problems of the cluttered data cube by aggregating trajectories via calculation of their time density. This paper is mainly a methodology article that proposes a new method for geographic visualization. As mentioned above, it is partly inspired by the methodology of home range/utilization distribution; we therefore chose to introduce it in the context of wildlife ecology. However, the above issues exist not only in the research field of wildlife ecology, but also in the many studies of moving objects (e.g., humans). Since this article endeavors to develop a common visualization method for most types of movement objects, we try to dilute the concept of home range/utilization distribution in the following sections of this work.

3. Method of Calculating 4D Time Density

This section describes the algorithm for time density. The simulated trajectories are used to demonstrate the generation, stacking, and normalization of time density in Matlab™ 8.5 software. In Voxler™ 4 software, volume rendering (direct volume rendering) and volume slicing are used to present the density volumes.

3.1. Algorithm for Calculating 4D Time Density of Trajectories

Time density is a measure of the frequency and intensity of space use of a species. Imagine a data cube, as shown in Figure 1, where the coordinate axis X, Y, and Z of the space Cartesian coordinate system serve as its three dimensions. The position thus at coordinate (x, y, z) is recorded as a point, and the movement of each object can be presented as a tilted 3D polyline (left view of Figure 1) [43]. If such a polyline is accurately mapped, we can find the corresponding trajectory segment in one of the 3D geographic spaces, i.e., the cube cell $\mathrm{Va}$, where the time density of the trajectory is calculated. The trajectory $P_1-P_2-\dots -P_m-P_{m+1}- \dots -P_n$ is thus divided and assigned to each cube cell. Here, line segments $P_m-P_{m+1}- \dots -P_{m+i}\left(i=1,2,3\dots \right)$ fall into the data cube $\mathrm{Va}$, where $P_m=\left(x_m, y_m, z_m\right)$ is the point coordinate and $x_m, y_m, z_m$ are the geographic coordinates. $P_{\left(\mathrm{Va}\right)1}$ and $P_{\left(\mathrm{Va}\right)2}$ are the intersections of the cube cell and trajectory. It is specified that the trajectory segments $P_m- P_{m+1}- \dots -P_{m+i}$ pass through the centroid (center) of the cube cell $\mathrm{Va}$, and all cube cells do not intersect and are parallel to the axes of the established Cartesian coordinate system.

The time density around each trajectory in 3D space is calculated as a volume (i.e., the 3D cube cells in Figure 1), the density of which indicates the frequency and intensity of the space use of the moving objects in the cube cell $\mathrm{Va}$. This is associated with the concept of home range/utilization distribution. The time density in the cube cell $\mathrm{Va}$ is given by Formula (1).

$${\check{d}}_{xyzt}\left(\mathrm{Va}\right)=\frac{{\check{L}}_{\left(\mathrm{Va}\right)}}{{\check{v}}_{\left(\mathrm{Va}\right)}} \frac{[\mathrm{m}/{\mathrm{m}}^3]}{\left[\mathrm{m}/\mathrm{s}\right]}=\left[\mathrm{s}/{\mathrm{m}}^3\right],$$

(1)

where ${\check{d}}_{xyzt}\left(\mathrm{Va}\right)$ is the time density of the cube cell $\mathrm{Va}$, and ${\check{v}}_{\left(\mathrm{Va}\right)}$ is the aggregated velocity. ${\check{L}}_{\left(\mathrm{Va}\right)}$ is the normalized path length in $[\mathrm{m}/{\mathrm{m}}^3]$ and is computed as Formula (2).

$${\check{L}}_{\left(\mathrm{Va}\right)}= \frac{L_{\left(\mathrm{Va}\right)}}{A_{\left(\mathrm{Va}\right)}},$$

(2)

where $L_{\left(\mathrm{Va}\right)}$ is the total path length within the cube cell $\mathrm{Va}$, and $A_{\left(\mathrm{Va}\right)}$ is the volume of $\mathrm{Va}$ in ${\mathrm{m}}^3$. The final output graphs can be explained as follows—high values indicate the areas that are frequently visited or where the objects spend more time; low values represent the areas that are not frequently visited or where the objects move faster.

There is another method to calculate time density. If the point records are sampled using an equal time interval $∆$t, and no observations are lost, the velocity is inversely proportional to the time $t_{\left(\mathrm{Va}\right)}$, as given by Formula (3).

$${\check{v}}_{\left(\mathrm{Va}\right)}=\frac{L_{\left(\mathrm{Va}\right)}}{t_{\left(\mathrm{Va}\right)}}=\frac{L_{\left(\mathrm{Va}\right)}}{{{\displaystyle \sum }}_{i=0}^m{n}_{i\left(\mathrm{Va}\right)}·∆t},$$

(3)

where $n_{i\left(\mathrm{Va}\right)}$ is the weight of the i-th point. This will produce a simple time-density estimation given by Formula (4); it produces the same result as Formula (1).

$${\check{d}}_{xyzt}\left(\mathrm{Va}\right)=\frac{{\check{L}}_{\left(\mathrm{Va}\right)}·{{\displaystyle \sum }}_{i=0}^m{n}_i\left(\mathrm{Va}\right)·∆t}{L\left(\mathrm{Va}\right)}=\frac{L\left(\mathrm{Va}\right)·{{\displaystyle \sum }}_{i=0}^m{n}_i\left(\mathrm{Va}\right)·∆t}{A\left(\mathrm{Va}\right)·L\left(\mathrm{Va}\right)}=\frac{{{\displaystyle \sum }}_{i=0}^m{n}_i\left(\mathrm{Va}\right)·∆t}{A\left(\mathrm{Va}\right)}\left[\mathrm{s}/{\mathrm{m}}^3\right],$$

(4)

Here, the time density is calculated by two methods, namely Formulas (1) and (2), and Formulas (3) and (4). Two methods are given and compared in the article. The first method uses Newton’s formula, $t=s/v$, to calculate the density indirectly, which derives 4D time density by dividing the normalized path length by the aggregated velocity.

We can clearly see that the time density actually indicates the time consumed by moving objects in the spatial region represented by the unit cube cell. Therefore, it is possible to directly sum the time differences between adjacent sampling points to obtain the total consumed time. The second method (i.e., Formulas (3) and (4)) is exactly such a derivation process. If the point records are sampled using a regular time interval $\Delta t$, the time density can be calculated by summing the weighted time within the cube cell via Formula (4), which considers the weight of time at each point in calculating the total consumed time. There is similar conclusion for irregular sampling intervals. The derivation of Formula (4) is simple and intuitive. However, the advantages of using Formula (1) instead of Formula (4) are as follows:

(1)

Since we use the path length and aggregated velocity to calculate the time density in Formula (1), the density should always be at the same scale, which allows us to compare the analytical output for different sampling intervals. Therefore, it is disadvantageous to create the time density directly through the time difference between adjacent points, because the calculated time density is nonuniform (i.e., not at the same scale) for the different sampling intervals in a trajectory datasets.

(2)

Individual modeling of path length and velocity can provide more possibilities to include various explanatory variables to explain the behavior of the moving objects. In kinematics, the velocity can be the function of various environmental conditions (e.g., elevation, terrain, season, wind, obstacles, road conditions, etc.) [44,45,46,47]. If the velocity and actual movement trajectory are modelled more accurately by using additional movement attributes, we should be able to generate a finer time density map for a given spatial region.

However, as this is the first expansion of the trajectory data considering the time dimension in 3D space, these advantages will not be discussed in detail herein. In the following sections, we present the detailed algorithm for calculating the time density. The pseudocode of the algorithm is shown in the following text.

Algorithm: 4D time density algorithm

arrange point records by time sequence and convert them to curves;   calculate cube cell size $D_{\left(\mathrm{Va}\right)}$ depend on the sampling resolution and spatial extent;   TotalTimeDensity = 0; AverageTimeDensity = 0;   for each trajectory    TrajectoryTimeDensity = 0;    allocate cube cells around the trajectory;     for each cube cell      calculate NormalizedLength of the trajectory curve;      calculate AggregatedVelocity by aggregating the velocities;      TrajectoryTimeDensity = normalized (NormalizedLength/AggregatedVelocity);     end    TotalTimeDensity = TotalTimeDensity + TrajectoryTimeDensity;   end   TotalTimeDensity = normalized TotalTimeDensity;   AverageTimeDensity = normalized (TotalTimeDensity/Number of trajectories);

3.1.1. Curve Fitting of the Trajectory

The trajectory of a moving object in 3D space is usually a continuous curve, but the actual sampling process produces a series of discrete points $P_1, P_2, \dots , P_m, P_{m+1}, \dots ,P_n$ [2,48]. The original continuous trajectory is thus usually approximated by a series of straight lines between the sampling points, which is not a real movement path. For our algorithm, we read the point records through the time sequence to generate a trajectory polyline $P_1-P_2-\dots -P_m-P_{m+1}-\dots -P_n$, and then use the cubic spline function in the Matlab™ software to fit the polyline into a continuous curve $S_{\left(\mathrm{Va}\right)}$. The algorithm reads each trajectory sequentially, i.e., one after the other. We use the cubic spline method instead of inserting multi-lines directly, based on the following considerations:

(1)

The fitted path passes through all the sampling points, and is maximally in accordance with the actual trajectory path.

(2)

As mentioned above, 4D time density is derived by dividing the NormalizedLength by the AggregatedVelocity. If the velocity and actual movement trajectory are modelled more accurately, we should be able to generate a finer time density map for a given spatial region.

3.1.2. Estimation of the Size of Cube Cell

The time density is derived in each cube cell. Consequently, we need to choose a suitable size $D_{\left(\mathrm{Va}\right)}$ of the cube cell to provide maximal information content while meeting the goals of spatial analysis. On the one hand, it is clear that the higher the resolution, the higher the accuracy of the rendered path. On the other hand, if the size of the cube cell is too small, there are not enough observations, and the output of the analysis will become noisy since it only emphasizes a part of the point, which is difficult to interpret visually [28]. Recall that we need more points in a cube cell to calculate the path length and velocity of the points accurately. To meet these two requirements, we propose the following geometric principle: the algorithm reads all trajectories sequentially, calculates the distance between adjacent points on the trajectories, and selects the third quartile as the cube cell size. This means that for approximately 75% of the points in the datasets, there will be at least two points within a cube cell. The remaining 25% are interpolated by dichotomy, until there are at least two points in the cube cell. According to this heuristic, we estimate the appropriate cube cell size $D_{\left(\mathrm{Va}\right)}$.

3.1.3. Calculation of the Normalized Path Length

Our algorithm allocates the cube cells around the trajectory from the first point of a trajectory. Once the size of the cube cell is estimated, we can derive the length of the trajectory curve within each cube cell $\mathrm{Va}$. It is first necessary to calculate the intersection coordinate $P_{\left(\mathrm{Va}\right)}$ of the trajectory $S_{\left(\mathrm{Va}\right)}$ and the cube cell $\mathrm{Va}$. The calculation of intersection uses the function on each segment, since the cubic spline function is a piecewise function. We calculate the length $L_{\left(\mathrm{Va}\right)}$ of the trajectory curve $S_{\left(\mathrm{Va}\right)}$, which yields the distance traveled by the object in the cube cell $\mathrm{Va}$. The NormalizedLength is then calculated by using Formula (2).

3.1.4. Calculation of the Aggregated Velocity

The next step is to calculate the AggregatedVelocity in the cube cell Va. Therefore, the velocities at the intersections $P_{\left(\mathrm{Va}\right)1}$ and $P_{\left(\mathrm{Va}\right)2}$ need to be interpolated by the adjacent points. Since this is the first realization of time density, we use the simplest interpolation function, linear interpolation, to insert velocities at the intersections. The AggregatedVelocity for each cell is subsequently obtained by averaging the velocities at all points within the cube cell. However, we can use other, more complex functions, such as kriging interpolation. In geostatistics, this interpolates predictions at unsampled locations according to the spatial correlation between adjacent geographic objects.

3.1.5. Derivation of the Time Density

In the final step, the estimation of the time density is obtained by dividing the NormalizedLength ${\check{L}}_{\left(\mathrm{Va}\right)}$ by the AggregatedVelocity ${\check{v}}_{\left(\mathrm{Va}\right)}$. The algorithm calculates the density of each trajectory separately, and in the meantime sums their corresponding densities to TotalTimeDensity—the time density for the entire set of trajectories. This process is implemented using 3D map algebra (i.e., voxel-by-voxel). A new trajectory is added once the moving object accesses the cube cell. Therefore, the TotalTimeDensity is the sum of the densities of multiple trajectories. When the number of trajectories is large, the corresponding TotalTimeDensity is large, due to the stacking effect of time density. The final step in calculating both AverageTimeDensity and TotalTimeDensity is to normalize the density to the range of [0, 1]. Section 3.2 demonstrates the generation of AverageTimeDensity and TotalTimeDensity in detail.

We propose a technique to estimate time density by separately modeling movement path and velocity in an established cube cell. In summary, time density can be derived from the following six steps:

• Generate trajectories from point records (arrange point records by time sequence, and convert them to curves);
• Calculate the appropriate cube cell size based on the sampling density and the spatial extent;
• Derive the total length of the trajectory in the cube cell, and normalize it;
• Interpolate the velocities at intersections by using linear interpolation, and aggregate the velocities to obtain the average velocity of each cube cell;
• Derive the time density by dividing the normalized path length by the aggregated velocity using Formula (1);
• Generate the total time density volume by using 3D map algebra, and normalize it.

3.2. 4D Time Density of Simulated Trajectories

Figure 2 shows the TrajectoryTimeDensity (i.e., 3D cube cells as described above) of four simulated trajectories. The grayscale value of each cube cell represents the time density of the trajectory. The higher the grayscale value, the higher the corresponding density. The density of each trajectory is gradually increasing along the positive direction of the Z axis, from 0.1 to 1. The internal trajectory curves are visualized by setting the transparency of the cube cells surface.

When the algorithm handles all trajectories, their densities are added to TotalTimeDensity. This process is implemented using 3D map algebra (i.e., voxel-by-voxel). Figure 3 shows the time density for four simulated trajectories. In Figure 3a, the TotalTimeDensity volume is shown using the same 3D cube cell as in Figure 2. The grayscale value of some cubes becomes larger, while others decrease. In Figure 3b, the trajectories are overlaid on the TotalTimeDensity volume. In these two views, there is a time density hot spot in the central region of the simulated area due to the stacking effect of density. The TotalTimeDensity is large in this region, which represents the trajectory convergence phenomenon and implies that many objects have visited the area. However, it should be noted that visitors have only accessed the area, but the time of their visit is not known.

In Figure 3a,b, the time density of the stacking region is the sum of the densities of multiple trajectories that entered the region. It can also be explained as the product of the AverageTimeDensity and the number of trajectories (i.e., the number of visits). The number of trajectories is calculated from the GPS positioning data of the trajectory dataset, which is the number of times that the object enters a certain spatial region (i.e., cube cell) during the statistical period. The AverageTimeDensity represents the overall average of the time that visitors spend in the region. The relationship between them can be expressed as Formula (5).

$${\check{d}}_{xyzt}\left(\mathrm{Va}\right)={{\displaystyle \sum }}_{i=0}^n{\check{d}}_{xyzti}=n · {\overline{d}}_{xyzt},$$

(5)

where ${\check{d}}_{xyzti}$ is the density of the i-th trajectory, $n$ is the number of trajectories in the density stacking region, and ${\overline{d}}_{xyzt}$ represents the AverageTimeDensity. A new trajectory is added once the object revisited the cube cell every time; therefore, the number of visits and number of trajectories mentioned later are the same concept. In Figure 3c, the TotalTimeDensity of the stacking region is decreased after it is normalized by the number of trajectories, and is even smaller than the density of the nonstacked trajectories in the upper half region of the simulated area. Therefore, in the two factors that affect the TotalTimeDensity, the number of visits occupies a larger weight, but the average time occupies a smaller weight; that is, the number of visits is large, but the average time spent by visitors is relatively short, resulting in the obvious contrast between Figure 3b,c. This is an interesting movement pattern which we need to explore. In Section 4, we present the visual and analytic capabilities of TotalTimeDensity and AverageTimeDensity through a real application case.

3.3. Visualizing the Time Density Volume

When studying the distribution of temperature, pressure, or humidity of a certain space in the physical system, it necessary to describe them only by the algebraic quantity in mathematics. The field determined by these algebraic quantities (i.e., scalar functions) is called a scalar field (https://en.wikipedia.org/wiki/Scalar_field). The time density of the trajectories is a volume in which the value of each point can be uniquely described by a scalar calculated from NormalizedLength and AggregatedVelocity. For one trajectory, the density of a certain cube cell is unique, but for the stacked trajectories, the density is the sum of the densities of multiple trajectories. Therefore, many density volumes will produce visual relationships such as inclusion, occlusion, and masking, making it difficult to visualize on a two-dimensional display [2,48]. Volume visualization is an important research field of information visualization; it is common in the medical domain (e.g., a 3D CT lesion distribution images and 3D ultrasound vision) [49] and in the fields of geology [50], mining (e.g., 3D map of the location, direction, length, and depth of an underground mine) [51,52], and meteorology (e.g., 3D representation of the atmosphere, such as cold air masses, hurricanes, and haze) [53]. The three main techniques of volume visualization are direct volume rendering, volume slicing, and isosurfaces [2]. In the following, we use the first two methods to visualize the time density volume in Voxler™ 4 software, which is a professional 3D data visualization system that provides a new way of visualizing volumetric data.

Direct volume rendering is a set of techniques for displaying 2D projections of 3D discrete sampled datasets (usually 3D scalar fields). The technique assumes that the volume is semitransparent, and maps the discrete 3D data field to the corresponding 2D image by controlling the color, light, and observation direction of the volume. We can visualize volume data from different angles by selecting different lighting models and viewing directions [54]. The classical algorithms for direct volume rendering include splatting, maximum intensity projection, shear-warp, volume-ray casting, etc., among which volume-ray projection is a classic algorithm. It emits projected light from each point of the projection plane, passes through the 3D data field, and subsequently calculates the intensity of the attenuated light by the light equation, and then draws it into an image [55].

Volume slicing or visualization with clipping planes (not necessarily a plane) is a surface that is colored, based on the value of the volume data in the region of the slice position [54]. The clipping plane is useful for detecting the interior of the volume data, which removes those unimportant regions to expose the internal patterns [56]. Volume slices are widely used in the medical industry. For example, the physician can easily display the 3D distribution of tumor lesions or of the brain’s blood stasis by using the 3D CT image [57].

The TotalTimeDensity volume and the AverageTimeDensity volume are visualized using these two methods in Figure 4. In Figure 4a–c, the time density volumes are displayed by direct volume rendering using the grayscale color and rainbow color schemes. Figure 4a shows the intersection of cube cells of the four trajectories of Figure 3. The non-intersecting parts are set to be invisible here. Obviously, the intersection parts are symmetrical about the center point. In Figure 4d,e, the density volumes are sliced in different directions, some of which are generated using two slicing planes in two different directions. These different volume visualizations can help visually identify the space usage patterns.

As mentioned above, the TotalTimeDensity is determined by the AverageTimeDensity and the number of visits. In different spatial regions, there may be significant differences between these two factors, which in turn reflect the differences (e.g., social, economic, cultural or ecological aspects) of different spatial regions. Examples include a central business district (CBD) and a residential area, an animal’s habitat and breeding ground, as well as a troposphere and stratosphere. Based on these two factors, we can therefore explore the behavioral characteristics and space usage patterns of moving objects in different scenarios. Potential applications across a wide range of fields include marketing promotion, advertisement serving, urban planning and construction, ecological and species conservation, etc. [58].

4. Application: 4D Time Density of the Aircraft Movement Trajectories

This section presents an application of time density on a real dataset that represents the movement trajectories of aircrafts at the Hong Kong International and Macau International Airports. We will subsequently explore the space use patterns and characteristics of the aircrafts of the two airports for the Pearl River Delta airspace.

4.1. Study Area

The Hong Kong International and Macau International Airports are located in the Pearl River Delta Metropolitan Region (PRD) in the southeast of Guangdong Province, China. PRD is the low-lying area surrounding the Pearl River estuary, where the Pearl River flows into the South China Sea. It has convenient land, sea, and air transportation, and is known as the “South Gate” of China [59].

The specific location of the Hong Kong International Airport (ICAO: VHHH; IATA: HKG) * is Chek Lap Kok, Lantau Island, New Territories, Hong Kong, with geographical coordinates: 22°18′32″ N and 113°54′52″ E. There are two runways in the airport with azimuths of 70 degrees and 250 degrees, respectively, both of which are 3800 m in length and 9 m in elevation. Hong Kong is the third largest financial center in the world, and an important international financial, trade, and shipping center. In 2016, VHHH handled 70.5 million passengers and 4.52 million tons of air cargo. At present, more than 100 airlines provide more than 1100 flights per day from the airport, to over 190 cities across the globe [60]. The geographical coordinates of the Macau International Airport (ICAO: VMMC; IATA: MFM) are 22°08′58f″ N and 113°35′30″ E, which is approximately 40 kilometers from the Hong Kong International Airport. There are two runways in the airport: Runway 16/34 is 3360 m long, and the taxiing runway is 1460 m long, both of which are 6 meters in elevation. Its designated passenger traffic is 6 million per year [61].

* ICAO: International Civil Aviation Organization; IATA: International Air Transport Association. VHHH and HKG are the ICAO airport code and IATA airport code for Hong Kong International Airport respectively.

4.2. Data Description

We collected trajectory data for 30 days and 1235 take off and landing flights at VHHH and VMMC. The research data was purchased and downloaded from the flight real-time tracking website. The data service agreement allows researchers to use these data for personal research, as opposed to commercial, purposes. The raw dataset contains movement attributes, such as the call sign, latitude, longitude, elevation, UTC, speed, and azimuth of the flight. The time resolution is approximately 20 s. However, the position (latitude, longitude, and elevation) and velocity attributes only are used for our experiments. According to the source/destination and departure/arrival conditions of the flights, the datasets is divided into different subsets to calculate the time density of different subsets of flight data. Here, we have removed the sampling points generated when aircrafts move on the apron, in order to avoid affecting the time density calculation.

Figure 5 shows the flight trajectory points in two different ways. Figure 5a is a 2D map of 1 day of flight data (12 December 2017) for the two airports. This display ignores the elevation information of the sampling points, and the 3D spatial patterns cannot be identified. We place these trajectory points in a three-dimensional data cube (Figure 5b), where we can clearly see the 3D spatial position of the points. There is also a map of the study area attached to the x-y plane of the cube, as a geographical background of the movement trajectories. The location of the two airports is also marked on the map, where position A represents VHHH and position B represents VMMC. The geographical coverage of the flight trajectory is 177 km $\times$ 186 km in the E-W and the N-S directions. Since the elevation difference of the flights’ movement range is too small compared to its length and width, the elevation difference is properly amplified to better observe the changes of the movement patterns in the vertical direction. However, this does not affect the overall view. Considering the spatial sampling resolution and the movement range of the flights, the appropriate cube cell size is estimated to be 800 m (see the principle in Section 3).

4.3. Results: Identification of Movement Patterns from Density Volumes

We apply direct volume rendering and volume slicing techniques to generate time density volumes for the flight trajectory. As expected, we are able to visually identify the specific movement patterns of aircraft movement in the density volume.

Figure 6 shows the visualization of different TotalTimeDensity volumes. Figure 6a,b are the TotalTimeDensity volumes of arrival and departure flights, shown with direct volume rendering, respectively. The display in Figure 5b is cluttered with a high amount of overprinting, especially in the airspace near the airport. Both these displays are clearer than the traditional data cube (Figure 5), and the spatial-temporal patterns are more visually recognizable. The runways of VHHH are approximately in the east-west direction; therefore, the flights followed a trajectory from the south to the airport, then turned right at position a (Figure 6a), and finally entered the airport from the west and landed on the runway. However, the departing flights leave the airport from the east. The runways of VMMC are approximately in the north-south direction. All flights enter the airport from the south and take off from the north (position b in Figure 6b). These movement patterns and characteristics are clearly shown in Figure 6a,b.

The hot spots of time density are even more eye-catching when the TotalTimeDensity volume is displayed using a horizontal clipping plane through the average elevation of the VHHH and VMMC (Figure 6c). It can be clearly seen that VHHH is busier than VMMC, because VHHH has a large number of daily flights, resulting in the high-density areas in red near the VHHH. Compared with VHHH, VMMC has fewer flights per day, so it is shown in Figure 6c as areas with low-density that are light green to blue. These patterns can also be observed in Figure 6d, which visualizes the TotalTimeDensity volume with multiple clipping planes in different orientations (mainly vertical clipping planes).

Figure 6e,f shows the density volume of the number of visits and AverageTimeDensity. The number of visits is the number of times that the flights enter the region $\mathrm{Va}$ within 30 days. In the airspace near the airport, all flights generally land and take off along the fixed air routes. In the meantime, flights have low speed and encounter heavy traffic, resulting in the accumulation of a large number of trajectories (the number of visits) and a high AverageTimeDensity. As the distance of the aircrafts from the airport increases, the air routes of the flight gradually become radial, and in the meantime, the flights move at a higher speed, resulting in spatial areas with low-density, i.e., low values of the number of visits and AverageTimeDensity values occur.

Figure 7 shows the TotalTimeDensity for all flights displayed from four different angles. All four panels show the same time density volume, where the top-left image sets the map transparency to 35% and looks up at the cube. In Figure 7, the time density by direct volume rendering is shown in red near the runway of the airport, on the approaching routes before landing, and on the departing routes just after takeoff, indicating that there are spatial regions with high time densities. As the distance from the airport increases, the TotalTimeDensity gradually decreases, and the corresponding volume rendering colors gradually change from orange to green to blue. These phenomena can be clarified by Figure 6e,f. According to Formula (5), TotalTimeDensity is regarded as the product of the AverageTimeDensity and the number of visits. In the airspace near the airport, the number of flights is large, but the flight speed is low; the corresponding TotalTimeDensity, therefore, is large. Similarly, in the airspace far from the airport, these two phenomena are exactly opposite, resulting in a low TotalTimeDensity value. There are two circles (c and d) in Figure 7a, where c is the trajectory of the planes that make an 8-shaped hovering maneuver, and d is the trajectory of the aircrafts that hover a circle. They are hovering arrival flights that are waiting for other planes to land or takeoff. These two movement features can also be confirmed in Figure 6a.

5. Discussion

We use the dynamic properties of the spatial trajectories to derive the 4D time density, and visually analyze and explore the movement patterns of the moving objects in 3D geographical space. Achieving this goal is becoming increasingly easier, as the new sensor and data collection infrastructure supports better access to the context properties (speed, acceleration, direction, curvature, barometric pressure, humidity, altitude, etc.) of movement [62].

There are several issues worthy of further consideration, which can be divided into three categories:

• Complexity of the algorithm. At present, our algorithm uses the calculated cube cell size to partition each trajectory. We specify that all trajectories pass through the cubes’ centroid (center), and the cube cells are parallel to the established Cartesian coordinate system (Figure 1). The time density is then derived by dividing the NormalizedLength by the AggregatedVelocity. Therefore, the generation of cube cells and the calculation of density are somehow complicated. Note that the time density calculation of vehicle movement in the real application case is computationally expensive, i.e., calculation time may take several hours on a standard PC, depending on the spatial extent of the area and the predefined data sampling resolution. In this case, it is not possible to use this method for real-time visualization. We plan to implement different types of density calculation methods to gradually develop this algorithm. For example, we can refer to Demšar and Virrantaus [2] to use a cylinder around the 3D trajectory as a calculation range for time density, while the radius of the cylinder can be adjusted. This is a less computationally demanding solution.
• Visualization techniques of volume data. Time density is a volume, a data type that is difficult to visualize and explore. A typical 3D dataset is a set of 2D slices obtained by MRI, CT, or 3D ultrasound scanners, which can be 1GB in size. Generally, these data are obtained in a regular pattern (e.g., one slice per millimeter), and typically have a regular number of image pixels in the regular pattern [63,64]. We are considering other volume visualization techniques, rather than the previously proposed methods. One alternative is hardware-accelerated volume rendering. Starting with programmable pixel shaders that originated in approximately 2000, people were increasingly aware of the power of multipoint parallel computing, and began to perform general-purpose calculations on graphics chips [65]. With such capabilities, any algorithm with parallel execution steps (e.g., volume ray projection or tomographic reconstruction) can be implemented with great acceleration [66]. Based on the current programmable graphics hardware technologies, the algorithm can use a programmable graphics processing unit (GPU) for hardware acceleration. Another alternative is the optimization technique of volume rendering, including early ray termination, volume segmentation, empty space skipping, etc. [65]. The main objective of the optimization is to skip as many voxels as possible and handle less information. The premise of using these methods is to meet the resolution requirement of the application without affecting the overall view. These methods produce volume visualizations in which patterns can be more pronounced and more quickly manifested, and that can therefore be used to explore our density volume.
• Applications and user-related evaluation. Last but not least, the time density algorithm and all these visualization techniques should be developed with intended users, and should be evaluated for strengths and limitations in a variety of situations. For example, visual validation aimed at different users, experts, and application scenarios (e.g., marketing, urban planning, rescue services, ecological protection, etc.) can be explored. Further research can focus on the empirical analysis of the resulting time density, in which the users and experts can approve or adjust these theories and technologies.

6. Conclusions

We present a 4D time density algorithm, apply the volume visualization techniques to visualize the resulting density volume, and explore the movement patterns of moving objects in space and time through a real application case. The method is based on the concept of a data cube. It regards the movement range of the objects as a 3D geographical space, and displays the resulting density volume in it. The 4D time density is derived by dividing the NormalizedLength by the AggregatedVelocity.

This is a new method for geographical visualization inspired by the aggregated space-time density of trajectories in the space-time cube [2] and the home range/utilization distribution concept in wildlife ecology; however we changed the density algorithm. The ideas of calculating home range/utilization distribution in 2D space (x-y) or 3D space (x-y + t) are utilized to derive 4D time density (x-y-z + t) in real 3D geographic space. As mentioned above, time density has the following geographical indications: high values indicate the areas that are frequently visited, or where the objects spend more time; low values represent the areas that are not frequently visited, or where the objects move faster. It is therefore associated with the concept of utilization distribution, and can be seen as a 3D analogy of the utilization distribution, since the space use intensity is highly positively correlated with the probabilities of occurrence.

The time density of trajectories is tested on a real datasets that represents the movement data of flights at VHHH and VMMC. We expect this method can be used for all types of trajectory analysis. In our test of movement datasets, the time density is suited for trajectories in regions with soft constraints (e.g., seasonal migration of animals in the ocean or sky). It can also be applied to the trajectories in areas with hard constraints (e.g., road traffic or people’s indoor activity).

In addition, time density indicates that the moving objects have visited the place, but the exact visit time (order) is not known. Actually, it presents the resulting view of the time function (Formula (1)), rather than reflecting the time (order) of visit. The space-time cube [2] could, however, embody the chronological order, since it explicitly employs time as its third dimension, representing the moment at which the movement occurs. We suggest two alternative solutions to the problem of the ignored time dimension in the data cube. The first solution is to use the time frame. The position of the moving object can be observed in the movement snapshot at each moment. In other words, we can see the distribution of moving objects at times $t_1$, $t_2$, …, $t_n$. Therefore, by changing the time frame, the sampling point at that moment, or the movement trajectory before that moment, can be displayed. We can manually pull the scroll bar (setting a scroll bar) or play an animation to show the trajectory changes over time in the data cube. Another solution is to introduce a classical space-time cube into our approach, where the space-time density is calculated for each trajectory based on the x-y plane position and time attribute of the trajectory points. Space-time density explicitly employs time as its third dimension, which represents the moment at which the movement occurs. Therefore, we can associate these two visualizations to show space usage at each specific moment. We can also calculate and display the time density directly in the space-time cube. For spatial analyses with explicit temporal patterns as their focus, this method can be used as an alternative.

This work seeks to make a scientific contribution to spatiotemporal analysis of trajectories in large movement datasets. The technology can potentially be used for the analysis of movement data in GIS science and related disciplines. Potential research objects include humans, vehicles, ships, fish or other wild animals, etc. Future work can focus on extending its application to a wider range of scenarios in geographic analysis and GIS science.