Visualization of the Strain-Rate State of a Data Cloud: Analysis of the Temporal Change of an Urban Multivariate Description

One challenging problem is the representation of three-dimensional datasets that vary with time. These datasets can be thought of as a cloud of points that gradually deforms. However, point-wise variations lack information about the overall deformation pattern, and, more importantly, about the extreme deformation locations inside the cloud. This present article applies a technique in computational mechanics to derive the strain-rate state of a time-dependent and three-dimensional data distribution, by which one can characterize its main trends of shift. Indeed, the tensorial analysis methodology is able to determine the global deformation rates in the entire dataset. With the use of this technique, one can characterize the significant fluctuations in a reduced multivariate description of an urban system and identify the possible causes of those changes: calculating the strain-rate state of a PCA-based multivariate description of an urban system, we are able to describe the clustering and divergence patterns between the districts of a city and to characterize the temporal rate in which those variations happen.


Introduction
One challenging problem in data analysis is the representation of three-dimensional discrete data [1]. This analysis becomes harder when a time-dependent change of the data takes place, introducing a new temporal dimension. In the present article, we explore formal approaches to quantify the temporal change of discrete three-dimensional data. Specifically, we build a methodology to assess the transformation of a data cloud that is derived from a Principal Component Analysis (PCA): a 13-year multivariate description in [2] that provides a reduced description of an urban system given only by the first three principal components. Since the points represent an abstraction of an urban system, one main goal is to understand the temporal variation of the multivariate description of the districts in order to analyze the behavior of the overall city in the time-span. Our main hypothesis is that these three-dimensional datasets can be thought of as a cloud of points that gradually deforms.
However, the challenging issue is that deformation between consecutive times cannot be visualized accurately. There are some methods to overcome this difficulty. One is the vector plot of three-dimensional displacements or velocities, which is typically used to visualize results in computational mechanics applications [3][4][5]. For example, these methods are used in the kinematic visualization of motion [6][7][8][9][10], but they are restricted to display the relative motion among the data, and are not able to identify the most dynamic regions of the dataset. Another is the parallel coordinate technique that successfully exhibits the temporal change of highly-dimensional statistical and information datasets [11]. However, multi-dimensional data are mostly segmented into two-dimensional subsets, which are easier to deal with, Similar to the computer tomographic scans of medical imaging [1, 12,13]. However, none of these methods are suitable for understanding the patterns of diversification or conformation, which are closely related to the temporal change of the differences between join data values: the maximum and minimum magnitudes of variation and the evaluation of their direction can be significantly helpful to identify differentiation patterns in the data [14]. Conversely, they can be beneficial for locating uniformity in a dataset that was previously differentiated.
One common approach to understanding the temporal change of the dataset is to use confidence ellipses in dimensionally-reduced representations (i.e., PCA) [15][16][17]. This type of approach has been applied to several time-dependent problems (e.g., [18][19][20][21]). Confidence ellipses quantify the dispersion of the results: the main directions and intervals of variance accounting for the spatiotemporal distribution are calculated first, and then, the ellipses (or ellipsoids in 3D) are drawn using those intervals as the main axes. A temporal analysis can be achieved by evaluating the geometrical characteristics of the ellipses, especially their area, to which the change of the distribution can be associated. Their main drawback lies in the failure to locate the most dynamic areas in the dataset [22].
The field of continuum mechanics provides a measure of the temporal variation of the distance in between points: the strain-rate tensor (see, for instance, [23]). The continuum mechanics theory-which arises from the classical Newtonian mechanics-analyzes the causes and effects of motion for a deformable media composed by an infinite group of particles. When a continuous media is being deformed in various directions at different rates, the strain-rate of a certain position in the medium cannot be expressed by a scalar value solely. It cannot even be expressed by using a single vector. Instead, the rate of deformation must be expressed by the rank-two strain-rate tensor with its components determined by the positional derivatives along each spatial dimension. Hence, the mathematical framework of tensors can determine exactly the deformation that is accumulated in a location inside the medium-which is typically subject to the imposition of displacements or loads. This tensor is commonly used to detail the amount of elastic energy in the physical descriptions of multiple materials, such as solids or fluids (see [24] for a complete mathematical exposition). Most of those models are formulated as the product of a constitutive tensor and the strain-rate tensor, giving the stress condition of the material that is balanced in the kinetic equations. In the present study, the calculation of the strain-rate tensor is not related to the kinetics of any material, and, thus, it can only be a mathematical tool that supports the examination of the deformation rates in the discrete statistical data.
However, the strain-rate tensor arises from the continuum assumption, and discrete displacements of points rather than continuous distributions take place in the variation of the data cloud. Typically, the issue of applying derivatives to discrete displacements of points is solved by several different approaches. Few statistical techniques use co-variance functions to represent directly the strain-rate field (see, e.g., [25]). The common approach is to compute a continuous version of the displacement-or velocity-field, so that the derivatives can be applied to the continuous displacements. Some methods in this line have interpolated the discrete displacements by minimizing the residual-or distance-between the continuous interpolation and the discrete version [26]. Other interpolation techniques weight the distance between an interpolated piece-wise continuous field and the discrete displacement field (see, e.g., [27]). This method results in a minimization technique where a continuous strain-rate field can be derived. For example, the piece-wise continuous field can be defined as being splines or the widely used linear polynomials in variational formulations [28]. These techniques have been applied in Earth science and medical imaging works [29][30][31], as well as in the strain-rate calculation of geodetic observations [32,33].
Another fundamental issue is the representation of the strain-rate state. One of the possible techniques that can help to visualize the deformation rate of the dataset is to plot the main components of the tensor using strain-rate diagrams, where concentrations of strain-rate patterns can be displayed as vector fields (see, for example, the ones in geodetical observations of the Earth's mantle [34][35][36]). The main drawback of strain-rate diagrams is that the strain-rate components are visualized as the projection of three-dimensional vector fields into the two-dimensional framework, and, therefore, the third-dimension component has to be necessarily neglected. Another method, more suitable to two-dimensional plots, is the contour graph of principal stresses, where the stress patterns in structural elements [37,38] and tectonics [39] are visualized using continuous lines that depend on the stress magnitude. That method overcomes the three-dimensional issue, but it does not give insights about the orientation of the principal stresses. A dual form of the contour graph is to calculate the family of curves that are instantaneously tangent to the extension and contraction components of the strain-rate tensor: the so-called trajectory curves in the continuum mechanics field [23]. The trajectories of each principal component of the stress can be depicted in a separated plot with the stress magnitude colored along the trajectory line such that stress patterns are completely shown in a two-dimensional framework.
Since a robust methodology that describes the temporal change of the urban system-represented by a multivariate dataset-has not been carried out before, we perform a quantitative analysis by including the strain-rate tensor as the fundamental metric. In this work, we calculate the strain-rate state of the discrete dataset without a priori assuming the mechanisms by which the system experiences transformation. To apply the continuum mechanics principles into the discrete dataset, we use interpolation methods, such as the ones applied in discrete variational formulations (i.e., Finite Element Methods (FEM), Particle Methods, Collocation Methods, Mesh-less methods, etc.). Specifically, we derive the three-dimensional strain-rate tensor from a FEM interpolation of the discrete velocity field, as introduced in [40][41][42]. We include the trajectory lines of the principal strain-rate components as the methodology that can be used in a computational (two-dimensional) framework for visualizing the main patterns of change in any time-dependent data cloud. This approach overcomes the three-dimensional representation problem by separating the strain-rate state of the cloud in several two-dimensional plots-one for each principal component of the strain-rate-which exhibit the magnitudes and orientations of the deformation patterns.
The remaining parts of this document are organized as follows. In Section 2, the methodology to compute the discrete version of the strain-rate tensor is presented. Since the main problem involves the calculation of the derivatives of discontinuous-discrete-velocities and the visualization of the strain-rate patterns, we extensively detail the numerical techniques that are adopted in the present approach. Next, in Section 3, we present the application of the methodology to the case study-the urban system of Barcelona-by visualizing its main strain-rate patterns, meaning the city's environmental, social and economic change. In Section 4, some conclusions of the proposed methodology close this article.

Methods
We begin this section with a review of the strain-rate tensor calculation provided a discrete three-dimensional data cloud. For doing so, the formal problem of the time-dependent dataset is introduced first. Then, we explain the numerical techniques that transform the discrete dataset into a mathematical framework by which the strain-rate tensor can be computed. Most of the ideas rely on the geometrical analysis of the discrete dataset by computing the spatial discretization of the dataset into geometric elements through a Delaunay Triangulation. After doing that, the computation of the strain-rate is performed with a FEM interpolation of the velocity field. Finally, we address the eigen-problem for the strain-rate tensor, such that the solution of the eigenvalues/eigenvectors gives the extrema strain-rates at each geometric element. The flow chart diagram of this methodology is abstracted in Figure 1, including the main outcomes of each step. The extended explanation of the methodology is developed along the present section.

Time-Dependent Three-Dimensional Dataset
Since the main objective of this work is to reveal the temporal transformation of a three-dimensional and time-dependent dataset, let us first introduce some notation in order to clarify the mathematical ideas to be used. Let us define the discrete time-dependent data to be the set of points P = {p i }, with i = 1, 2, . . . , m, being m the total data. The values in each one of the three dimensions can be seen as scalar coefficients for a set of basis vectors. This tuple of components compose the vector that we call the position or coordinate x i = [x i,1 x i,2 x i,3 ] , with the superscript denoting the transpose operation, the first subscript referring to the point i and the second to the dimension. Hence, we call X n (t) = {x 1 (t) , x 2 (t) , . . . , x i (t) , . . . , x n (t)} ∈ R 3 the positions of the points in a certain time t. Let us consider a uniform partition of the time span t ∈ [t d , t g ] in a sequence of discrete time-steps t d = t 0 < t 1 < . . . < t n < . . . < t N = t g , with δt > 0 the time-step-size defining t n+1 = t n + δt for n = 0, 1, 2, ..., N. Thereby, we use the superscripts to denote the discrete time-steps, with the only exception of denoting the transpose operation with the superscript .
Since the time-dependent dataset of the case study comes from a PCA reduction of a higher-dimensional multivariate dataset Y n (t) = {y 1 (t) , y 2 (t) , . . . , y i (t) , . . . , y n (t)} ∈ R k , with k >> 3 and t ∈ [t d , t g ], into a lower-dimensional one X n (t), t ∈ [t d , t g ] that possesses only three independent dimensions, namely Principal Component 1 (PC1), Principal Component 2 (PC2), and Principal Component 3 (PC3), we use the Cartesian coordinate system straightforwardly with each principal component being a dimension. That is, Hence, the discrete time-dependent dataset can be thought of as a cloud of points in the three-dimensional space that deforms gradually over time.

Finite Element Method Interpolation
The main idea of the present approach is to transform the discrete cloud of points into a mathematical framework-similar to a deformable medium-by which the strain-rate tensor can be computed. To do so, we generate a mesh T h (t) = {K} from the set of points P that is composed of non-overlapping and conforming geometrical elements K of diameter h. There are several methods to generate a mesh from a set of points, all which are studied in the computational geometry field [43]. Here, we apply Delaunay triangulation DT (P ) for several reasons: (1) the aspect ratio of the triangulated elements produce a high-quality mesh; and (2) fast Delaunay triangulation algorithms have been developed recently (see, for example, the one in [44]).
The result of applying the Delaunay triangulation over the set of points is a discrete mesh T h := DT (P ) that possess the following characteristics: it covers exactly the convex hull Ω of the point set, no point p i is isolated from the triangulation, and all the elements {K} are four-point tetrahedrons, which are completely defined by the position of their four corner points K := x j , with j = 1, 2, 3, 4. The generated mesh T h = DT (P ) can be seen as a-material-domain Ω that suffers deformations from the displacements of the points between consecutive time-steps. Since only discrete displacements between consecutive time-steps are known for the set of points, we now explain how the continuous velocity field inside the mesh is calculated. Even though the FEM has been used to perform interpolation using the point-wise data (see [41,42] for applications in one and two-dimensional data), in this work, we apply this well-known method in a three-dimensional setting. In FEM, the finite interpolating space V h is defined as made of continuous piece-wise polynomials N(x) in the mesh T h , where the discrete approximation F h (x, t) ∈ V h of any multi-dimensional function F(x, t), x ∈ Ω, can be written as We use the simplest finite element: the tetrahedron with linear polynomials and four nodes. Again, some notation is required to define the polynomials inside the element. The set of normalized coordinates η 1 , η 2 , η 3 , η 4 in each tetrahedron K is such that the value of η i is one at the point p i ∈ K, zero at the other three corner points, and varies linearly from that point to the opposite edges. This set of coordinates has the property that the sum of the four coordinates (each belonging to one tetrahedron point) in any location inside the tetrahedron is identically one: Hence, the shape functions inside each linear tetrahedron are defined to be these coordinates: N i (x i ) = η i (x i ), with i = 1, 2, 3, 4 denoting the corner points. The FEM interpolation in Equation (1) of a three-dimensional vector function, e.g. F (x, t), can be defined inside each linear tetrahedron K as The way the tetrahedral coordinates η i , i = 1, 2, 3, 4, are defined is by means of the previous interpolating relation together with the summation constraint. That is, when one aims to define the tetrahedron geometry and calculate any position inside the tetrahedron in order to obtain the tetrahedral coordinates system where the coefficients of the inverted matrix are given by Here, the abbreviation x ij = x i − x j is used, and the volume υ can be calculated with the expression At this point, it is possible to calculate the spatial derivatives of any interpolated function ∂ ∂x F (x, t) in terms of the tetrahedral coordinates as The way to calculate the continuous stress-rate tensor field is through the derivation of a continuous version of the velocities. Hence, we calculate the continuous velocity field by means of the FEM, in which linear piece-wise polynomials are used to interpolate the velocity at any spatial position inside the mesh. First, let us explain how to calculate the discrete version of the velocity from the cloud of points. We suppose that the displacement s i of point p i in the time interval t n , t n+1 can be defined-without loss of accuracy-as infinitesimal, in the sense of With this supposition in hand, we use the Taylor expansion: in order to calculate the discrete velocity where the second (and higher) order terms are neglected. The previous result then generates the continuous version of velocity by replacing it in Equation (2).

Elemental Strain-Rate Calculation
Having defined the continuous space of velocities, we can calculate the derivatives along each of the spatial directions and derive the strain-rate tensor field.
Following the continuum mechanics concepts in [23] and assuming small deformations, the strain-rate tensor is calculated as with ∇v the gradient of velocity. Each component of the 3 × 3-tensor is developed in Cartesian coordinates as Moreover, the six independent components of the strain-rate tensor can be arranged using Voigt's notation into a six-component strain-rate vector as follows: where are the shear-rate strains.
With this notation in hand, the strain-rate tensor can be calculated as by defining the matrix operator of derivatives over the velocity field. In the case of the right hand side velocities, we can arrange a node-wise vector of discrete velocities in the tetrahedron K, as Using the definition of the finite element interpolation of any function in Equation (2) together with its partial derivatives in Equation (4), and replacing those in Equation (8), we obtain Now, the operation Hence, E (K, t n ) can be calculated as the product of the matrix S (K) and the vector V (K, t n ). That is, with x ∈ K and the discrete matrix S (K, t n ) defined as Thus, this last matrix can be computed solely in terms of the coordinates of the nodes.

Principal Strain-Rates
Up to this point, we have demonstrated how to calculate the elemental strain-rate. Now, our purpose is to identify the data cloud transformation throughout the visualization of the strain-rate patterns. That is, we need to identify the extrema strain-rates and their orientations. In a formal sense, this is the well-known eigenvalues and eigenvectors problem, which is stated as: if T is a linear transformation from a vector space V over a field F into itself, and v is a vector in V that is not the zero vector, then v is an eigenvector of T if T(v) is a scalar multiple of v. Knowing that by definition the second-order strain-rate tensor is a linear operator from a vector field into another first-order tensor field, the previous definition applied to the strain-rate tensor leads to: where I is the 3 × 3 identity tensor; n (K, t n ) ∈ R 3 is a normalized (non zero), i.e., unit, vector called eigenvector; and λ (K, t n ) ∈ R is the eigenvalue associated with the eigenvector. In other words, an eigenvector is a vector that changes by only a scalar factor when the strain-rate tensor is applied to it, resulting in a vector parallel to itself. By solving Equation (12) ,one obtains three different eigenvalues λ 1 (K, t n ) , λ 2 (K, t n ) , λ 3 (K, t n ) , and three eigenvectors n 1 (K, t n ) , n 2 (K, t n ) , n 3 (K, t n ) , associated with each eigenvalue.
The eigenvalues and eigenvectors describe the principal magnitudes and orientations of the strain-rate tensor: since the diagonal components of the strain-rate tensor E 11 (K, t n ) , E 22 (K, t n ) , and E 33 (K, t n ) have different values in different reference systems, one finds with the set of eigenvalues the extreme-maximum and minimum-possible values that any of these components may take. Indeed, the maximum and minimum stress-rates-and their orientations-are related to the maximum and minimum eigenvalues. In this work, we follow the notation in which positive values for the eigenvalues represent the extension-rate and negative values represent contraction-rate. Hence, λ 1 (K, t n ) is the maximum and positive eigenvalue meaning extension-rate, λ 3 (K, t n ) is the minimum and negative eigenvalue meaning contraction-rate, and λ 2 (K, t n ) is either extension or contraction rate, but in smaller magnitude.
Hence, with the extrema strain-rates at the elemental level, we can reveal the deformation trend of the data cloud, and, above all, locating which regions suffer the most abrupt change in the time-span. We also propose to draw the family of curves-trajectories-that are instantaneously tangent to λ 1 n 1 (K, t n ) , λ 2 n 2 (K, t n ) , and λ 3 n 3 (K, t n ) in the complete mesh Ω. This results in a graph plot for each of the λn (K, t n ) fields that can be easily represented in a two-dimensional framework. These tangent lines can be colored according to the eigenvalue magnitude such that they completely illustrate the main patterns of change inside the data cloud. Note, nevertheless, that λn (K, t n ) is a composition of a vector using tensor components. Those differ in formal definition, but we use this concept merely for visualization purposes.
Nevertheless, we also perform some temporal statistics of the strain-rate states, where the principal strain-rates for each time-step λ i n i (K, t n ) are accounted as the temporal events: each strain-rate state is accounted as a single observation. Specifically, we calculate the time-average of the eigenspace components λ i (K) and n i (K), the maximum extension-rate in the time span L ∞ (λ 1 (K)), and the maximum contraction-rate in the time span L −∞ (λ 3 (K)). In this way, we calculate representative deformation results for the complete temporal interval of analysis, as well as the location of extreme values. Besides, trajectory lines can be plotted using these statistical results.

Results
In the present section, we demonstrate the application of this methodology to quantify the temporal change of an urban multivariate system First, we cite the case study that includes the multivariate description of the ten districts of Barcelona, and whose reduced three-dimensional dataset is used as the starting point. Then, we derive the strain-rate state of the dataset, pursuing the extension and contraction patterns visualization. Finally, we close this section with insights about the city transformation implied in the strain-rate state of the data cloud.

Time-Dependent Data Cloud from an Urban Multivariate Description
The time-dependent data cloud comes from the PCA output of a multivariate description of the city of Barcelona in [2]. Since 1987, the city has been divided into 10 administrative districts, which are the largest territorial units of the city and can be compared with neighborhoods in a common metropolitan area: Ciutat Vella, Eixample, Gràcia, Les Corts, Sarria, Sant Andreu, Sant Marti, Horta, Sants, and Nou Barris. Barcelona has a population of approximately 1.6 million inhabitants living in 10,216 ha. The inclusion of all 10 districts in the multivariate description was aimed at representing the city at its overall scale and allowing comparisons between them.
The raw multivariate description-from which the PCA is calculated-comprises the data of 40 environmental, economic, and social indicators for the ten districts in the time span of t 0 = 2003 ≤ t n ≤ 2015 = t N , n = 0, 1, ..., 12. Hence, the case study data cloud comes from a PCA reduction of the higher-dimensional multivariate dataset Y n (t n ) ∈ R 40 , into a lower-dimensional one X n (t n ) ∈ R 3 that possesses only three independent dimensions: PC1, PC2, and PC3. The dimensionally-reduced dataset from the application of the PCA is presented in Appendix A. Hence, the three-dimensional and time-dependent data cloud is composed of the coordinates X n (t n ) of the ten p i points defined in the sequence of N = 12 time-steps from 2003 to 2015, with the time-step size of δt = 1 year. These points are displayed in Figure 2, where all observations-districts each year-in the time-span are included. As the first step of our methodology, we applied Delaunay Triangulation (DT) to the data cloud. Specifically, we calculated the DT to the set of coordinates at each time-step X n (t n ). This resulted in a mesh T h (t n ) composed of nel = |K| non-overlapping tetrahedrons. Table 1 expands the resulting triangulation for year 2003, with the vertices information for the nel = 17 tetrahedrons. This triangulation is also depicted in Figure 3, where the following notation for the districts in the vertices is used: C is Ciutat Vella, E is Eixample, G is Gracia, H is Horta, L is Les Corts, N is Nou Barris, A is Sant Andreu, M is Sant Martí, S is Sants, and R is Sarria. Since the position x i (t n ) of a given point p i at a later time-step can surpass the initial tetrahedron's circumscribed sphere, we recalculated the mesh triangulation at each time step t n , n = 1, ..., 11. Table 1. Tetrahedral elements derived from the Delaunay Triangulation of the set of points.

Principal Strain-Rates
We computed the strain-rate tensor of each tetrahedron with the interpolated version of the velocities for the case study, such that linear piece-wise polynomial functions defined inside each tetrahedron were used in the FEM interpolation. Certainly, we supposed that the velocities came from an infinitesimal analysis, in which the higher-order terms of the displacement were neglected. The gradients inside each tetrahedron were also considered to be constant since the polynomial functions were of first order. Applying Equation (10), we computed the strain-rate tensor of every tetrahedron, E(K, t n ) for time-steps n = 0, . . . , 11. Note that displacements could not be calculated for the last year t 12 = 2015, and that the strain-rate tensor units are year −1 (for the case study).
We were interested in the magnitude and orientations of the principal strain-rates-extension and contraction-at the elemental level. Hence, the next step was to solve Equation (12) and obtain the eigenspace components (eigenvalues and eigenvectors) of the strain-rate tensor. For the sake of conciseness, we list in Table 2 the results of the principal strain-rates for the year 2003 solely. The application of this methodology to the case study is displayed graphically in Figure 4, beginning with the map of the ten districts of Barcelona as the abstraction of the multivariate and time-dependent dataset. The three-dimensional coordinates arising from the PCA output are displayed next. We also present next the triangulated mesh at the initial year 2003. It is clear from the visual inspection that the quantitative analysis of the temporal transformation is greatly justified, thus we calculated the strain-rate tensor over the FEM interpolation of discrete velocities and computed its principal components.

Trajectory Patterns of the Principal Strain-Rates
In favor of the analysis, we display the principal strain-rate components in a graphical way. One approach is to illustrate the patterns of extension-rate and contraction-rate using a vector representation, to what is referred as the strain-rate diagrams [29]. In that approach, the centroid of the tetrahedron serves as the location from which the principal components of the strain-rate tensor give a representative result inside the element. We draw the strain-rate diagram of the year 2003 in the sixth step of Figure 4, where extension-rate is represented by symmetric blue vectors λ 1 n 1 pointing out the centroid, and contraction-rate is represented by the red vectors λ 3 n 3 pointing in. However, it is hard to visualize the distribution of the principal strain-rates and their three-dimensional orientations using this type of illustration.    Our approach to ease the visualization and understanding of the strain-rate state is to draw the trajectories of the principal strain-rate components, as used for displaying stresses in beams and columns in [45]. In the following, we demonstrate our findings of the strain-rate state at the year 2003 using the trajectories visualization. In Figure 5, we display the principal components trajectories, where the lines are colored by the magnitude of the principal strain-rate and those are parallel to its orientation. From these representations, we can understand the magnitude and orientation of each principal component of the strain-rate tensor. More importantly, the trajectory patterns overlapped with the coordinates of the districts (in Figure 3) provide information about local regions of extension and contraction rates inside the urban description, where extension-rate patterns indicate differentiation and contraction-rate patterns indicate clustering-or homogenization. This can be observed from the color intensity: locations of high strain-rate are colored with great intensity, while the locations of low deformations are not visible since the white color represents null strain-rates. 1 Figure 5. The trajectory patterns of the principal strain-rates at the year 2003: first principal strain-rate (top); second principal strain-rate (middle); and third principal strain-rate (bottom). The color scale in these graphs determines the strain type: the red color is used to represent extension, while the blue color depicts contraction.
In the case of the first strain-rate component, which is shown in Figure 5 (top), we observe that the larger magnitude of extension-rate is localized in the middle of Nou Barris Sant Andreu, Sant Marti, Horta and Sants, and that it decreases near Eixample, Les Corts, and Ciutat Vella. Therefore, the main transformation is located at the group of nearby districts: Nou Barris, Sant Andreu, Sant Marti, and Horta. The extension-rate patterns are oriented from this group apart to Ciutat Vella, suggesting that there is a divergence of Ciutat Vella from the grouped districts. Indeed, the main extension pattern is oriented along the PC3 dimension and covers the grouped districts. The pattern which comprises the districts of Nou Barris, Sant Marti, and Sants and ends at Gracia and Eixample is of lesser importance.
Contraction-rate, on the other hand, is expressed by the third principal strain-rate component, which by definition is orthogonal to the first and second principal strain-rates. The third principal strain-rate component is shown in Figure 5 (bottom), where we can appreciate this orthogonality by noticing that the trajectories of the third principal strain-rate are perpendicular to the extension-rate pattern. We observe that the contraction-rate trajectories are mostly homogeneous, with a minor importance among Sant Marti, Sant Andreu, Nou Barris and Horta districts, and completely declining at Sants and Gracia. This direct relation between extension and contraction is found in solids deformations, where it is ruled by the conservation of mass-or Poisson ratio [23].
Apart from the extension and contraction patterns of the mesh, locations of smaller strain-rates are represented by the second principal component. Considering Figure 5 (middle), we recognize that the orientation of this strain-rate component is concentrated in the middle of Les Corts, Sarria, Horta and Nou Barris, and that it is directed towards Eixample, fading at Ciutat Vella. This principal strain-rate component is certainly orthogonal to the first and second components, but it implies a strain-rate pattern that is two orders of magnitude smaller.
In the previous lines, we have demonstrated the application of the trajectories diagrams of the principal strain-rate components as a powerful visualization technique of the three-dimensional strain-rate state of a data cloud. The strain-rate patterns can be used to analyze the system's development, for example, with the identification of regions with a special behavior: although there are some grouped districts in the case study, they are all separated at a high rate in the first and third principal components (spatial dimensions). Hence, those are differentiating themselves in the PC1 and PC3 description. On the contrary, low strain-rates can be an indication of stagnation, and, thus, an expression of inactivity where an abrupt change is unlikely to occur. That is especially the case of the Ciutat Vella district, which is separated from the grouped districts but is neither diverging nor converging to them.
One final remark is that our visualization approach is mesh independent. This means that the principal strain-rate trajectory plots coincide regarding the triangulated mesh: different triangulations will produce different positions, magnitudes, and orientations of the principal strain-rate components; nevertheless, the trajectory lines that join them coincide for any type of strain distribution.

Temporal Statistics of the Principal Strain-Rates
Plots of the principal strain-rate components trajectories can be completed for the remaining years of the time span, t n , n = 1, 2, ..., 11. Readings of the strain-rate streamline patterns for those years can be completed straightforwardly, as discussed in the paragraphs above. Nevertheless, we performed some temporal statistics of the strain-rate states, where the principal strain-rates calculated for each time-step were considered as the temporal events: each strain-rate state was considered as a single observation.
The first statistics that we performed is the time-average of the principal strain-rate components, separated as the first, second, and third principal strain-rates. Table 3 presents the time-averaged results of the principal strain-rates. In addition, in the last schematic of Figure 4, we plot the strain-rate diagram of the time-averaged principal strain-rates at the time-averaged centroid of tetrahedron elements. This figure gives insights about the orientation and magnitude of the first and third principal strain-rates, again by plotting the symmetric arrows pointing out for extension and pointing in for contraction. However, we complete our analysis by drawing the trajectory curves of those time-averaged strain-rate diagrams in Figure 6.  1 Figure 6. The trajectory patterns of the time-averaged principal strain-rates: averaged first principal strain-rate (top); averaged second principal strain-rate (middle); and averaged third principal strain-rate (bottom). In these graphs, the red color is used to represent extension, while the blue color depicts contraction.
With the aid of Figures 3 and 6, we analyzed the trajectory patterns of the time-averaged strain-rate state of the data cloud. In the case of the first and third strain-rate components, those are comparable to the ones described in the previous paragraphs for the year 2003. We observe a high extension-rate behavior among Nou Barris, Sant Andreu and Horta in the PC1 dimension. The contraction-rate, on the other side, is oriented in the PC2 dimension and it is mostly located in between Gracia and Eixample and decreases near Sants. In the case of the contraction-rate in the PC3 dimension, it is mostly present between Nou Barris and Horta, and in smaller magnitude between Horta and Sant Andreu. The contraction-rate between the remaining districts is negligible in all dimensions. Both the extension and contraction rates demonstrate trajectory patterns that are directed from the grouped districts towards the separated district of Ciutat Vella. However, the magnitude of the time-averaged strain-rates is much smaller than the ones obtained for year 2003. In the case of the time-averaged second strain-rate, its magnitude is greater for Eixample, Les Corts, and Sarria nodes than for the grouped districts. The orientation of the trajectories involving this second strain-rate component is parallel to the one linking Eixample and Les Corts to the grouped districts.
The second temporal statistics we calculated is the L ∞ -norm of the temporal strain-rate distribution. That is, we calculated the year the maximum extension and contraction rates occurred within each tetrahedron. The results of the application of the L ∞ -norm to the case study are presented in Table 4. We observe that the maximum strain-rates occurred in year 2003: either expansion or contraction. In addition, the important contraction-rate magnitudes took place between the years 2003 and 2010. This is not the case of the extension-rate magnitudes, which are more prevalent after year 2008.

Conclusions
In this study, we quantified the temporal change of a time-dependent and three-dimensional dataset. Contrary to other approaches [1, 7,11,14], we calculated the three-dimensional strain-rate state of the dataset based on the interpolation of discrete point-wise displacements-or data variations. We applied a technique in continuum mechanics using a FEM interpolation of non-overlapping linear tetrahedral elements that covers the three-dimensional dataset.
The methodology was demonstrated to exhibit regions of major deformation-rate. Departing from the calculation of the numerical strain-rate values, we introduced some data-visualization techniques that help to locate the magnitudes and orientations of the strain-rate state in the three-dimensional framework. This was the case of the principal strain-rate trajectories, which where demonstrated to be more detailed than other possible visualization techniques, e.g., strain-rate diagrams [32,36] or confidence ellipsoids [18][19][20][21][22]. The main difference with strain-rate diagrams is the ability of the present approach to display a continuum version of the strain-rate inside the dataset, and to separate the analysis into each principal component of the strain-rate state, while being mesh independent. Compared to the method of confidence ellipsoids, the kinematic basis of the proposed methodology makes an accurate description of the deformation patterns in the complete three-dimensional domain, rather than in a compact and partial region that contains the data distribution. It also allows a detailed temporal description of the change between one instant and another. This is specifically not possible using the method of the ellipsoids since that method accounts for the complete temporal distribution in order to evaluate change.
The calculation of the strain-rate state shows that the methodology is suitable for quantifying the temporal change of a reduced three-dimensional dataset describing the social, economic and environmental state of the city of Barcelona. In particular, high strain-rates are associated with the localized deformation of regions that represent the districts of the city in the time span of 13 years. The method not only portrays the similarities and differences-distances-between the districts of the city, similar to cluster analysis [46][47][48] or distance measures [49,50], which explore similarities and groups in data. The distinctive attribute of the present work is the feasibility to quantify the districts' differentiation with time: the strain-rate tensor provides quantitative information about local regions of extension and contraction, where extension-rate patterns indicate differentiation and contraction-rate patterns indicate clustering-or homogenization. In this regard, the strain-rate methodology can be potentially used as an extension of clustering analysis, since it reveals contraction and extension patterns of clusters. Conclusions about the divergence or clustering of districts in time can therefore be stated. For example, it reveals the time, location, and orientation of pressures affecting the inhabitants of certain districts of the city, essentially those which are rapidly diverging from the rest (e.g., the case of Ciutat Vella for the case study). This methodology locates regions where detailed action is necessary, as well as the foretelling of possible ruptures in the system.
Finally, we mention that this method is not limited to the study of the data change, but can also be applied to other descriptions: a natural sequel to the present article is the study of the time-dependent data cloud as a deforming elastic solid under equilibrium. Solving the inverse problem, namely the identification of the constitutive moduli of the deforming material, emerges as the first step for predicting the system's future state given its historical data.
Author Contributions: The authors equally contributed to the elaboration of this article.
Funding: This research received no external funding.
Acknowledgments: Lorena Salazar-Llano greatly acknowledges the financial support received by the pre-doctoral grant "Beca Rodolfo Llinás para la promoción de la formación avanzada y el espíritu científico de Bogotá" from Fundación CEIBA.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. Coordinates (seen as the component loadings of the PCA analysis) for the ten districts in the temporal span.