{'fraction_non_alphanumeric': 0.03980206540447504, 'fraction_numerical': 0.016320383575116795, 'mean_word_length': 4.10267133671389, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
GraCT: A Grammar-based Compressed Index for Trajectory Data

We introduce a compressed data structure for the storage of free trajectories of moving objects (such as ships and planes) that efficiently supports various spatio-temporal queries. Our structure, dubbed GraCT, stores the absolute positions of all the objects at regular time intervals (snapshots) using a $k^2$-tree, which is a space- and time-efficient version of a region quadtree. Positions between snapshots are represented as logs of relative movements and compressed using Re-Pair, a grammar-based compressor. The nonterminals of this grammar are enhanced with MBR information to enable fast queries. The GraCT structure of a dataset occupies less than the raw data compressed with a powerful traditional compressor such as p7zip. Further, instead of requiring full decompression to access the data like a traditional compressor, GraCT supports direct access to object trajectories or to their position at specific time instants, as well as spatial range and nearest-neighbor queries on time instants and/or time intervals. Compared to traditional methods for storing and indexing spatio-temporal data, GraCT requires two orders of magnitude less space, and is competitive in query times. In particular, thanks to its compressed representation, the GraCT structure may reside in main memory in situations where any classical uncompressed index must resort to disk, thereby being one or two orders of magnitude faster.


Introduction
More than two decades after it emerged, the field of moving object databases is still an active area of research. The data collected from the GPS positions of large sets of cars, ships, planes, smartphones and wearable devices, has lead to a growing interest in applications that exploit trajectory information, for example with data mining purposes.
As these datasets grow in size and the applications managing moving objects become more sophisticated, more space/time-efficient storage techniques are required.
A trajectory is the path followed by a moving object through space as a function of time. Due to storage requirements and the limitations of the devices used to acquire the object positions, the continuous movement of an object is usually approximated with discrete samples of spatio-temporal location points. The more samples taken, the more accurate the trajectory. High sampling rates, however, result in large amounts of data, which increase storage, transmission and processing needs. Even when storage, network and processing capacity grows rapidly, the collected data grows even faster, and thus it is necessary to aim for reduced trajectory representations [49].
Traditional methods for compressing trajectories include line generalization (or simplification) techniques, which keep only some of the trajectory points and discard the rest. This approach results in some loss of information on the real trajectory. A lossless strategy to obtain compression is the use of delta compression, where each new position is stored as the difference with the previous one. This idea exploits the fact that consecutive positions are expected to be closer to each other, and that smaller numbers can be stored using fewer bits. Extracting a whole trajectory with this arrangement is easy. Efficiently accessing the position of an object at a given time, instead, requires sampling some absolute positions at regular time intervals, which introduces a space/time tradeoff.
Another issue is how to index the trajectory data in order to answer queries other than just retrieving a whole trajectory or finding the position of an object at a given time instant. A number of indexes have been proposed since the 90's to handle a rich set of queries on trajectory data. Most indexes were modifications of the R-tree [13] augmented with another dimension to deal with the time. None of those works, however, compress the data. Rather, they are designed to work on disk, which is much slower than the main memory.
With the increasing gap in the access time of main memory versus disk, compressing the trajectories in order to query them in main memory is an attractive option. Some recent proposals following this trend [9,43] build on delta compression, coupled with an encoding that favors small numbers. The optimal codes for delta compression can be obtained with a statistical encoder that exploits frequency bias (typically, smaller numbers are more frequent).
In this article we introduce a new way to compress trajectories and demonstrate that it is much more effective than statistical compression on real-world data. Instead of exploiting the higher frequency of smaller differences, we exploit the fact that, in many applications, trajectories tend to be similar to others, wholly or piecewise. We resort to grammar compression [17], which is a method from the family of dictionary-based compression methods. Grammar compression generates a context-free grammar that generates (only) a given input sequence. Our input sequence is the concatenation of all the trajectories in differential form. Grammar compression exploits the similarities between trajectories: if many trajectories are similar, then a small grammar can be found that generates the whole dataset.
Our index, called Grammar-based Compressed Representation of Trajectories (GraCT), is an in-memory representation that combines grammar compression of trajectories with additional data structures to support efficient queries of various kinds on the dataset, without the need of decompressing it. On one hand, snapshots of the positions of all the objects at regular time intervals are created, and those positions are indexed in a 2 compressed quadtree representation called k 2 -tree [7]. Such spatial index allows us to restrict the set of candidate objects when answering queries. The movements of the objects between two consecutive snapshots are represented in a sequence called the log of movements, which is represented in differential form and compressed with Re-Pair, an effective grammar compression method [19]. To avoid decompressing and sequentially traversing the logs to answer queries, the nonterminals of the grammar are enhanced with additional information, most notably the minimum bounding rectangle (MBR) of all the (relative) positions represented by the nonterminal. MBRs allow us to skip whole nonterminal symbols without the need to decompress them, thereby speeding up the traversal of logs. This arrangement enables GraCT to efficiently answer time-slice, time-interval and k-nearest neighbor queries.
GraCT requires discretizing the data in order to obtain compression. When the data must be used with full precision, GraCT can be used as an in-memory cached index using a sufficient degree of discretization. The candidate set of answers delivered by GraCT are then read from disk and verified.
We evaluate our techniques on three real-life datasets: ships at sea, commercial aircrafts, and city taxis. We show that the grammar compression of GraCT is highly effective, reducing the raw data to 4%-7% of its original size. This is even less than the space achieved by p7zip, a powerful dictionary-based compressor that does not support direct access nor querying. It is also less than half the space used by a baseline index we implement that replaces grammar compression by statistical symbolwise compression. Compared to a traditional method for storing and indexing spatio-temporal data (concretely, the MVR-tree [40]), the compressed representation of GraCT occupies two orders of magnitude less space. In relation to query times, GraCT is up to 4 times faster than the statistically compressed index, except for retrieving whole trajectories, where it is slightly slower. GraCT is also faster than the MVR-tree running in main memory (which, as explained, is feasible for small datasets only) for all queries except on nearest-neighbor queries, on which it can be clearly slower. GraCT is one or two orders of magnitude faster, however, when the MVR-tree cannot be maintained in main memory. Finally, a simple combination with a relational storage of the full-precision trajectory points shows that an in-memory cache based on GraCT is an order of magnitude faster than the MVR-tree operating on disk for the basic time-slice queries.
The paper is structured as follows. Section 2 covers related work. Section 3 introduces the GraCT structure and Section 4 explains in detail the compression of the logs. Sections 5 and 6 describe how the various supported queries are processed in GraCT. Section 7 presents the experimental results. Finally, Section 8 presents our conclusions and directions for future research.

Related Work
This section describes the previous work on indexing and compressing trajectories. Section 2.1 describes the related work on indexing trajectories, Section 2.2 on compressing trajectories, and Section 2.3 on combining both. Readers willing to go directly to the description of our new index may jump directly to Section 2.4, where the basic concepts of compact data structures are covered.

Indexing trajectories
Indexes that handle free trajectories of points in the space are also called spatiotemporal indexes. One type spatio-temporal index builds on a classic multidimensional spatial index, in the form of a temporally augmented R-tree [13]. While the R-tree uses two-dimensional Minimum Bounding Rectangles (MBRs) to enclose the spatial objects stored in the database, the 3DR-tree [42] uses Minimum Bounding Boxes (MBBs), where the third dimension is time, to enclose the segments of the trajectories. The problem is that the three-dimensional space covered by an MBB can be large, even in small segments, resulting in high levels of overlap and limited power of discrimination. The STR-Tree [32] is an extension of the 3DR-tree designed to overcome this problem, by modifying the insertion/split strategy that produces the MBBs. The same work also proposes an index called TB-tree, where several segments of a trajectory are bundled into partial trajectories that are inserted as MBBs of an R-tree.
A second kind of index is a versioned R-tree, which creates an R-tree for each timestamp and a B-tree to select the relevant R-trees. Creating an R-tree for each timestamp requires large amounts of space. To overcome this, instead of storing the complete R-tree for each timestamp, these techniques store only the part of the R-tree that is different from the R-tree in the previous timestamp. Examples of indexes of this type are MR-tree [46], HR-tree [27], HR+-tree [39] and MV3R-tree [40].
A third family of indexes, called grid-based indexes, partition the space and build a temporal index for each of the spatial partitions. The Scalable and Efficient Trajectory Index (SETI) [8] divides the spatial dimension into static, non-overlapping cells. Within each cell, the trajectory segments are indexed by time using an R * -tree. If a segment crosses the boundary between two cells, it is clipped and stored twice in the R * -trees of both partitions. Other examples from this family include Multi Time Split B-Tree [50], Compressed Start-End tree [44], and PIST [3]. More recently, GCOTraj [47] presented a work along the same line, which also divides the space into cells, and all the subtrajectories in a cell are stored in the same disk page. Pages are ordered on the disk following a space filling curve to optimize I/O times.
PA-trees [29] use a completely different approach to avoid MBBs. The index is built with linear and order-k polynomial approximations of the trajectories. Of course, the approximation deviates from the real trajectory, thus the index also stores the maximum deviation with respect to the actual trajectory to avoid false negatives.
A recent active research line is the management of trajectories in the distributed computing framework. PRADASE [22] and CloST [38] are trajectory storage systems built on top of Hadoop and extended with a spatio-temporal index. MD-HBase [31], R-HBase [14], and GeoMesa [15] use distributed key/value stores. TrajSpark [48] is an in-memory system based on Spark and equipped with a two-level spatio-temporal index. While we do not address distributed computing in this article, we remark that, in general, a compressed storage enables the distribution of spatiotemporal data across fewer computers, thereby reducing the communication costs.
The closest predecessor of our GraCT index is SEST-Index [12,45]. This was the first structure considering snapshots and logs of events, but it did not consider compression and was oriented to storing "events" (objects that appear in or disappear from an area).

Compressing trajectories
A simple way to reduce the size of a trajectory is to produce a new trajectory that approximates the original one by selectively removing some of its points. This approach, called trajectory simplification, is based on so-called line generalization techniques. This strategy is typically used as a preprocessing step that is applied before the trajectories are compressed or indexed [49]. Therefore, these methods are compatible with other compression and indexing techniques.
The simplest trajectory simplification method is uniform sampling [33], which selects points at regular intervals (e.g., 5th, 10th, 15th, etc.), and discards the rest. Though simple to apply, the resulting trajectory may lose too much precision.
More evolved algorithms analyze the points in order to identify those containing more information, while discarding more redundant surrounding points. One of the best-known algorithms of this type is Douglas-Peucker [10]. Other methods in this family include top-down time ratio [23], sliding window [16], SQUISH [24], and OPERB [21].
Trajectory compression techniques based on speed and direction include points in the approximated trajectory provided they represent a significant change in the course of the trajectory. Examples of methods from this family include dead reckoning [41], threshold-guided sampling, and STTrace [33].
Finally, some methods use the knowledge of an underlying network to decide which points are retained (e.g., those at bus stops or at street intersections) [35,37].
Independently of whether and how trajectory simplification is applied, the resulting trajectories can be compressed. The typical alternative is delta compression, which stores the first point in the trajectory in absolute form and all the others as the difference with the previous point. This applies to all three components (x, y, and timestamp).
Trajic [30] uses a strategy related to delta compression. It uses a predictor of the next point of a trajectory and an encoder of the residuals, that is, the difference between the prediction and the actual point. If the residual is small, the encoder uses few bits.
Finally, a somewhat related research line faces the compression of traffic data over networks, mainly flow at intersections and average speed, using low-dimensional models such as principal component analysis or low-rank approximation [20,1]. Basically, these techniques approximately represent a matrix as the product of smaller matrices. Similar techniques are used in 3D animated video to compress the trajectories formed by the positions of the vertices of animated synthetic figures [34], and could probably be applied to compress GPS trajectories.
GraCT does not discard points, but discretizes the positions according to an underlying grid, whose size is parameterized by the needs of the applications (e.g., it makes no sense to store ship positions with a precision of centimeters). It also uses a discrete temporal scale, interpolating the measures so that each trajectory has a point at each time instant. GraCT is unique in using grammar-based compression on the trajectories.

Trajectory compression and indexing
There are just a few systems that combine trajectory compression and indexing. TrajStore [9] is similar to SETI [8], as it also slices the trajectories into subtrajectories that fit in the cells into which the space was divided. The difference is that the data of each cell is compressed and that the cells are adapted to the distribution of the data. TrajStore uses delta compression and a cluster-based compression, which groups several similar trajectories and stores only one of them; thus it is a lossy method. In addition, the cells are spatially indexed by a global quadtree, and each cell has a temporal index of the disk pages that make up the cell.
SharkDB [43] is an in-memory trajectory storage system that combines a columnoriented data structure and compression. SharkDB divides the time into frames of a given length (e.g., one minute) and stores one point for each trajectory within that frame. All of the points within the frame are stored as a column of a column-oriented database management system. To reduce data size, it also uses delta compression. SharkDB uses multithreading to speed up query processing.
PRESS [36] is a system that stores compressed trajectories and supports queries, but the trajectories must be paths in an underlying network (e.g., train stations, bus stops). We only consider free trajectories in this article.

Compact data structures and compression methods
Compact data structures [28] are a new family of data structures that combine query support with compressed data representation. They often use space close to that of a plain compressed representation, while supporting queries on the data in time competitive with classical data structures. Our index makes use of various compact data structures and compression methods.

Rank and Select
The following two operations over bit sequences (or bitmaps) are basic components of most compact data structures: rank b (B, p) counts the number of occurrences of bit b in bitmap B up to position p, whereas select b (B, n) returns the position in bitmap B where the n-th occurrence of bit b is located.
Some data structures (see [25], for example) allow these operations to be solved in constant time, using n + o(n) bits of total space (in practice, approximately 5% extra space over the original bitmap).

k 2 -tree
The k 2 -tree is a compact data structure originally designed to represent Web graphs within reduced space, allowing them to be navigated directly in compressed form [7]. In general, a k 2 -tree can be used to represent the adjacency matrix of any graph, as well as binary matrices.
Conceptually, the k 2 -tree is a k 2 -ary tree built from a binary matrix by recursively subdividing the matrix into k 2 submatrices of the same size. It starts by subdividing the original n × n matrix into k 2 submatrices of size n 2 /k 2 . The submatrices are ordered from left to right and from top to bottom. Each submatrix generates a child of the root node whose value is 1 if there is at least one 1 in the cells of that submatrix, and 0 otherwise. The subdivision proceeds recursively for each child with value 1 until it reaches a submatrix full of 0s, or until it reaches the cells of the original matrix (i.e., submatrices of size 1 × 1). Figure 1 shows an example.
Instead of using a pointer-based representation, the k 2 -tree is compactly stored using bitmaps T and L (see Figure 1). T stores all the bits of the k 2 -tree, except those in the last level. The bits are placed according to a level-wise traversal: first the k 2 binary values of the children of the root node, then the values of the second level, and so on. L stores the last level of the tree, consisting of cell values of the original binary matrix.  Any cell, row, column, or region of the matrix can be obtained efficiently via rank and select operations over bitmaps T and L. For example, assuming a value of 1 at position p in T , its k 2 children start at position p children = rank 1 (T, p) · k 2 of T . If the children of a node return a position p children > |T |, the actual values of the cells are retrieved by accessing L[p children − |T |]. Similarly, the parent of position p in T : L is q − (q mod k 2 ), where q = select 1 (T, p/k 2 ), and q mod k 2 indicates the submatrix of p within that of its parent.

Re-Pair
Re-Pair [19] is a grammar-based compression method. Given a sequence of integers I (called terminals), the method proceeds as follows: (1) it obtains the most frequent pair of integers ab in I; (2) it adds rule s → ab to a dictionary R, where s is a new symbol not present in I (called a non-terminal); (3) every occurrence of ab in I is replaced by s, and (4) steps 1-3 are repeated until all pairs in I appear only once (see Figure 2). The resulting sequence after compressing I is called C. Every symbol in C represents a phrase (a sequence of one or more of the integers in I). If the length of the represented phrase is 1, then the phrase consists of an original (terminal) symbol; otherwise, the phrase is represented by a new (non-terminal) symbol. Re-Pair can be implemented in linear time and a phrase may be recursively expanded in optimal time (i.e., proportional to its length).

DACs
Directly Addressable Codes (DACs) [4] are a variable-length encoding scheme for sequences of integers, which support fast direct access to any given position in a sequence. In other words, they allow the i th integer to be decoded without the need to decompress the preceding integers. If the sequence of integers has many small numbers and few large ones, then DACs obtain a very compact representation.
Given a sequence of integers X = x 1 , x 2 , . . . , x n , DACs take the binary representation of that sequence and rearrange it into a level-shaped structure as follows: the first level, B 1 , contains the first (least significant) n 1 bits of the binary representation of each integer. A bitmap C 1 is added to indicate whether the binary representation of each integer requires more than n 1 bits (1) or not (0). In the second level, B 2 stores the next n 2 bits of the integers with a value of 1 in B 1 . A bitmap C 2 marks the integers that need more than n 1 + n 2 bits, and so on. This process is repeated for as many levels as needed. The number of levels and the number n l of bits at each level l, with 1 ≤ l ≤ , is calculated in order to maximize compression. Each value x i is then retrieved using less than rank operations on the bitmaps C l and extracting chunks from the arrays B l .

Permutations
A permutation of the integers [1.
.n] = {1, . . . , n} is a reordering of the values {1, . . . , n}. The two main operations over the reordered sequence are π(i) and π −1 (j). The first yields the number at position i ∈ [1..n] in the sequence, while the second identifies the position in the sequence where j ∈ [1.
The simplest way to store a permutation is an array P [1.
.n], where P [i] = π(i), so that π(i) is answered in constant time. In order to answer π −1 (j) efficiently, we can double the space by storing a second array. We instead use an intermediate option that uses (1 + )n instead of 2n cells and answers π −1 (j) queries in time O(1/ ) [26].

Discretization
GraCT represents moving objects that follow free trajectories in space. The method assumes that the positions of all the objects are synchronized and stored at regular time instants (e.g., every minute). Since GPS devices may report the positions of objects at times that do not match the required time instants, the positions in the desired time instants are obtained from the raw timestamped positions using interpolation.
A raster model is used to represent the space, which is divided into cells of a fixed size, and objects are assumed to fit in one cell. The size of the cells and the length of the period between represented time instants are parameters that can be adapted to the specific domain. The shorter the length of the period is and the smaller the size of the cell is, the more accurate the trajectory data will be, though achieving less compression. This is a tradeoff between space usage and precision in the storage.
Alternatively, we may regard this discretization as a way to obtain a more compact in-memory index, while the full data is maintained on disk and accessed only to filter out false-positives from a usually small set of candidate answers. We will also experiment with this variant of the index, showing that this is better than directly operating on the full data on disk.

Snapshots
Every d time instants, GraCT uses a data structure based on k 2 -trees to represent the absolute positions of all the objects. These data structures are called snapshots. The distance d between snapshots is another parameter of the system, which trades space for query time. Between two consecutive snapshots, the trajectory of each moving object is represented as a log, which is an array of movements represented in differential form.
Let S h denote the snapshot representing the position of all the objects at time instant k. The components of S h are: • The time instant represented by the snapshot.
• A k 2 -tree storing the positions in the space (i.e., the cells of the raster) where there are objects.
• An array of integers perm storing the identifiers of the objects at each cell in the space where there are objects.
• A bitmap Q, which is an auxiliary data structure of perm that is used to determine the correspondence between positions in the space and the positions of perm.
The k 2 -tree represents a binary matrix where a cell set to 1 indicates that the cell contains one or more objects. However, we need to know what objects are in the cell set to 1. Each 1 in the binary matrix corresponds to a bit set to 1 in bitmap L of the k 2 -tree. The list of object identifiers corresponding to each of those bits set to 1 is stored in perm, where the object identifiers are sorted according to their order of appearance in L. This array turns out to be a permutation of all the object identifiers. Bitmap Q is aligned with perm; a 0 at a given position of Q means that the object identifier aligned in perm is the last object in the cell, while a 1 signals that the next object identifier in perm is in the same cell. Figure 3 shows an example. The left side shows the matrix representing the space and the objects placed in the cells (these cells are the same containing a 1 in Figure 1). The four arrays on the right are the actual data structures that represent the snapshot information. Arrays T and L are the same as those in Figure 1 and mark the positions having objects. The objects in those positions are stored in arrays Q and perm. For instance, the object identifiers corresponding to the second 1 in L (which is at position 5 in L and corresponds to the cell at position (9,5)) are stored starting at position 3 in perm, since the first 0 (signaling the end of the entries corresponding to the first 1 in L) is at position 2. In order to find out how many objects there are in the corresponding cell, we access Q starting at position 3 and search for the first 0 after that position, which is at position 4; this shows that there are two objects in the inspected cell. By accessing positions 3 and 4 in perm, we obtain the object identifiers 4 and 5. The object identifiers corresponding to the third 1 in L start at position 5 in perm, and so on.
With these structures used to represent the absolute positions of all moving objects at the snapshots, we can answer two types of queries in time O(log k s + output), where s × s is the grid size and output is the size of the query output: • Find the objects in a given cell : First, we traverse the tree from the root until we reach position n in L corresponding to that cell. Next, we count the number of 1s in the array of leaves L until position n; this gives us the number of leaves with objects up to the n-th leaf, x = rank 1 (L, n). Then we calculate the position of the (x − 1)-th 0 in Q, which indicates the last bit of the previous leaf with objects, and we add 1 to obtain the first position in our leaf. That is, p = select 0 (Q, x − 1) + 1 is the position in perm of the first object identifier corresponding to the searched position. From p, we read all the object identifiers aligned with 1s in Q until we reach a 0, which signals the last object identifier in that leaf.
• Find the position of a given object in the space. First, by using the π −1 operation (see Section 2.4.5), we obtain the position k in perm of the searched object. Next, we have to find the leaf in L corresponding to the k-th position of perm. For this, we calculate the number of leaves with objects before the object in position k of perm; that is, we count the number of 0s until the position before k, y = rank 0 (Q, k − 1). We then find in L the position of the (y + 1)-th 1; that is, select 1 (L, y + 1). With that position in L, we can traverse the k 2 -tree upwards in order to obtain the position in the space of that cell, and thus the position of the object.

Log of relative movements
The movements of objects between two consecutive snapshots are recorded in sequences of consecutive object positions called logs, with one separate sequence per object. Logs actually represent the positions in relative form with respect to the previous one, where the first absolute position can be obtained from the preceding snapshot. Therefore, logs can be thought of as the sequences of consecutive movements of the objects along time. As we will see soon, this model must be extended to consider cases where objects appear and disappear.
Objects change their positions along the two Cartesian axes, so every movement in the log can be described with two integers. In order to save space, our method encodes the two values using a single positive integer. The cells around the actual position of an object are enumerated following a spiral in which the origin is the initial position of the object, as shown in Figure 4 (left). As an example, assume that an object moves one cell to the East and one cell to the North with respect to the previous known position. Instead of encoding the movement as the pair (1,1), it is encoded as an 8. Figure 4 (right) shows the trajectory of an object starting at cell (0,2). Each number indicates a movement between two consecutive time instants. Since most relative movements involve short distances, this technique usually produces a sequence of small numbers.
Between two consecutive snapshots S h and S h+d , each object is represented by a log, L h,h+d (id), where id is the identifier of the object. Figure 5 presents a situation in which there is a snapshot every eight time instants. Time instants t 0 and t 8 are represented with snapshots, while the time instants in the intervals [t 1 , t 8 ] and [t 9 , t 16 ] are represented with logs. Note that, although t 8 is represented with a snapshot, it is also included in the preceding logs. This duplicity allows us to traverse the logs in both directions, which accelerates certain queries, as seen later. The meaning of the other components of the figure is also explained in the sequel.

Managing appearances and disappearances
Because the model does not explicitly store the time instant of each recorded position, a mechanism is required to signal which objects are active at each time instant and to determine to which time instant each log entry corresponds.
For this purpose, the log manages three events: when an object starts emitting its positions, when it stops, and when it stops emitting signals for periods of time. In fact, the three cases are homogeneously managed by the following events, each identified in the log with a special codeword: • Disappearances (D). The occurrence of codeword D in L h,h+d (id) means that object id stopped emitting its position at the corresponding time instant of that portion of the log, and that id did not restart emitting (at least) until the time instant of the next snapshot (S h+d ).
To enable backward traversal of the logs, each disappearance must store the time instant and the absolute coordinates of the object at the moment where it disappears. 11 • Absolute appearance (AA). The occurrence of codeword AA in L h,h+d (id) means that id started to emit its position in the corresponding time instant of that portion of the log, and that id had not emitted its positions (at least) since the previous snapshot (S h ). AA may signal the first appearance of an object or its reappearance after having disappeared in a preceding log.
To enable forward traversal of the logs, an absolute appearance stores the time instant and the absolute coordinates where the object appears. Note that, even if the object is reappearing, it might do so far away from the position where it had disappeared in a previous log.
• Relative disappearance. The occurrence of a codeword indicating a relative disappearance in L h,h+d (id) in an entry corresponding to time instant t i means that id stopped emitting positions at that time instant, and that id restarted emitting its positions in a time instant between t i and the time instant of the next snapshot (S h+d ). In other words, a relative disappearance occurs when the object disappears and reappears within the same log. 12 A relative disappearance must store the number of time instants in which id was not emitting its positions and, depending of the type of relative disappearance, also a movement: -Relative disappearance with movement (RM ) means that id appears in a different position with respect to its last known position. This requires the log to store the relative movement (using the spiral encoding) with respect to that last known position.
-Relative disappearance without movement (RN M ) means that id appears in the same position where the signal was lost, therefore there is no need to store a movement.
The extra information is stored in two auxiliary arrays associated with each log L h,h+d (id): • Array D h,h+d (id) stores: in the case of relative disappearances, the time the disappearance lasted; in the case of an absolute appearance or a disappearance, the absolute time instant of that event.
• Array P h,h+d (id) stores: in the case of relative disappearances with movement, the relative movement with respect to the last known position; in the case of an absolute appearance or a disappearance, the absolute position of that event. Figure 5 shows two snapshots; S 0 is the snapshot of Figure 3. The right-hand side shows the logs representing the movements of objects in the time instants after those snapshots, thus representing the positions of objects between time instants t 0 and t 16 . Object 1 disappears at t 2 and then reappears at t 5 (see L 0,8 (1)). Since the time instants occur between the same pair of snapshots, this is a relative disappearance. The number of time instants that the disappearance lasts is stored in array D 0,8 (1). In the first position of D 0,8 (1) there is a 3, which means that the first disappearance of L 0,8 (1) lasted three time instants. In the same portion of the log, there is another relative disappearance between time instants t 6 and t 7 , therefore in the second entry of array D 0,8 (1) there is a 2, indicating that the disappearance lasted two time instants. The relative disappearance that starts at t 2 is of type RM , so the log must store the position where object 1 reappears. Array P 0,8 (1) stores this position using spiral encoding with respect to the last known position. In the example, a 45 means that object 1 reappeared three cells to the North.
In the entry in L 0,8 (3) corresponding to time instant t 5 , codeword D means that object 3 disappears at that time instant and does not reappear in L 0,8 (3). Therefore, D 0,8 (3) stores a 5 in its first position (there are no prior relative reappearances), indicating that the object disappeared at t 5 , and in P 0,8 (3), the only entry stores the coordinates (7,9), which denote that the object was in that position when it disappeared. Object 3 reappears in L 9,16 (3), in the entry corresponding to t 12 . This reappearance is signaled with 13 codeword AA aligned with time instant t 12 , but this is only to improve the readability of the figure. In fact, codeword AA is the first entry in L 9,16 (3), so, in order to find out to which time instant codeword AA corresponds, D 9,16 (3) stores a 12. Additionally, in order to find out the position of object 3 at t 12 , P 9,16 (3) stores the absolute coordinates of the object at the time instant of reappearance, which in our example are (4,8). Finally, to help manage disappearances and appearances during queries, two additional arrays are stored with each snapshot S h : • Dis, a list of the objects that were active in the previous snapshot and stopped emitting before the time instant represented by S h .
• App, a list of objects that are not present in S h and appear (or reappear) before the next snapshot. Figure 5 shows that S 8 .Dis contains objects 3, 5 and 6, meaning that they were active in their portions of the log immediately before S 8 , but are missing at time instant t 8 . Similarly, S 8 .App contains objects 3, 5, 6 and 7, indicating that they are missing at t 8 , but appear in the portion of the log immediately following t 8 .

Compressing the Log
In many applications, objects spend most of their time either stopped or moving along a specific course at a fixed speed. This generates long sections of the log with numbers representing the same or contiguous values of the spiral. For example, the moving object in Figure 6 follows a NE trajectory, moving one or two cells in the time elapsed between two consecutive time instants. Its log represents the series of relative movements 2,9,2,9,8,7,9,8,7,9,9,2,9,2,9,9,9 (see array I in Figure 2). These series of similar movements are compressed very efficiently using a grammar compressor.
In order to query trajectories without completely decompressing them, the grammar rules in GraCT not only include the symbols to be replaced, but are also enriched with additional information. Specifically, each rule in R will have the following information: s → a, b, #t, x, y, M BR, where: 14 • s → ab is the normal rule of Re-Pair, • #t is the number of time instants covered by the rule, • (x, y) are the relative coordinates of the final position of the object after the application of the rule (that is, the total displacement in both coordinates), and • M BR is the minimum bounding rectangle enclosing the movements of the rule. MBR is represented as (x1, y1, x2, y2), where (x1, y1) corresponds to the bottomleft corner and (x2, y2) to the top-right corner.
For example, the rules in Figure 2 are enriched as follows. The first rule of R is W → 2, 9, 2, (3, 0), (0, −1, 3, 0): 2 and 9 are the substituted symbols; the next 2 indicates that the rule represents a sequence of two movements; (3,0) indicates the position of the object after the application of the rule if we start at (0, 0), and the last four values are two corners (bottom-left and top-right) defining a rectangle that encloses all the movements encoded by the rule.
With this additional information, most of the non-terminal symbols of a compressed log file L do not need to be decompressed in order to obtain the position of an object at a given time instant between two snapshots. Consider, for example, the compressed log L = Z, Y, 9, Y, X, Z, X of Figure 6 (this is the array C of Figure 2). Assume that we want to find the position of the object at the 5th time instant, that is, (7,2). From the preceding snapshot we obtain that the absolute position of the object at the beginning of the log is (0,1). Next, we inspect the log L from the beginning. The first value is a Z. The enriched rule indicates that this symbol represents four time instants, after which the object is displaced six columns to the East while remaining in the same row. Since it started at (0, 1) the object will be at position (6, 1) at t 4 . Since our target time instant is t 5 > t 4 , we can skip the whole nonterminal Z. The next symbol is Y , which lasts two time instants. This would take us to t 6 , which surpasses our target time instant. Therefore, and only in the last step of the search, we have to decompress the rule and process its components: Y → 8 7. The 8 is a terminal symbol that lasts one time instant, and thus it is enough to reach our target time instant. An 8 moves one column to the East and one to the North, which applied to the current position (6, 1), takes us to the correct position (7,2).
The M BR component also helps speed up other queries, as seen later. These additional elements that enrich the rules are compressed with DACs. To obtain better compression, the times of all of the rules are compressed using one DAC, separately from the three pairs of coordinates of all of the rules, which are compressed using a separate DAC. Figure 7 shows the resulting data structures that make up the GraCT index for our running example in Figure 5. Figure 8 shows a flow chart that summarizes the whole process to create a GraCT data structure.

Query Preliminaries
This section presents a series of definitions for later use. • S (x,y) , an absolute position (x, y) obtained from a snapshot.
• NT/T, a grammar symbol encoding the movements of O, where NT denotes a nonterminal and T denotes a terminal.
• D t,(x,y) , a disappearance code, where t is the time instant of the disappearance and (x, y) is the object position at that time instant.
• RN M # , a reappearance without movement code, where # is the duration of the disappearance.
• RM #,value , a reappearance with movement code, where # is the duration of the disappearance and value is the spiral code representing the movement.
• AA t,(x,y) , an absolute appearance code, where t is the time instant of the appearance and (x, y) is the object position at that time instant.
If object O is missing at t b , then m 1 will refer to the first time instant after t b containing information about the position of O. In the same way, if the object is missing 16  For example, in our running example in Figures 5 and 7, for object O 6 , the compressed trajectory between 0 and 13 is CT O6 (0, 13) = (S (10,13) , 1, V, D 4, (10,15) , AA 12, (3,6) , Y ). Note that the last code (Y ) covers time instants 13 and 14, since it is a non-terminal symbol. In contrast, in CT O5 (11, 15) = (AA 11, (7,2) , W, D 13, (8,4) ). The last code (D 13,(8,4) ) corresponds to time instant 13, since there are no positions after that time instant in the requested time interval.
, p c is the absolute position of object O at time instant t c , and t b ≤ t c ≤ t e , then operation moveJ(p c , t c , t e , m i ) returns either: • t nc , p nc , where p nc is the position resulting from the application of the movements encoded in m i to p c , and t nc ≤ t e is the time instant of that position, or • t e , p e , where p e is the position at t e , if m i encodes movements surpassing t e .
The name moveJ is the abbreviation of "move jump". If m i is a non-terminal code, it uses m i .(x, y), the relative coordinates of the enriched information included in the definition of the non-terminal symbol m i , which, when applied to the previous position, gives the new position of the object after the application of all the movements encoded by m i . For example, moveJ((9, 5), t 1 , t 3 , Z) uses the relative (x, y) coordinates of the definition of Z → 4, 5, 2, (−2, −1), (−2, −1, 0, 0), that is (−2, −1), which, when applied to (9,5), obtain the new position (7,4). Note that the process never obtains the position at t 2 but jumps from t 1 , (9, 5) to t 3 , (7,4) .
Note also that if we issue moveJ((9, 5), t 1 , t 2 , Z), where the application of Z would produce the position at t 3 > t 2 = t e , the operation has to decompress Z into its components (4 and 5) and apply them one by one until it reaches t e . In this case, therefore, only code 4 is applied, to move the object to (8,4). The decompression of non-terminal codes in moveJ only occurs in this situation.
In this example we fully decompress the nonterminal that exceeds t e , but in general we continue advancing recursively on the terminals and nonterminals of the right-hand 17 compressed trajectory of O between time instants t b and te moveJ(pc, tc, te, mi) returns the pair tnc, pnc resulting from applying the movements in mi to position pc, but if that application produces a position surpassing te, then the operation returns te, pe . RmoveJ(pc, t b , tc, mi) returns the pair tnc, pnc resulting from undoing the movements in mi to position pc, but if that application produces a position previous to t b , then the operation returns t b , p b . moveS(pc, tc, te, mi) returns ( (tc+1, pc+1 , . . . , t h , p h ), the list of time instants/positions resulting from applying mi to pc, except if mi encodes movements that surpass te, in which case h = e. ERr(t b , te) expansion of region r from t b to te using the maximum speed expanded region of object O from t b to te using the maximum speed pipj distance in cells of the shortest trajectory from point pi to point pj pr distance in cells of the shortest trajectory from point p to region r side of its rule. That is, if the application of the next nonterminal Z would exceed t e , and the rule associated with Z is Z → Z 1 Z 2 , then we replace Z by Z 1 Z 2 and try again. This means that, if applying Z 1 also exceeds t e we will continue recursively with Z 1 , and otherwise we will process the whole Z 1 in constant time and continue recursively with This is a reverse version of moveJ, and can be implemented within the same time complexity.
Definition 5.4. If m i is an element of CT O (t b , t e ), p c is the absolute position of object O at time instant t c , and t b ≤ t c ≤ t e , then moveS(p c , t c , t e , m i ) returns the list of pairs ( t c+1 , p c+1 , . . . , t h , p h ) where t h ≤ t e , which is the list of time instants/positions resulting from applying m i to p c . If m i encodes movements that surpass t e , however, the decompression ends at t e .
The name moveS is the abbreviation of "move step by step". Unlike moveJ, if m i is a non-terminal code, moveS decompresses the symbol until it obtains terminal codes, then applies those terminal symbols to p c and returns the positions of each symbol.
For example, moveS((9, 5), t 1 , t 3 , Z) decompresses Z to obtain the terminal codes 4 and 5, and thus returns t 2 , (8, 4) , t 3 , (7, 4) . In contrast to moveJ, moveS does obtain the position at t 2 . This operation can be implemented in time O(h − c + log t) on a balanced grammar.
Definition 5.5. Let r = [x 1 , y 1 ] × [x 2 , y 2 ] be a rectangular region in the two-dimensional space. Let t b and t e be two time instants, where t b < t e , and let M sp be the maximum speed of the dataset. Then ER r (t b , t e ) denotes the expanded region of r from t b to t e using M sp : The maximum speed of the dataset is the maximum speed of any object during the period of time covered by the dataset. ER r (t b , t e ) is the rectangle resulting from expanding r in all directions (North, South, West and East) so that any object within r at t b is bound to be within ER r (t b , t e ) at t e . In other words, to obtain ER r (t b , t e ), r is enlarged by M sp · (t e − t b ) cells in all four directions.
Definition 5.6. Let t b < t e be two time instants, and let M sp be the maximum speed of the dataset. Then ER O (t b , t e ) denotes the expanded region of object O from t b to t e using the maximum speed That is, ER O (t b , t e ) delimits the area where object O is bounded to belong at t e even if running at maximum speed M sp in any direction, starting at time instant t b .
Definition 5.7. Given a point p i and a point p j , the distance between p i and p j (p i p j ) is the length of the shortest trajectory (in number of cells) from p i to p j , using Euclidean distance.
Definition 5.8. Given a point p and a region r, the distance between r and p (pr) is the length of the shortest trajectory (in number of cells) between r = [x 1 , x 2 ] × [y 1 , y 2 ] and p = (x, y). This is zero if Otherwise, it is the minimum Euclidean distance from p to the four corners of r. Table 1 summarizes these definitions.

Supporting Queries
We now describe how each of the supported queries are answered using the structures of GraCT.

Obtain the position of an object in a given time instant
Algorithm 1 shows the pseudocode to answer this query. The algorithm first identifies a time instant close to t q where the position of the object is represented in absolute form and, from that time instant, applies the relative movements in order to obtain the absolute position of the object at time instant t q .
The first option to find that time instant is the snapshot S h that is nearest to the queried time instant (line 1 of the algorithm). If t q coincides with the time instant of S h , then the position of the object in the snapshot is returned (line 28; the position could be null if the object is not present in S h ). The other possible cases are handled in lines 4-15 (when t q follows the snapshot) and 16-27 (when t q precedes the snapshot).
In the first case, t c stores the time instant and p c the absolute position at any step of the log traversal. In line 4, t c and p c store the information of the object at S h . If the object is missing at S h (i.e., p c = null), line 6 checks if O id appears in the immediate next portion of the log (L h,h+d (O id )), as signaled by the presence of O id in S h .App. In line 7, the algorithm obtains the time instant of the appearance of O id from the first entry of D h,h+d (O id ). If the appearance is prior to or occurs exactly at t q , then the algorithm obtains the absolute position from P h,h+d (O id ). Otherwise, the object does not have a position at t q and the algorithm returns null.
Once we have the absolute position of the object, lines 12-15 traverse the compressed trajectory using moveJ operations from time instant t c (either the closest snapshot or an absolute appearance) to t q . If the algorithm is able to reach t q , the position is returned.
Lines 16-27 perform the same process as lines 4-15 but in reverse, for the case where t q occurs before the snapshot.
If there are t time instants, n objects, the grid is of size s × s, k is the arity of our k 2 -tree, d is the number of time instants between two snapshots, and our grammar is balanced, then this query takes time O(d + log t + log k s + log n). If the flow reaches line 13, the algorithm has found a suitable timestamp t c ≤ t b for which the absolute position of the object (p c ) is known. If t c < t b , lines 13-17 skip entries of the log(s) until t c = t b and p c is the position of the object at t b . Finally, lines 18-26 traverse from instants t b to t e , this time collecting all the positions in a sequence Answer, which is finally returned.

Obtain the trajectory of an object between two time instants
Compared to the complexity of Algorithm 1, we have an additional cost of O((t e − t b + 1)(1 + (log n + log k s)/d)), which can be made the optimal O(t e − t b + 1) by choosing a sufficiently large d = Ω(log n + log k s). This is, however, a worst-case complexity that ignores the fact that the log is compressed: on average, the term O((log n + log k s)(t e − t b + 1)/d)) must be multiplied by the compression ratio.

Time slice query
Given a time instant t q and a window rectangle r = [x 1 , x 2 ] × [y 1 , y 2 ] of the space, this query returns the objects lying within r at time t q , and their positions. Algorithm 3 shows the pseudocode.
The process starts by obtaining the snapshot nearest to t q , S h . The algorithm considers three cases. Firstly, if t q is the time instant of S h , the algorithm just needs to access the snapshot in order to obtain the objects within r. The second case occurs when t q is between the snapshots S h and S h+d . In this case, the algorithm has to follow the trajectory of the objects from S h until time t q , to see if they end up within r. Instead of following all of the objects in the snapshot, the algorithm tracks only those within ER r (S h .time, t q ), ignoring objects that have no chance of being in the answer. The algorithm must also include in the tracked objects those that are not present in S h but appear at a later time instant, not exceeding t q , and are capable of reaching r at t q .
Once this list of candidates has been obtained, the algorithm follows the log to obtain the trajectory of the tracked objects. However, the algorithm prunes the tracking further as it processes the log: a candidate object may follow a direction that takes it away from region r, so the algorithm rechecks the condition after applying each code of the compressed trajectory, discarding the objects as soon as they become incapable of reaching r at time t q .
The third case is analogous, and occurs when t q is between S h−d and S h . To illustrate the idea using the running example shown in Figures 5 and 7, Figure 9 outlines the process of a time slice query where t q = t 10 and r = [7, 10] × [3,4]. Since the target time instant is t 10 , the closest snapshot is S 8 . Figure 9(a) shows the target region r and the extended region ER r (t 8 , t 10 ) superimposed. In this example, M sp = 1 (that is, an object can only move one position per time instant), so ER r (t 8 , t 10 ) is r enlarged by two cells in all directions.
Line 1 obtains the snapshot closest to t 10 (S 8 ). Since t q is not the time instant of the snapshot, the algorithm retrieves in line 4 all the objects in the extended region ER r (S h .time, t q ) = [5,12] × [1,6], which are stored in the variable Candidates; in our example, these are the objects O 2 and O 1 . Next, in lines 7-11, for each object in S h .App, the algorithm obtains the time and position of its appearance (t c and p c ); if that time is t q or earlier, and that position is within ER r (t c , t 10 ), then the object is added to the candidates to be tracked. In our example, objects O 3 , O 6 and O 7 appear after t 10 , therefore they are discarded. Instead, O 5 appears at t 9 , so it is added to the list of candidates, because it is also within ER r (t 9 , t 10 ) (see Figure 9(b)). The resulting list of candidates in our example comprises O 1 , O 2 and O 5 .
In lines 12-20, the algorithm follows the movements of each candidate object using the log, until it reaches t q or until the object can be discarded. There are two ways to preempt the log traversal of an object. One is that we reach a nonterminal that already covers t q while its MBR is disjoint with r; this is discarded in lines 16-17. If this does not happen, the log entry is processed in line 18 and we consider the second condition: at this point, we can ensure that the object will not reach the region r. This is verified in line 19. Finally, line 20 includes the object in the result if we reach t q and the p c is in r. Lines 22-35 process the logs backwards from the following snapshot.
Assuming we first process O 1 , Step 1 of Table 2 shows the initial state of the algorithm. The algorithm traverses the trajectory of O 1 and checks if it can reach r at t q after applying each m i of CT O1 (t 9 , t 10 ) = 7, V . The state after applying the first code (7) is shown in Step 2. The algorithm checks if the new position, (9, 6), lies inside the extended region corresponding to t c = t 9 . As shown in Figure 9(b), O 1 is outside the region, so the object is discarded.
The algorithm now processes O 2 . Its initial position, p c = (10, 3), and its time instant, 22  t c = t 8 , are obtained from the candidate list; this state is shown in Step 3. The state of Step 4 is obtained by applying the first code of CT O2 (t 9 , t 10 ) = 5, W . Since (9, 3) is within ER r (t 9 , t 10 ), as shown by Figure 9(b), the algorithm applies the second code of the trajectory: m i = W (a non-terminal symbol). The algorithm checks that the application of the non-terminal symbol reaches or exceeds t q (t c + value.#t ≥ t q ). In our example, t c = t 9 and W.#t = 2, so the rule would lead us to t 11 . The MBR of W , however, is [0, 1] × [1,2], which added to the current position (9, 3) is [9, 10] × [4,5]. This is not disjoint with r, so O 2 cannot be discarded. The algorithm then has to decompress W , to determine if the object is within r after applying the part of W that spans exactly up to t q . Since moveJ((9, 3), t 9 , t 10 , W ) yields that the position of O 2 at t 10 is (9, 4), as shown in Step 5, the object O 2 is reported. Finally, Step 6 shows that we start processing O 5 at t 9 , where it appears. The case of this object is similar to O 2 , and it is also reported in Step 7. Table 2: Trace of Algorithm 3 with the example of Figure 9. An average-case analysis of the time complexity of this query can be made assuming that the objects distribute uniformly in the space and without taking advantage of the compression of the log. Let the query region r be of dimensions x × y , and the query instant be at distance t from its closest snapshot. The spatial query on the k 2 -tree, for a range of size ( x + 2t) × ( y + 2t) and retrieving c candidates, takes O( x + y + t + (c + 1)k log k s) time [28, Sec. 10.2.1]. On average, c = ( x +2t)( y +2t) x y · o, where o is the actual output size. Further, we traverse all the objects that are within ranges of size ( x + 2t) × ( y + 2t) (at the snapshot time) decreasing by one until size x × y (t instants later or earlier). The sum of all those costs is, on average, o· t i=0 ( x +2i)( y +2i) x y = O(o·(t+t 2 ( x + y +t)/( x y ))). In total, the query has a fixed cost of O( x + y +t+k log k s) plus an average cost of O(k log k s + t + (t( x + y ) log k s + t 2 ( x + y + t + k log k s))/( x y )) per object reported.
The analysis shows that, for this query, it is especially worthwhile to start from the closest snapshot, be it preceding or following t q , so that t is on average d/4. In addition to having to process fewer log entries, the distance t between S h .time and t q directly impacts the size of the snapshot area from where we have to collect the candidate objects, and therefore the average number of candidates. The analysis also shows that the effect of t is reduced when the query area is larger.

Time interval query
Time-interval queries return the objects that were within r at any time instant during time interval [t b , t e ]. An easy solution to this query is to use a time-slice algorithm 24 ProcessPortion is a modification of the time-slice algorithm. In lines 3-10, it collects the candidates to verify from the snapshot and the appearances that follow, avoiding the objects already reported. Lines 11-26 treat the trajectory of each candidate. A novelty with respect to time-slice queries is in lines 18-19: since we only report the objects and not their positions, we can directly report an object if the MBR of its next nonterminal is fully contained in the query region r. Otherwise, we must expand the nonterminal (lines 20-21) unless its MBR is disjoint with r, in which case we can process it as a whole, even if it exceeds t last . Note that, since we manage the expansion of the nonterminals in line 21, line 23 always treats a complete symbol of the trajectory, knowing that it cannot exceed t last (or that, if it exceeds t last , it is because its MBR is disjoint with r). Finally, lines 24-26 process a single step, adding the objects found inside r and skipping them when they are shown to be unable to reach r within the time left until t last .
We have omitted other possible speedups for simplicity. For example, let r = [x 1 , x 2 ]× [y 1 , y 2 ] and pc + value.M BR = [x 1 , x 2 ] × [y 1 , y 2 ]. Then, in line 18, we can also directly add O j to the answers if t c + value.#t ≤ t e and, either x 1 ≤ x 1 ≤ x 2 ≤ x 2 and [y 1 , y 2 ] ∩ [y 1 , y 2 ] = ∅, or if y 1 ≤ y 1 ≤ y 2 ≤ y 2 and [x 1 , x 2 ] ∩ [x 1 , x 2 ] = ∅. In both cases, we know that O j must be inside r at some moment not exceeding t e .
The analysis for these queries is similar to that for time-slice queries, now taking t as the distance between t e and the closest snapshot preceding t b (and ignoring the help of intermediate snapshots). While omitted for simplicity, we can also process this query from the snapshot following t e with the aim of reducing t.

Nearest neighbor queries
Nearest neighbor (knn) queries return the K objects nearest to a given query point p q at a given time instant t q . To solve this query, we obtain the candidate objects from the snapshot (S h ) that is nearest to t q , follow their trajectory up to t q , and retain the K candidates that are closest to p q at time t q . The problem is that, unlike in time-slice queries, we cannot easily determine which objects cannot be part of the solution. This condition changes dynamically as we find better and better candidates for the answer set: given the K closest objects to p q found up to now, the distance from p q to the K-th closest object in t q can be used to compute an upper bound to the distance between 25 Algorithm 5: ProcessPortion(r, t b , t e , S h , Answer) 1 tc ← S h .time; 2 t last ← min(te, S h+d .time); 3 Candidates ← ; // The candidates to be considered, as a list 4 foreach O j , tc, pc ∈ S h .f indObjectsInRegion(ERr(tc, t last )) do 5 if O j ∈ Answer then Candidates.add( O j , tc, pc ); if tc ≤ t last and pc ∈ ERr(tc, t last ) and O j ∈ Answer then  else if tc < t last and pc ∈ ERr(tc, t last ) then break; // Skip candidate ; a candidate and p q in S h , so that if the candidate is farther away, it has no chance of getting closer than the current K-th closest object at t q .
Since this upper bound becomes more restrictive as we improve the answer with closer objects, it is important to obtain objects close to p q as soon as possible, as this allows us to discard other objects earlier. For this purpose, we will use the hierarchical space partitioning induced by the k 2 -tree of S h , and will traverse the regions prioritized by their potential relevance to the query (i.e., those closer to t q will be processed first).
The traversal is carried out with a priority queue, which stores k 2 -tree regions r prioritized by increasing values of d(r) = p q r. Initially the whole grid, corresponding to the k 2 -tree root, is inserted in the priority queue, and then we iteratively extract the most promising node v (i.e., covering the region r with minimum d(r) value). Let {r 0 , r 1 , ..., r k 2 −1 } be the subregions of r induced by the children of v. The nonempty subregions r i corresponding to children v i of v are then reinserted in the queue with priority d(r i ) = p q r i . When the children of the region r that we extract are leaves of the k 2 -tree, the corresponding objects O j are processed in a different way.
Those candidate objects O j we find are inserted in a second priority queue, using a similar ordering on p q p c , where p c is the position of O j ; we also record the time instant t c of the snapshot S h associated with the objects. In this second queue we initially accept every object that does not disappear in its log by time instant t q , until we have K elements. From this point, objects O j are inserted only if at time t q they can get closer to p q than the K-th object closest to p q in the second queue. More precisely, let O, which was known to be at position p in time instant t, be the K-th object closest to 26 Let O, t, p, d be the K-th smallest distance in Q C ; 32 dmax ← pqp + Msp · (tq − t) 33 return Q C p q in the priority queue (this is easily maintained by keeping the K closest objects in a separate queue prioritized by decreasing distance). Then the maximum distance to p q it can reach at time t q is d max = p q p+M sp ·(t q −t). Let O j be at position p c in time instant t c . Then, at instant t q , O j cannot be closer to p q than d min (O j ) = p q p c − M sp · (t q − t c ). Therefore, we insert O j as a candidate only if d min (O j ) < d max . When we populate the queue of candidates, we must also consider the objects O j in S h .App appearing at times t c ≤ t q (and not disappearing before t q ). We prioritize the candidates O j using d min (O j ). This condition to insert an object as a candidate gives also a criterion to stop the processing of the priority queue of regions. If we extract region r from the queue, no object from r can be, at time instant t q , at a distance to p q below d min (r) = p q r − M sp · (t q −t c ). Therefore, if d min (r) ≥ d max , there is no need to consider any point from region r, nor from any region yet to be extracted from the first priority queue. At that point we can discard the first queue and start processing our queue of candidates. Algorithm 6 gives the code (for simplicity, we consider only the nearest snapshot preceding t q ).
Let us assume that k = 2, M sp = 1, K = 1, t c = t 8 , t q = t 9 , p q = (10, 0), and the grid S 8 of Figure 10. We number the regions with a sequence of subindices in {0, 1, 2, 3} considering the consecutive levels. The children of the root at r 0 to r 3 , the children of r 0 are r 00 to r 03 , and so on. From the root, we insert r 2 and r 3 in the queue (r 0 and r 1 are empty). We then extract r 3 , the closest one to p q (indeed, containing it), and reinsert r 30 and r 32 . The next most promising region is r 32 , from whose children we only reinsert r 321 , the nonempty one. From the current set, {r 2 , r 30 , r 321 }, the most promising one is r 321 , at distance p q r 321 = 2. So we extract it and insert object O 2 as our first candidate, at distance p q O 2 = 3. Thus, at time t 9 , O 2 could be at maximum distance p q O 2 = 4. Since the next most promising area, r 2 , is at distance p q r 2 = 3, it is still worth considering, so we reinsert its nonempty child r 23 . This is still more promising than r 30 , so we extract it and reinsert r 230 , its nonempty child. Now our candidate regions are {r 30 , r 230 }, and the most promising is r 30 . We extract it and reinsert r 302 , its nonempty child, which becomes the first region to extract. Since p q r 302 = √ 17 ≈ 4.1, a point inside r 302 might reach distance ≈ 3.1 to p q at t 9 , and thus it might still be interesting. We then extract it and obtain O 1 , which is at distance p q O 1 = √ 26 ≈ 5.1, and can then get as close as distance ≈ 4.1 in t 9 . Thus, there is no point in considering O 1 , because it cannot get closer to p q than O 2 . Similarly, there is no point in considering the remaining region r 230 (nor any other region that would have remained to process), because it is at distance p q r 230 = √ 29 ≈ 5.4. In this case, we only have one candidate. From the priority queue of candidates, the algorithm extracts the first one (with position p c , at time t c , with minimum d min ) and calculates the next position of the object using moveJ. The object is then reinserted in the queue with its new position p c and time instant t c , and with d min recomputed accordingly. As we iterate, the candidates advance in their log L h,h+d , with the most promising objects advancing sooner towards t q . At some point, we will start extracting objects with time t c = t q . For those objects, d min is their exact distance to p q at time t q . We store them in a third priority queue, of maximum size K and sorted by decreasing d min . When we have K answers in this third queue, the insertion of a new result may displace the currently K-th result (the first in the priority queue). Further, we can use the value d min (O) of the K-th result O Algorithm 7: Knn(K, p q , t q ,) 1 S h ← f indN earestP reviousSnapshot(tq); 2 Q C ← ObtainCandidates(S h , K, tq, pq); 3 Q F ← ∅; // Priority queue of pairs (obj,dist) sorted by decreasing dist, max size K 4 dmax ← +∞; // Infinite until obtaining K results 5 if S h .time = tq then return first K objects of Q C ; 6 while Q C is not empty do Algorithm 7 shows the pseudocode of the knn query. To illustrate it, we use the example from Figures 5 and 7 and obtain the nearest object (K = 1) with respect to location p q = (10, 0) at t q = t 9 . After obtaining the previous closest snapshot to t 9 (S 8 ), Algorithm 7 issues the call ObtainCandidates(S 8 , 1, t 9 ) to obtain the candidates from the snapshot. As we have already seen, the traversal of regions yields only one candidate, O 2 , with minimum distance d min = p q O 2 − M sp · (t q − t c ) = 2. Therefore, the queue of candidates contains only Q C = { O 2 , t 8 , (10, 3), 2 }. However, objects O 3 , O 5 , O 6 , and O 7 appear in the log after S 8 . From those, only O 5 appears before or at t 9 , and thus it can be added to Q F . Indeed, O 5 appears at t 9 , in position (7, 2). Its minimum (and exact) distance to p q in t q is d min = p q (7, 2) = √ 13 ≈ 3.6. Since d max = 4, this object cannot yet be discarded, so the final queue of candidates is Q C = { O 2 , t 8 , (10, 3), 2 , O 5 , t 9 , (7, 2), 3.6 }.
Once Algorithm 7 has all the candidates, it starts processing them in lines 6-16. It first extracts O 2 , t 8 , (10, 3), 2 . Since t 8 is not the target time, the algorithm considers CT O2 (t 9 , t 9 ) = 5 and applies moveJ on the first element, 5. This converts t 8 into t 9 and (10, 3) into (9, 3). It then recomputes d min = √ 10 ≈ 3.2 and reinserts O 2 into Q C with these new values. The queue of candidates now contains Q C = { O 2 , t 9 , (9, 3), 3.2 , O 5 , t 9 , (7, 2), 3.6 }. The next best candidate is again O 2 , but now its time is t 9 , the same as the query, and this means that we have our first result, Since now |Q F | ≥ K, we can compute d max = 3.2. In the next step, when we extract O 5 from Q C , its d min value is 3.6, larger than d max , and thus disregard it. We have exhausted Q C and the answer is in Q F .

Experimental Evaluation
In this section we experimentally evaluate the space and time performance of GraCT and compare it with related work. It has been very hard to obtain implementations of previous indexes. Further, compression-only systems developed for trajectories like Trajic [30] and SQUISH [24] did not compress well on our data, reducing the space of a plain text representation, but not the space of a binary representation of the integer coordinates. We have, nevertheless, compared GraCT with two baselines: scdcCT: This structure is similar to GraCT with a different representation of the log structure. The sequence of spiral codes that differentially encode the movements is compressed with a statistical zero-order compressor, namely byte-oriented (s, c)-Dense Codes [5] (SCDC) optimized to s = 183 and c = 73. SCDC combines good compression ratios with fast decoding, thereby supporting fast scanning of the logs. We use scdcCT as a representative of the various systems that use statistical compression, instead of the grammar-based compression used by GraCT. We aim to demonstrate that grammar compression is more effective, and that the ability of GraCT to skip over whole nonterminal symbols makes it faster than statistical encoders that, no matter how fast, must process each movement individually. Note that sampling absolute values at regular intervals to speed up the processing of logs in scdcCT is equivalent in space to shortening the sampling period d of the snapshots.
MVR-tree: This is the only classic spatio-temporal structure for which we could obtain code. The complete structure is the MV3R-tree [40], which combines an MVRtree with an auxiliary 3DR-tree, but only the MVR-tree part was available. The 3DR-tree is used to speed up the most basic queries, like obtaining the position of an object at a given time instant or retrieving its trajectory. We aim to show that a classic structure uses orders of magnitude more space than GraCT, without being much faster (actually, the MVR-tree is slower for some queries). We used the MVR-tree implementation provided by the spatialindex library 1 . This is a diskbased structure, but we make it run in main memory for a direct comparison with GraCT. We also compare both structures on disk, as explained later.
We used an implementation for supporting rank operations that adds 6.25% space on top of the bit sequence to accelerate queries. The sample value 1/ of the permutations was set to 5, which increases the size of the data structure by 20% but yields fast π −1 queries. Our grammar compressor was the Re-Pair implementation provided by G. Navarro 3 , using the balanced version. The extra information associated with the Re-Pair nonterminals was encoded using DACs with an unlimited number of levels and without a predefined chunk size. Arrays P and D, in contrast, are represented using DACs configured to use chunks of one byte and with a maximum of two levels.

Datasets
We used real-world and pseudo-real data from three sources.
• Ships, a real dataset storing the movements of 3,654 boats sailing in UTM Zone 10 over one month in 2014, obtained from MarineCadastre 4 .
• Planes, a real dataset composed by trajectories of aircrafts from 30 different airlines and flights between 30 European airports. This data was obtained from the OpenSky Network 5 ; we discarded the altitude component and retained latitude and longitude.
• Taxis, a larger synthetic dataset simulating trajectories of taxis in New York City during 2013. We download the data from NYC Taxis: A Day in the life 6 . Since this includes only the start and end location of each taxi trip, we computed the trajectory as the shortest route between these two points.
GraCT stores positions in a discrete grid, therefore each position emitted by an object is discretized into a matrix. We chose cell sizes that seem appropriate in each application, in the sense that more detail is probably useless: 50 × 50 meters in Ships, 5000 × 5000 meters in Planes and 10 × 10 meters in Taxis. In addition, our structure deals with object positions at regular timestamps, but in the real datasets, the frequency of signals varied. For this reason, the signals were preprocessed and discretized, in order to obtain the position of each object at regular timestamps: 1 minute for Ships and 15 seconds for Planes and Taxis.
After this normalization step, some location signals were observed to be incorrect, because they implied object movements at an extremely high speed (these kind of errors are not rare in GPS measurements). We set a maximum speed for each dataset, 800km/h for Planes and 234km/h for the other datasets. If the signal of an object at a timestamp implies exceeding the maximum speed, that signal is deleted. This amounts to deleting 0.01% of the signals from Ships and 0.83% from Planes. In addition, the frequency of signals is variable. For example, in Ships the frequency of signals from most boats was regular when they were sailing, but it stopped for long periods, typically when they were at a port. This produces disappearances and appearances, which amount to 0.21% of the entries in Ships and 0.96% in Planes. Finally, since our setup requires the position of each object at regular timestamps, we interpolated the two closest signals if they were less than 15 time instants apart, since an object cannot move far during a short period of time; otherwise we consider it a disappearance/appearance. Each dataset is stored in a plain text file comprising four columns: object identifier, time instant, coordinate x and coordinate y. Each value in the columns is stored as a  string. In order to obtain a better compression comparison, the data was alternatively stored in binary form, with the minimum number of bytes. For example, in Ships, two bytes were used to represent object identifiers (max value 3,654), two for the time instant column (max value 44,642), two for the x-axis (max value 12,719) and three for the y-axis (max value 368,186). Table 3 shows the resulting sizes.
To illustrate how compressible the data is, the plain data files were compressed with p7zip 7 with default settings. As we see, p7zip compresses the data to 9.5%-11.8% of its binary representation.

Compression
The first experiment compares the compression ratio of GraCT with scdcCT on the two real datasets, Ships and Planes. GraCT and scdcCT data structures were built on both datasets using snapshot intervals of d = 120, 240, 360 and 720 time instants. The construction time for the complete structure was approximately one minute. Table 4 shows the results of compression. The first row gives the size of both structures for each snapshot period d, whereas the next two rows separate the space due to the snapshots (i.e., the k 2 -trees) and to the compressed logs. It can be seen that the space of the snapshots is much lower than that of the logs, even with the highest sampling rate. Still, the compression of the logs themselves improves for larger d. In this aspect, Re-Pair improves faster than SCDC: the quotient of the log space with d = 720 versus d = 120 is 0.68 on Ships and 0.47 on Planes, whereas with scdcCT they are 0.78 and 0.60, respectively.
The last two rows give the compression ratios obtained with respect to the plain and binary representations, respectively. The compression achieved by scdcCT with d = 720, 13%-15% of the binary representations, shows that statistical compression of the differences is competitive with state-of-the art compression like the one offered by p7zip. The grammar-compression used by GraCT, on the other hand, obtains less than half the space of scdcCT, and less than 60% of the space used by p7zip.
We remark that, within this space, GraCT not only represents the trajectory data, but it offers efficient access and queries on it. Devoided of the extra information it stores on the nonterminals to speed up queries, the space used by the grammar compression  of the logs would be halved or so. This shows that grammar compression is much more efficient than statistical compression to exploit the redundancy of real-world trajectory data.

Query performance
This experiment studies the response times of GraCT and scdcCT, with different snapshot periods, on the real datasets. Figure 11 shows the average time (in milliseconds) per query for the following seven types of queries: • object: this query obtains the position of a given object at a given time instant t q .
We averaged over 20,000 queries for random objects and time instants.
• trajectory: this query returns the trajectory followed by an object during an interval [t b , t e ], where t e −t b is fixed at 2,000 time instants. We averaged over 10,000 queries for random objects and time instants t b .
• time-slice S: this query obtains the identifiers and positions of the objects lying within a small region (40 × 40 cells) at a given time instant. We averaged over 1,000 queries for random region positions and time instants.
• time-slice L: this query obtains the identifiers and positions of the objects lying within a larger region (320 × 320 cells). We averaged over 1,000 queries for random region positions and time instants.
• time-interval S: this query obtains the objects present in a small region (40 × 40 cells) at any time instant between t b and t e , where the interval size is t e − t b = 100 instants. We averaged over 1,000 queries with random region positions and instants t b .
• time-interval L: this query obtains the objects present in a larger region (320 × 320 cells) over a longer interval of size t e − t b = 500 instants. We averaged over 1,000 queries with random region positions and instants t b .
• knn: this query obtains the K nearest neighbors to a given position at a given time instant, where K is a random value between 1 and 50. We averaged over 1,000 queries with random objects and time instants.

Object
As seen in Figure 11(a), GraCT outperforms scdcCT in this query, answering within 2-4 microseconds. The logs in GraCT have far fewer entries than those in scdcCT, and in many cases the extra information stored on nonterminals allows processing of each GraCT log entry in constant time, just as a simple scdcCT log entry. Further, as previously demonstrated, Re-Pair compression improves faster as the logs increase, which makes GraCT times less sensitive than scdcCT to the increased distance between snapshots. The improvements of GraCT over scdcCT range from 5.9% with a snapshot period of d = 120 to 46.3% with d = 720 on Ships, and from 12.3% to 100.4% on Planes. Figure 12 shows the relative time spent on the k 2 -tree of the snapshot versus the time spent on the logs for the different queries. Since the times spent on the snapshots are identical, it follows that GraCT processes the logs way faster than scdcCT. This demonstrates that, provided with the proper enhanced data, the use of grammar compression also yields good query times, outperforming the statistical encoding of trajectories, because it processes fewer entries at about the same speed per entry.

Trajectory
In trajectory queries, as shown in Figure 11(b), GraCT also outperforms scdcCT, but the difference in response times is small: 12%-20% on Ships and 8%-13% on Planes. The majority of the time in this query is spent in decompressing the trajectory (of 2,000 entries), which is done at about the same rate in statistical (scdcCT) and grammar-based (GraCT) compression: 25-45 nanoseconds per decompressed datum on Ships and 15-20 on Planes. Since the trajectories to recover are significantly longer than the snapshot periods, snapshots mostly disrupt the processing of the logs, which explains why times improve with fewer snapshots. Figure 12 shows the case of shorter trajectories, of lengths 100 (Trajectory 100) and 1,000 (Trajectory 1000), on Ships. As expected, more time is spent on the logs to retrieve longer trajectories or with longer sampling periods.

Time-slice S and time-slice L
In time-slice queries, the difference in performance between GraCT and scdcCT becomes much more significant. As shown in Figures 11(c) and 11(d), GraCT performs 2-4 times faster than scdcCT in both datasets, running each query in less than a millisecond.
This query can be regarded as running an object query on each of the candidate objects, but as shown in Figure 12, the algorithm spends more time working on the log than for object queries, because the k 2 -tree amortizes its traversal cost when retrieving many points. Thus both indexes spend more time in the traversal of the log, where GraCT outperforms scdcCT, and less on the k 2 -tree, where they perform similarly. As a result, GraCT can be up to 4 times faster than scdcCT on time-slice queries, whereas in object queries it was up to twice as fast.

Time-interval S and time-interval L
As observed in Figures 11(e) and 11(f), time-interval queries are very similar to timeslice queries with respect to the comparison between GraCT and scdcCT, although this time GraCT is only up to 3 times faster than scdcCT. On longer time intervals, the difference is smaller because both structures may have to decompress the movements along the interval [t b , t e ]. The difference is larger than on a trajectory query, however, because GraCT can avoid decompressing nonterminals if their MBRs are disjoint or contained in the query range. Table 5 evaluates the improvement we obtain due to the use of MBRs of nonterminals on Ships. As it can be seen, the use of MBRs has a relatively small impact on time-slice queries (below 5%), but a very significant one on time-interval queries: increasing with the sampling period, the speedup factor goes from 5.6 to 12.7 on time-interval S and from 4.9 to 7.4 on time-interval L queries. Without using MBRs, indeed, GraCT would be slower than scdcCT on time-interval queries.

Knn
Nearest neighbor queries also require processing the logs of many candidate objects. However, the process done on each log entry is significantly heavier than on previous queries: objects are extracted from and reinserted in priority queues, their d min distances are recalculated, etc. In this scenario, the fewer number of entries processed by the grammar-compressed logs of GraCT, which can skip whole nonterminals, is even more   valuable. As a result, GraCT is up to an order of magnitude faster than scdcCT for these queries, answering them in less than 3 milliseconds. Figure 12 shows that the prioritized traversal of the k 2 -tree consumes a significant part of the query time on GraCT, but on scdcCT it is counterweighed by the more expensive processing of the log.

Scalability
In order to evaluate the scalability of our structure in compression and query times, we create five prefixes of Taxis with sizes 500 MB, 1,024 MB, 5,120 MB, 10,240 MB, and 61,240 MB. In each longer prefix, more objects are incorporated to the dataset. We built GraCT and scdcCT on each dataset, with the period of snapshots set to 720. Over the GraCT indexes, we ran queries with the same settings of Section 7.3.
As shown in Figure 13(a), the compression ratio of GraCT improves as the dataset grows, from 9% to 4%. This is a consequence of exploiting repetitiveness via grammar compression. As the number of objects and trajectories of the same kind increases, finding similar trajectories is more likely. Instead, as it can be seen, statistical compression is much less sensitive to repetitiveness, improving only from 12% to 10%. This shows another advantage of grammar compression when handling very large datasets.
With the growth of the data, on the other hand, the depth of the grammar used by GraCT also grows, logarithmically in case the grammar is balanced. This increases the cost of navigating inside a nonterminal, for example. This logarithmic growth is apparent in the object queries (Figure 13(b); since the x axis is logarithmic, a logarithm looks like a straight line), though in trajectory queries (Figure 13(c)) the time is amortized over the 2,000 output positions (those are decompressed in constant time each from the grammar, with only a logarithmic additive penalty for the whole query). The other queries grow linearly because, as more objects are added on a fixed grid, the output size and/or number of candidates to consider grows proportionally to the data size. This is apparent in Figure 13(d), where the query times divided by dataset size are essentially constant.

Comparison with a classical uncompressed index
As anticipated, we compare GraCT with the MVR-tree [40], a classic spatio-temporal structure based on the R-tree. The MVR-tree comprises multiple versions of the same data, each associated with a different interval of time, and an R-tree representing the positions of the objects during that interval. The intervals of the versions are disjoint. The MVR-tree exploits the repetition of nodes along consecutive R-trees, by sharing common nodes between versions. It is important to note that this spatio-temporal structure is efficient at solving time-slice, time-interval, and knn queries only, which are solved by identifying versions that overlap with the queried time interval. The algorithm traverses the R-tree of each such version, following the nodes involved in the spatial query range. Queries object and trajectory, instead, are costly for the MVR-tree. These are speeded up in the full structure, the MV3R-tree, by using an auxiliary 3DR-tree structure [40] (which requires even more space). We therefore omit these queries in our comparison.
The MVR-tree was configured to run in main memory. Both structures were built on the Ships and Planes datasets described in Section 7.1, and GraCT was configured with the same snapshot distances used in the previous experiments (120, 240, 360 and 720). The size of the MVR-tree in Ships was 12.16 GB and in Planes it was 11.72 GB, while the maximum-space configuration of GraCT used 36.40MB and 51.63MB, respectively, two orders of magnitude less (more precisely, GraCT uses 342 times less space on Ships and 232 times less space on Planes). The MVR-tree is also 30-40 times larger than the raw data represented in binary form.
We now study if this gigantic difference in space induces a significant impact on query time performance. Figure 14 shows the times of both structures averaged over 1,000 random queries of each type: time-slice S, time-slice L, time-interval S, time-interval L, and knn queries. Instead of the sampling d = 720, we include a new one, d = 60, where GraCT uses 50.70 MB on Ships and 83.14 MB on Planes. This is still 245 and 144 times smaller than the MVR-tree, respectively, and is sufficient for GraCT to be as fast as MVR-tree on most queries.
Time-slice queries are the most efficient ones for the MVR-tree. Since they retrieve objects within a given region at a single time instant only, the MVR-tree needs to traverse just one R-tree. When the traversal of the tree reaches the leaves, it can retrieve more than one object per leaf, which is faster than traversing one log of GraCT per candidate. Figures 14(a) and 14(b) show that, even in this case, the maximum-space configuration of GraCT matches the performance of the MVR-tree on Ships, and outperforms it on 38  Planes. It is likely that the higher locality of reference when accessing a structure that is two orders of magnitude smaller plays in favor of GraCT. This advantage of the MVR-tree vanishes on time-interval queries. In an interval of time, the MVR-tree may have to traverse a number of R-trees, which makes it much slower and always loses to GraCT, as seen in Figures 14(c) and 14(d). Figure 14(f) studies the turning point on Ships by varying the length of the queried time interval, using the smallest space configuration of GraCT (snapshot distance 720). As the time span of the queries increases, both structures slow down. The MVR-tree times, however, increase much faster, becoming slower than GraCT when the time interval exceeds 140 time instants on time-interval S and 60 on interval L. GraCT is thus less sensitive to the length of the queried interval, outperforming MVR-tree already in short time intervals even on its minimum-space configuration.
Nearest neighbor queries refer again to a single time instant, and therefore the MVRtree can handle them within a single R-tree, even if it usually has to traverse many more nodes than for a time-slice query. In this query, GraCT also has to traverse a large number of k 2 -tree nodes and verify many candidates because, unlike the MVR-tree, it does not have a snapshot of the precise query time instant. As it can be seen in Figure  14(e), the MVR-tree clearly outperforms GraCT on Ships, though not on Planes, where GraCT always improves and the MVR-tree always degrades with respect to Ships.

Disk and full precision
With the minimum-space configuration (snapshot sampling of 720), GraCT uses nearly 600 times less space than the MVR-tree. We can then use GraCT to index in main memory datasets for which the MVR-tree would by far exceed the available memory, and would then have to operate on disk. As we show next, even a time-slice query (where the MVR-tree most prominently excels) requires around 50 milliseconds to be solved by the MVR-tree on disk, which is 50 times slower than the time required by GraCT in main memory. If we use a snapshot sampling of 60, then GraCT becomes more than 300 times faster, and is still 150-250 times smaller than the MVR-tree.
The compression of GraCT is achieved, however, by means of discretizing the data and therefore losing precision. In most real systems, full precision is anyway unnecessary. For example, using eight decimals for the GPS position of an object discriminates micrometers, which is much finer than what most physical devices can acquire, and only useful in very specialized fields like tectonic plate mapping.
Still, if full precision is required, GraCT can be used as an in-memory cached index of the full-precision data stored on disk, exploiting its very small main memory footprint. To answer a query, GraCT is first used to obtain a list of candidate answers, which satisfy the query up to where the precision of GraCT can discriminate. Then, in a second step, those candidate answers are verified by reading the full-precision data from disk.
As a proof-of-concept of the feasibility of this approach, we designed a simple arrangement for time-slice queries, where the MVR-tree performs best. We use the GraCT data structures of the Ships dataset with 720 time instants between snapshots. In addition, we store the full-precision data in a relational table of PostgreSql DBMS v9.1.24. 8 The table stores tuples object, time, x, y to represent all the trajectories of all the objects. To efficiently verify the answers to time-slice queries, a hash index is set on the first two columns (together). For simplicity, we still assume that the time instants are discretized; a B-tree can be used for the second column in case they are real numbers as well.  Table 6: Response times of time-slice queries on disk, in milliseconds.
As a baseline, we use the MVR-tree data structure on disk. GraCT, including the fullresolution data and the hash indexes, occupies 4.18 GB, whereas the MVR-tree occupies 78.61 GB. That is, GraCT still requires 18 times less disk space than the MVR-tree.
To emulate a scenario where the full-resolution data does not fit in main memory, we configured our machine to use 2 GB for the kernel, so that neither the MVR-tree nor the PostgreSql table can be completely cached by the operating system or the DBMS.
The experiment begins in cold state, with all the data structures residing on disk and the caches of the operating system and the DBMS completely empty. The first row of Table 6 shows the average time of 1,000 executions of the first time-slice query in cold state. That query requires completely loading the GraCT data structures into main memory. In the case of the MVR-tree, the execution only loads the nodes of the R-tree traversed by the query. The second row shows the times resulting from running 1,000 time-slice queries after running the first one, that is, in warm state with the GraCT index already loaded in memory. We independently tested the second, third, and so on, seeing no significant differences.
Even in cold state, GraCT is 10% faster than the MVR-tree for time-slice S queries, but once in warm state, the response time of GraCT drops to 3.57 milliseconds, whereas the MVR-tree requires around 45 milliseconds, which is about 12 times slower. GraCT is much slower than the MVR-tree for time-slice L queries in cold state, on the other hand, but again, once we enter in warm state, GraCT becomes about 8 times faster. That is, a disk-based system on regime, with full precision using GraCT as an in-memory cached index is an order of magnitude faster than the MVR-tree and uses much less disk space.
Other queries, like time-interval, require a more sophisticated arrangement to use GraCT as an in-memory cache. For example, we can store the consecutive points of each trajectory in contiguous form, indexed by time with a B-tree, and a hash-indexed table of objects pointing to their corresponding trajectories. This enables an efficient traversal of the trajectories of candidate objects in the time intervals where GraCT estimates that they may intersect the query range. Those arrangements also support the simpler queries object and trajectory. Our experiment suggests that it is much more efficient to verify potential answers on disk than to find them using a disk-based structure.

Conclusions
This article introduces GraCT, a compact data structure for the representation of moving objects. In our experiments on real-world trajectory datasets, GraCT compresses the data to 4%-7% of its raw representation size. Within this space, GraCT supports not only direct access to the data, but also various queries looking for objects in spatiotemporal windows. The trajectory data is manipulated directly in compressed form, with only small portions requiring decompression in order to answer queries. The space also includes the storage of spatial indexes representing snapshots at regular time intervals, which speeds up queries. The main algorithmic novelty of GraCT is the representation of the logs of movements between snapshots using grammar-based compression. On our real-life datasets, grammar-based compression captures the repetitiveness arising from recurrent trajectories, outperforming statistical compression by a factor of 2 and improving as the datasets become larger. A second novelty is to enrich the grammar by associating information with the nonterminals, such as their relative MBR. This enables GraCT to traverse the logs very quickly without decompressing nonterminals in most cases, thereby answering queries within a few milliseconds. A comparison with the MVRtree, a classical index, shows that GraCT is two orders of magnitude smaller and still competitive in query times when both run in main memory. When the MVR-tree cannot fit in main memory, GraCT is one or two orders of magnitude faster. This includes an arrangement where GraCT is used as a low-precision main-memory cached index that must verify its result on disk, where the full-precision data is kept.
Future work in this area will explore other compression mechanisms that exploit repetitiveness, looking for ways to process the trajectory data without decompressing it. Some alternatives are the use of block trees [2] or relative Lempel-Ziv [18]. We also plan to further use GraCT as a main-memory cached index, extending our arrangement to handle all the queries considered in this article.