Temporal Topology Density Map

Spatiotemporal data visualization is of paramount importance in several applications, especially for analyzing urban mobility data encoded on topologies. This paper introduces a novel solution, called Temporal Topology Density Map (TTDM), to represent 2D discrete spatial data with temporal variations into a 2D continuous spatial space constrained by a topology. The proposed approach combines topological density maps with Change Frequency Heatmap (CFH) to convey visual information regarding temporal changes, leading to a new visualization method. Two case studies related to the analysis of response time associated with emergency services and walkability changes over time in specific areas of interest demonstrated the effectiveness of TTDM in challenging scenarios. The proposed solution provides an intuitive visualization for supporting the accurate analysis of spatiotemporal data changes over time using topology density maps.


I. INTRODUCTION
Advances in terms of processing power and storage capacity of devices have led to the creation of large data collections in various domains.In particular, a significant amount of spatiotemporal data is being generated and consumed, primarily because of the wide use of sensing technologies or the application of algorithms (data-driven or simulationbased).The proper analysis of these data, associated with identifying and understanding their relevance based on trends and patterns, is crucial for facilitating informed decision-making.In particular, it is essential to analyze and comprehend changes related to multiple variables (e.g., attributes) over space and time.
In the context of spatial data analysis, density maps are commonly employed to visually explore areas surrounding groups of points of interest [1].Such maps encode the continuous distribution of scalar fields in a 2D space.
The associate editor coordinating the review of this manuscript and approving it for publication was Haipeng Yao .
These approaches have been employed to visualize the projected spatial patterns associated with various data types.By dividing the data space into an arbitrary number of density fields and then visualizing them, distinct patterns can be discerned within specific datasets [2], [3].Density map visualizations have been successfully explored in several applications, particularly in the context of reasoning based on urban data.These approaches include support analyses related to traffic conditions [4], [5], pollution distribution [6], [7], and demographic evolution [8], [9].
Despite the success of existing methods in encoding spatial distributions (e.g., kernel-based methods [3], [10]), few initiatives have been dedicated to presenting density maps that incorporate topological information found in networks (e.g., road networks [11]).Recently, Feng et al. [11] introduced a promising approach for computing Topology Density Map (TDM).Their method extends the density estimation process from a 1D network to a 2D space, resulting in correct and intuitive visualizations.The correct visualization means an accurate reflection of the directional and topological road network, as well as the dynamic path costs directly linked to traffic conditions.The intuitive aspect refers to providing density fields in 2D planar fields in a clear and comprehensible manner.Figure 1, created using one digital twin tool [12], illustrates a density map (left) and a Topology Density Map (TDM) (right) associated with the analysis of spatial data.The color intensity in the figures is proportional to the density at each position.The TDM result on the right side of the figure depicts the 1D road network with colored tapered lines, which greatly assists users in visualizing density variations along with road segments.Representing the density intensity changes solely through density map calculation (left) can be challenging.
The algorithm proposed by Feng et al. [11] shows great promise and has been validated in the context of compelling applications related to urban mobility analysis.However, their solution did not account for changes in the topology density maps over time.In fact, the proper analysis of urban data requires not only an understanding of the spatial distribution of scalar fields, but also their temporal variation.In several applications, comprehending the trends and patterns of spatial data over time is a key element in supporting better-informed decision-making.To the best of our knowledge, the visual encoding of temporal changes associated with topology density maps remains a problem overlooked in the literature.
Several approaches have been proposed in the literature to represent and characterize change patterns [13], [14], [15], [16].Notably, Change Frequency Heatmap (CFH) [17] has been established as a promising method for encoding temporal changes using customized metrics and behavior patterns specific to an application.It has been successfully employed in the analysis of multivariate temporal data associated with plant phenology [17], [18], including data obtained from direct observation of individuals [17] over time and variations defined by sequences of images [18].
This paper introduces a novel visualization scheme for encoding temporal changes observed in topology density maps.Our solution relies on the computation of Change Frequency Heatmap [17] and its proper integration with topological information.
The main contributions of this paper are summarized as follows: 1) we address the problem of encoding temporal changes of topology density maps by introducing TTDM (its components and algorithms).2) we introduce a new approach based on the Change Frequency Heatmap for encoding changes in topology density maps.3) we propose a novel method for visualizing temporal changes associated with topology density maps.4) we investigate the use of CFH in the context of urban data analysis.5) we demonstrate the feasibility and usability of the proposed approach in the context of the analysis of real spatiotemporal urban (network) data based on the visualization of changes over time.The remainder of this paper is organized as follows: Section II briefly presents relevant background concepts related to our proposal.Section III provides an overview of the related work.Section IV introduces our method, that is, the Temporal Topology Density Map TTDM.Section V evaluates the proposed algorithm through quantitative and qualitative assessments.Section VI describes and discusses two use cases related to analyzing traffic conditions and walkability indicators over time using the TTDMs.Section VII concludes the paper and highlights the directions for future research.

II. BACKGROUND CONCEPTS
This section presents background concepts used in the proposed formulation for computing temporal topology density map (Section IV).First, a recently proposed method for computing topology density map [11] is introduced in Section II-A.Next, Section II-B presents the Change Frequency Heatmap and how this representation encodes temporal changes.

A. TOPOLOGY DENSITY MAP
A density map is an effective visualization method, that can present a continuous distribution of scalar fields in a 2D planar space by assigning a specific color to each scalar vertex.
Recently, Feng et al. [11] introduced a new method for computing density maps called Topology Density Map (TDM).The objective was to create more correct and intuitive density maps in the context of urban data visualizations.Their approach involved utilizing a directed acyclic graph (DAG) to propagate nonlinear scalar fields along 1D road networks.Next, the formulation extends the density calculation from the scalar fields to a 2D planar surface by identifying points-ofinterest (POI) nodes and computing density scalar fields for all points in this 2D space.
Figure 2 illustrates their pipeline, which comprises three modules: encode network data, compute accessibility data, and surface mapping.The first module encodes the network data, represented by a graph G = (V , E), where V denotes the node set and E denotes the set of directional edges.The node set V = {v 1 , v 2 , . . ., v i , . . ., v nodeNum }, where v i is a vertex and nodeNum is the number of nodes in V .To each node, it is assigned attributes, such as name, location, and whether it is a POI.The edge set E = {e 1 , e 2 , . . ., e i , . . ., e edgeNum }, where e j = (v x , v y ) is an ordered pair of vertices ((v x , v y ) ∈ V 2 , and v x ̸ = v y ) and edgeNum is the number of edges in E. A weight may be associated with the edge (v x , v y ), representing the cost from v x to v y , such as time, distance, or other customized scalars.In this research, the cost is typically the access time rather than the physical distance (e.g., Euclidean distance). 1he output of the first module is a wrapped structured DAGs data (G DAG ), computed from different POIs.Using the Dijkstra shortest path algorithm, the second module calculates a set of shortest path costs between each non-POI node and POI nodes based on G DAG .Then, it constructs the accessibility data G cost , which contains all the relevant information regarding accessibility for each non-POI node, including the shortest path cost from the nearest POI and the corresponding POI name.The accessibility data G cost serves as the input of the last module, surface mapping.This module is responsible for computing (G TDM ), which represents the estimated density field values for a 2D planar surface.

B. CHANGE FREQUENCY HEATMAP
Mariano et al. [17] presented a novel image-based representation, named Change Frequency Heatmap (CFH), which encodes the frequency of occurrence of temporal patterns associated with multivariate numerical data.
The CFH computation algorithm comprises three steps as illustrated in Figure 3: encode temporal data, compute temporal binary pattern, and compute histogram.In the first step, the temporal multivariate data S = {X 1 , X 2 , . . ., X n }, a set with n elements X i =< x i,1 ,x i,2 , . . ., x i,m > with m dimensions, is encoded.This step produces a stack of matrices M, where M =< M 1 , M 2 , . . ., M T >, and M t is an n × m numerical matrix at timestamp t ∈ [1, T ] composed of n lines and m columns.Each element m t of the matrix M t represents the pattern information at the corresponding position in the space S(n, m) and at timestamp t.The space S(n, m) = <(1, 1), (1, 2), . . ., (1, m), . . ., (n, 1), (n, 2), . . ., (n, m)> denotes the discretization grid of a continuous 2D space.The second module computes temporal binary patterns by applying customized functions that encode the temporal change profiles associated with different matrix cells.This step results in a new stack of matrices D, where matrix at timestamp t ∈[1,T-1] composed of n lines and m columns.Therefore, any cell d t in D t is at the same position as m t in M t .Equation 1 defines one example of an encoding function, where number 1 indicates a change across consecutive neighboring timestamps, while number 0 means the absence of temporal changes.
In the final step, the algorithm utilizes patterns of interest p within the sequences of numbers '0's or '1's (e.g., p 0 = 0, p 1 = 010, CFH h can be viewed as a matrix calculation algorithm that encodes the pixel change pattern of a stack of matrices and generates a binary motion histogram matrix D CFH .The operation CFH h counts the number of occurrences of the pattern of interest (binary pattern) in any element h of the matrix D CFH , defined as follows: where 1 ≤ x < y ≤ T − 1, y = l − 1 + x, and {d} is the set of all sub-parts of d.

III. RELATED WORK
This section provides an overview of the work related to our study.scalar field data in one 2D space with map rendering.Density map visualization is one promising approach to density-based field schematization to construct the visual coding of field data by density value, emphasizing geographic accuracy over visualization.
Due to their effectiveness in supporting space-oriented assessments, density maps have been successfully applied to urban analysis, especially concerning services for transporting people and goods.Danese et al. [19] explored the application of density maps for seismic risk analysis to improve civil protection planning.Xie andYan [20] presented a new density estimation method for traffic accidents and validated the visual effects with real datasets.Scheepens et al. [2], [3] proposed a novel framework for using density maps to visualize multiple attributes by multivariate trajectories and validated it in vessel traffic conditions analyses.More recently, Delso et al. [21] evaluated the density map application for the impact of the obstacles on pedestrian walkability assessment.Ren et al. [22], in turn, designed an interactive visual analytic system on the density map technique to demonstrate the effectiveness of an air pollution distribution assessment.
An effective density map computation method reflects the attribute value distribution by the density intensity.Because the data observation is limited to the nodes and edges of the topological network, the computation of density values for each position in a 2D space is the most critical process of density map visualization.This process is called density estimation.Various methods for density estimation have been proposed.For example, Borruso [23] validated the Quatic Kernel Density Estimation method in the urban immigrant population from the spatial distribution and tendency aspect.Krisp et al. [24] presented an adaptive directed Kernel Density Estimation (AD-KDE) to recognize the underlying dynamics of the vehicles for traffic condition analysis.Nie et al. [25] designed a kernel density method called NKDE-GLINCS, which integrates Network-constrained Kernel Density Estimation and Network-constrained Getis-Ord Gi* (GLINCS), to detect the abnormal status of road segments.Yuan et al. [26] proposed a new Quad-tree-based Fast and Adaptive Kernel Density Estimation (QFA-KDE) algorithm to compute the aggregation patterns more efficiently.
In the context of urban analysis, these density estimation methods can be distinguished into two categories.One category, known as planar KDE refers to the estimation of the density value based on the Euclidean distance between positions without any constrained network.The typical applications include the trail analysis of vessels [2], [3], trains [27], and flights [28].The other category, named NKDE, emphasizes the density estimation along a constrained network.It is widely adapted for the analysis of road-constrained events like traffic conditions, etc.
Topology Density Map (TDM) [11] (described in Section II-A) proposed by Feng et al., is a novel method integrating the advantages of planar KDE and NKDE.Their formulation utilized real traffic datasets and road networks to validate the effectiveness of the method in terms of density estimation in 2D space and a network.Compared with other research and existing methods [29], [30], [31] regarding spatiotemporal data analysis, TDM provides an intuitive visualization for supporting decision-making.
In our research, a new visualization method was explored for encoding and visualizing temporal changes in TDMs.This algorithm was validated in the context of mobility applications.Both formulations aim to leverage the added value of using density maps by decision-makers when exploring spatiotemporal datasets.Based on the available information, this study represents the first dedicated effort to encode the temporal changes associated with topology density maps.

B. TEMPORAL CHANGE VISUALIZATION
In the spatiotemporal context, temporal changes represent attribute value changes over time, while spatial changes are related to variations observed in distinct locations [32].Recently, Fang et al. [33] presented a comprehensive survey of the methods for time series data visualization.Their work divided most of the temporal change visualization methods into two categories.One is the visualization of time attributes, such as Spiral diagram, Calendar view, ThemeRiver view, and Dynamic visualization.The other is the visualization of highdimensional time series data, for instance, Parallel coordinate methods.However, all of them focus on the visualization of the changes in attributes over time, not on detecting the change patterns for multiple attributes simultaneously.
With the increasing collection of urban context data, it is challenging to analyze the temporal aspects associated with spatial data.One promising direction is to analyze and visualize the temporal changes in an explicit way.Krukowicz et al. [34] provided a comprehensive analysis of animal-vehicle road crashes by the temporal analysis with Calendar view and KDE of the number of the accidents.Liang et al. [35], in turn, explored temporal changes in population in a geographical area using an eigendecomposition method.Ziwen et al. [36] focused their work on the analysis of the spatial and temporal patterns for the tourist source market using the travel demand index pattern evolution method.Many other methods have been proposed to represent temporal changes in the literature [13], [14], [15], [16], but their applications are limited to specific scenarios or applications, i.e., they are not generic enough to be tailored to other applications.
In the context of urban data, Zheng et al. [37] summarized the existing visualization techniques from temporal, spatial, and other aspects.One of their focuses was on data exploration and pattern interpretation, which are highly relevant for decision-making.According to them, for multiple attribute visualization, pixel-based techniques interpreted by a matrix form are popular techniques.However, these techniques pose challenges when their use requires filtering some patterns of interest explicitly for all pixels simultaneously.
On the other hand, popular density estimation methods, such as KDE, have been extended to handle the temporal dimension.The spatiotemporal kernel density estimation (STKDE) method is an example [38].A predictive hotspot mapping to represent the risk factor by considering the temporal dimension is a sample application of this type of density estimations [39].Based on our current knowledge, none of them are suitable to visualize the temporal changes for TDMs.
To encode the temporal changes of multivariate attributes, Mariano et al. [17] proposed the Change Frequency Heatmap method (CFH) -see Section II-B.It may utilize ST datasets as a stack of matrices to produce a matrix recording the occurrences of the interested change patterns for each position and has been validated in plant phenology analysis.Their method has not been utilized for the visualization of time-related urban data.
In this work, the combination of CFH with topology density maps was explored and investigated.The main motivation relies on the fact that the computation of CFH is essentially a matrix-based operation.The temporal changes in the results of TDM are encoded using CFH with applicationoriented customized metrics and behavior patterns.While the traditional visualization of CFH is typically presented as a heatmap, this paper explores an alternative approach by encoding the CFH results into a 3D map view combined with density maps.

IV. TEMPORAL TOPOLOGY DENSITY MAP
This section describes the steps for the computation of the Temporal Topology Density Map (TTDM).

B. COMPUTE REPRESENTATIVE TOPOLOGY DENSITY MAP
This step aims to compute a graph to serve as a representative of the whole set G. This representative graph (G rd ) will be used as input to the TDM algorithm.The goal is to produce maps to be used for the visualization of the 1D network and the 2D space.Given the input M, this module processes all graphs in G based on a predefined function, effectively acting as a filter.The default function is a weighted computation based on edge values defined as M rd =(1/T ) t=T t=1 w t M t , where M rd is a matrix that encodes the representative graph G rd , w t are weights and (1/T ) t=T t=1 w t = 1.The default weighting scheme implemented assigns the same weight value for all matrices, i.e., w 1 = w 2 = w 3 = w 4 = 1.
The other module (Compute TDM), introduced in Section II-A, computes the shortest path cost and the source POI label for each node.The cost and label information can then be used to determine the accessibility of nodes (the first output).Also, this module outputs the cost and label maps that will be used later to produce the final visual representation related to TTDM.

C. ENCODE TEMPORAL CHANGES
To encode temporal changes, TDMs for each timestamp matrix in M need to be computed.This step produces sequences of label, cost, and density maps.The module Compute CFH utilizes any of these three maps as input.It then uses a metric f m and behavior pattern S bp to create a change frequency map, which reflects the number of times the behavior pattern occurs for each position of the map.In the following, these two steps are detailed.

1) COMPUTE TOPOLOGY DENSITY MAP (TDM)
This step concerns the computation of the TDM for each timestamp t, i.e., it computes TDM for each given input graph defined in terms of its edge cost matrix M t (t ∈ [1, T ]).The Compute Accessibility Data module in the TDM algorithm (Section II-A) generates the shortest path cost (Cost nr ) and 110600 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

2) COMPUTE CHANGE FREQUENCY HEATMAP
This step concerns the computation of the Change Frequency Heatmap, as described in Section II-B.This component receives a metric function that defines the changing pattern of interest as a parameter.Another parameter is the change binary pattern.The input of this module is a set of temporal maps (e.g., label, cost, or density maps), marked as After computing the temporal binary pattern with the metric function f m , the algorithm creates a group of change binary maps.Finally, it outputs one change frequency map CHM f with the number of occurrences of the given behavior pattern.
Algorithm 1 outlines the steps of temporal change encoding.It begins by calculating binary change elements d t for each location in one map, for example, M c 1 (Lines 3-11).The algorithm then constructs a binary change string d and computes the corresponding element value in the change frequency map CFH f (Lines 12-13).These calculations can be performed independently for all locations to obtain the output change frequency map CFH f , used for the visual representation.
Figure 6 shows a more detailed example of the temporal changes with label maps and colored nodes.In the example, the center region of the label maps is highlighted.The corresponding graph is presented in the middle of the figure for different timestamps.It is assumed that the color of H 1 is blue and the color of H 2 is red.At timestamp 1, the point of interest H 2 is the closest to node C, and H 1 is the closest  to other nodes.At timestamp 2, the nearest point of interest for node B changes from H 1 to H 2. At timestamp 3, node B is more accessible from H 1 again, while node D falls into the influence area of H 2. Finally, at timestamp 4, there are no changes in the nearest point of interest among the nodes.Counting these changes reveals two status changes for node B and one for node D. Therefore, the temporal variation value for node B is two, and for node D is one.

D. COMPUTE THE VISUAL REPRESENTATION
The module Compute the Visual Representation provides the resulting visualization related to the execution of the TTDM algorithm.As Figure 7 shows, it computes a 1D network with the module (Network Color Mapping) and 2D space with the module (Density Field Color Mapping).Both of them rely on the outputs produced by the module Compute Representative TDM (Section IV-B).Moreover, it utilizes the output of the module Encode Temporal Changes (Section IV-C) to calculate a height map in the module Height Mapping.The last module, Visual Integrator, interpolates the heights for the vertices in the 3D mesh with the given interpolation scale m and integrates all visual layers into a 3D representation.The final visualization combines the maps related to temporal changes and the representative density field map simultaneously.The following sections describe these modules.

1) NETWORK COLOR MAPPING
The network color mapping aims to create a visual representation of the 1D road network with colored nodes and tapered edges.Algorithm 2 outlines the key steps for network color mapping.The first step is to draw all the nodes with the associated POI color (Lines 2-5).The other step is to draw all edges (Lines 6-13) by computing the normalized cost values along the edges (function f norm in the algorithm) of the path cost separately for the start and end nodes first (Lines 6-11).After that, it draws the edge, considering the different width values along the edge (Line 12).The color of the edge is the same as that of the start node (Line 9).

2) DENSITY FIELD COLOR MAPPING
Density field color mapping utilizes the density map and the label map L to create the output color map L c .Algorithm 3 outlines the key steps for the density field color mapping calculation.For each pixel in the color map L c (Lines 2-6), the 110602 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Figure 8 illustrates the results of the density field color mapping and its integration with the network color mapping.It is calculated based on the average cost result of the timestamps.Therefore, we can find the color consistency of 1D network (network color mapping) and 2D space (density field color mapping).The region defined by ABCD (with light blue and red colors) is easily accessible through H 1 and H 2 while the region (H 1A, only blue color) is only accessible by H 1.

3) HEIGHT MAPPING
The TTDM algorithm encodes temporal variations related to the costs between the nodes or the accessibility of nodes to the POIs.As discussed in Section IV-C, those temporal variations are encoded in a change frequency map.If we choose the label maps as the input for encoding the temporal changes, the change frequency map reflects the source POI label changes for each pixel.The height mapping module aims to compute a height value for each pixel p.This height value of the pixel p is calculated based on its Euclidean distance to the nearest non-POI node using Equation 3.
where p is a pixel in the height map, K is the scaling factor, CHM f (p) is the change frequency map value in the same pixel position p, and d(N nr , p) is the distance between p and the nearest non-POI node N nr .Figure 9 is the height mapping result for the sample graph used, with H 1 and H 2 representing the POI nodes.The black marked circles represent the nodes with more frequent temporal changes.The pixels close to the nodes B (CHM f (B) = 2) and C (CHM f (C) = 1) have higher change frequency map values, which are proportional to the elevations there.It also provides the visual effects integrated with 1D network (network mapping) and 3D mesh (height mapping).

4) VISUAL INTEGRATOR
This module integrates the visualization by taking into account an ''interpolation scale,'' referred to as m.The number of vertices in the 3D mesh determines the complexity of the calculation, which defines the number of details encoded in the visualization.For the existing 3D mesh, let V orig be a set of mesh vertices composed of s v rows and s v columns, and s v be the number of mesh vertices.The result of mesh interpolation is a new set of mesh vertices W composed of s w rows and s w columns mesh vertices.The required 3D mesh may be created by the W directly.The relationship between s w and s v is defined as s w = s v + (m − 1) × (s v − 1), Figure 10 presents the typical visual results of the module visual integrator for the sample topology example.It combines multiple 3D element layers in one 3D space.For instance, the 1D network layer, 2D space layer, and 3D mesh layer represent the visual effects of the network color mapping, density field color mapping, and height mapping modules, respectively.For a 2D visualization of the TTDM, the 3D mesh layer was not used.

V. EVALUATION
This section addresses the quantitative and qualitative assessment of the proposed visual structure and associated algorithms.The overall running time of the visual structure computation is compared according to three aspects: number of POIs, resolution, and number of timestamps.We also assess the scalability of the proposed method by considering four topologies with different numbers of nodes and edges.Furthermore, the qualitative assessment of the produced label maps, cost maps, and density maps with various resolution scales indicates the appropriate choice of the scale parameter m.

A. EVALUATION PROTOCOLS
Four topologies (road network) are utilized in the performed assessments.Table 1 presents their main properties.Topology 1 is the simplest one and contains 6 nodes and 11 edges.This topology is also the primary running example discussed in the previous sections.Therefore, it has no corresponding actual network type and geographic region coordinates.The other two topologies are real networks, downloaded from OpenStreetMap. 2 Topologies 2 and 3 refer to geographic regions associated with a walkable road in the Ålesund (Norway) center.Topology 4, in turn, is composed of the main drive roads and intersections in Ålesund municipality.
With regard to the efficiency assessment, we compare the computational time under different conditions, such as resolution scales, different network sizes, and the number of POI nodes.The qualitative assessment relies on the comparison of cost, label, and density maps with various resolutions.We also assess the fidelity of density maps, i.e., to what extent a particular density map differs from a reference.The indicator (Percentage difference) %D is used to compare the current density map with the reference result, which is the density map ref when the highest resolution scale (m = 10) is used.The percentage difference %D for any pixel p in the map is the absolute value of the difference divided by the original value multiplied by 100 as shown in Equation 4.
The visual effects are encoded into a PGM P5 image. 3 It displays the grayscale colors with the maximum value (white color) and the minimum one (black color) for each pixel.The same format and method are also suitable for the label and cost maps.
There are no predefined POI nodes in the performance assessment.The experiments were carried out on a Lenovo Xiaoxin Air 14 (2021) 2.4GHz i5-1135G7 with one Nvidia GeForce MX450 graphs card and 16GB RAM memory.

B. RESULTS
2 presents a detailed comparison of the running times for our proposed method based on several parameters, including the number of timestamps, points of interest (POIs), image resolution, and topology complexity.The running time is about 30 minutes (1956.2s)for Topology 4 when we set 100 timestamps as well as the resolution of 360 × 270.To emphasize the effects of the first two aspects, we present two figures.Figure 11a shows that the computation time increases almost linearly with the number of timestamps, indicating that the calculation of the TDM for different timestamps consumes most of the time, rather than the CFH algorithm.Furthermore, the computation time is significantly influenced by the network complexity and the number of POI nodes, as Figure 11b shows.Those results suggest that the running time is less costly when the number of POIs increases in some cases.For instance, when we randomly select 50 and 100 POIs for Topology 3 (see light green bars), the corresponding number of non-POI nodes is 70 (130-50) and 30 (130-100).The experiment using Topology 2 (light blue bars) presents a decrease in the running time when the number of POIs increases from 2 to 20.
Given the running time, the scale value m = 10 is considered as the reference resolution for all these four topologies.Taking density maps of Topologies 1 and 2 as examples, Figure 12 utilizes a gray-level image to present the TTDM calculation result when m = 1, 4, 10 (highlighted in red, blue, and black borders).The larger figures (highlighted in pink and cray) present different percentages for m = 1, 4 compared with the reference resolution.Topology 1 has a higher difference around the POI nodes.For the real road network, under the various resolutions, a great difference is found in the boundary regions with a larger gradient value.Furthermore, m = 4 can greatly diminish the percentage difference for the borders compared to the image computed with m = 1.Similar results were observed for other topologies, and for label and cost maps.Therefore, we set m = 4 for the case study as it provides a good trade-off between the running time and qualitative assessment.

VI. CASE STUDIES AND DISCUSSION
This section presents compelling case studies concerning the use of TTDM in two different mobility-based applications.Section VI-A introduces the datasets used.Section VI-B describes the first application, which refers to improving the walkability of the relevant regions in Ålesund, Norway.Next, Section VI-C describes the second case study that concerns the speedup of the response time associated with emergency services also in Ålesund.Finally, Section VI-D discusses alternative designs and limitations of the proposed approaches.

A. DATASETS
This section presents the datasets used in the case studies.It refers to the real and simulated data.For both cases, we used the snow depth historical data at the weather stations located at Ørskog (SN60800), Eide På Nordmøre (SN62900), and Refvik (SN59250) from the Norwegian Climate Service Center. 4 It covers the time period from 01.01.2021 to 31.03.2023 and is separated into nine quarters (2021Q1-2023Q1).Any quarterly dataset can be filtered according to weekdays (Monday, Tuesday, Wednesday, and Thursday), weekends (Friday, Saturday, and Sunday), or all days (Monday -Sunday).By default, the data related to all days in 2021Q1 works as our primary data source.

2) SIMULATED DATA
Our research assumes that the average moving speeds of the people and vehicles vary linearly with the snow depth on the road, following the assumption by Wang and Liu [40].Therefore, when the snow depth is greater, the edge cost (time) for each edge increases at a slower average speed.
The snow data collected from the three weather stations were converted to raster data with spatial interpolation.Inverse distance weighting (IDW) [41] is a common method used for this and we utilized the implementation available in QGIS.It has been widely used in spatial interpolation of weather and air quality data [42], [43].The output raster image resolution is 80 × 21 using pixel sizes of 0.0036768475 • and 0.003505219 • for longitude and latitude, respectively.We assumed the snow depth D snow on this road would be computed using the average value of all pixels passed by, as defined in Equation 5: where p i is the snow depth on the target road in any pixel interpolated using the IDW method.
The average moving speed on this road S avg is calculated by Equation 6, where S ref is the reference speed, k r is the assumed rate of speed decrease, and S min is the slowest speed. 4https://seklima.met.no/(As of May 2023). 5https://qgis.org/en/site/(As of May 2023).S max_road is the maximum vehicle speed allowed on each road.Table 3 summarizes the parameter settings for different case studies.
The cost of a road was computed as the distance divided by the average moving speed S avg .This variation of these cost values based on the snow data will determine the temporal changes that will be encoded by TTDM.

B. CASE STUDY 1: WALKABAILITY ANALYSIS
This case refers to mobility analyses towards improving the walkability of relevant regions of a city.Specifically, in this case, we assess how weather conditions affect the walkability in different areas in Ålesund, Norway.Walkability is a representation of how easy it is to move in a path [44], [45]; thus we assume that the most walkable road/region is the one with the lowest access time, i.e., the path with the largest values in a density map.It is accepted that people will choose shorter paths, which may be especially important for those with reduced mobility, such as the elderly or mobilityimpaired persons.
In the considered scenario, high temporal changing scores associated with a particular location, which can be encoded in a change frequency map (Algorithm 1), indicate that the walkability of that region changes very frequently.In this situation, people may decide to change their destination more often or may hesitate to walk through that region.
Suppose that urban planners need to make a decision concerning the definition of the location of a new public service center in the region (it can be seen as a new POI).In this scenario, urban planners would prefer to ensure that smaller values of the Change Frequency Map (fewer temporal changes) are associated with target regions after the 110606 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.construction of the new facility.This is the kind of analysis for which the use of the TTDM is demonstrated in this section.
Figure 13 details the target area considered in the case.There are two existing indoor public service places considered POI nodes: H 1 (highlighted in blue) and H 2 (red).Planners may decide to retrofit some malls to build a new public service center and try to select a suitable location from two options.In the figure, they refer to H 3 (green) and H 4 (purple) as POI nodes.Therefore, they will choose between Solution 1 (H 1, H 2, H 3) and Solution 2 (H 1, H 2, H 4). Furthermore, they are more concerned with the impact in the sample area (highlighted with a blue dotted line border in the top-left corner).Three small core regions in this area will be explored in the analyses: Region The rest of this section is organized as follows.The overall visual effects will be presented with all POIs (H 1, H 2, H 3, H 4) at first.After that, we discuss the results of TTDM after computing the representative TDM (Section IV-B), which reflects the average walkability affected by weather conditions in Quarter 1, 2021.Next, TTDM (Section IV-C) results with different binary pattern strings are discussed.Finally, we discuss results related to the use of TTDM in the encoding of changes in label maps defined in terms of filtered data by date.

1) OVERALL VISUAL EFFECTS
Figure 14 presents the overall visual effects for the two solutions when using the density map for color mapping, and cost maps for temporal changes with the binary pattern string ''01'' and threshold of 60s.This analysis refers to the occurrence of walkability changes in a one-minute window.As we can observe, the differences in height are concentrated in the regions far from the POI nodes H 3 (green) and H 4 (pink).This means these regions have a higher frequency of cost changes.The variation in color reflects accessibility.The colors of pixels close to white indicates greater time costs for access.For example, the density values in Regions 1 ⃝, 2 ⃝, and 3 ⃝ in Solution 1 (Figure 14a) are significant larger than that in Solution 2 (Figure 14b).That means the population in these regions could have a shorter access time when the new public service center is built on H 3 instead of H 4. Therefore, Solution 1 is more suitable in general.

2) COMPUTING A REPRESENTATIVE TDM
Figure 15 shows the results of computing representative TDM for the defined sample regions.Region 3 ⃝, which has the most pixels highlighted in green, is affected only by H 3. This means the people in Region 3 ⃝ have the highest walkability to the new public service center in location H 3. Furthermore, the pixels and edges close to node N 2 are green in Solution 1, which means that people will probably choose to walk to H 3. These visual layouts further support the claim that Solution 1 ⃝ is more suitable.

3) ENCODING OF TEMPORAL CHANGES
The cost maps serve as the initial data for encoding temporal changes, specifically representing the variation in access time over a period.In this context, a value of ''1'' in the binary pattern string indicates that the consecutive day has experienced an increase beyond the defined threshold value of 60 seconds.Conversely, a value of ''0'' represents the opposite scenario, indicating a lack of significant increase.Figure 16 presents the TTDM layouts for two binary string

4) SUPPORT TO COMPLEX DECISION-MAKING
With the proper choice of dates (using filtering options), TTDM can also support decisions regarding the time period to operate the new public service center.For the weekdays and weekends illustrated in Figure 17, Solution 1 led to a more uniform distribution of changes across the different regions compared to Solution 2, for which we can observe less temporal change frequencies, i.e., the maps have bars with lower average heights.The differences in terms of high are more evident for intersection N8 in these two solutions.Temporal evolution regarding POI labels relates to how frequently the closest POI node changes.In this case, it encodes to what extent people located in these areas would change their preferred public service center.The computation of temporal changes regarding the density maps now takes label maps as inputs instead of cost maps.Figure 18 shows that road N4-N3-N7 with a higher height.Actually, changes happen 8 times for the pattern ''01'' in Quarter 1, 2021 (total 90 days).Solution 2 performs worse in these regions with higher change frequency values compared to Solution 1.That means that there are substantial temporal variations in walkability.

C. CASE STUDY 2: EMERGENCY SERVICES
In this case, the objective is to speed up the response times associated with emergencies in Ålesund, Norway.Ålesund, one of the most beautiful cities in Norway, is often selected as a destination for cruise ships.A large number of visitors may increase the demand for emergency services like ambulances, fire departments, or police.Since Ålesund is a small city, the emergency services are limited; thus, there is a risk that these services will collapse.Therefore, planning tools to arrange the existing resources of emergency services are critical for securing their response time.
An examination of the ambulance services provided in the Sentrum area of Ålesund is considered as a scenario for the use of TTDM. Figure 19  ⃝ is a bus terminal nearby.The change in weather conditions like daily snow depths (Section VI-A) affects the access time of the service by vehicles (Section VI-A2).Similar reasoning may be used to model the presence of tourists over summer conditioned by other weather conditions (e.g., temperature).The decreasing pattern of density value directly reflects increasing access time.
Initially, the computed overall visual effects will be illustrated, especially considering the area of interest.Subsequently, the representative TDM and the encoded temporal changes will be presented and discussed.At last, an analysis comparing the results related to different time periods and using label maps for temporal changes is presented.
110608 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

1) OVERALL VISUAL EFFECTS
Based on the snow depth data in Quarter 1, 2022, Figure 20 presents the overall visual effects for the whole road network and the area of interest.In this case, we consider the downtown area (zoom in), when using density maps for temporal changes with the binary pattern string ''01.''Recall that, even though the whole road network is large (Figure 20a), we could select any small geographic area like downtown to compute the TTDM.This downtown area is used to present the results related to this case.Figure 20b shows that both Region 4 ⃝ and Region 5 ⃝ have similar pretty light red colors, close to white.This means that using only the TDM algorithm is not enough to support the analysis of the scenario considered in this case.Anyway, we can observe that the average elevations in Region 5 ⃝ are significantly higher than the ones in Region 4 ⃝.
2) COMPUTING A REPRESENTATIVE TDM ⃝ for all cases.As we can observe, the representative map differs a lot from the ones related to the investigated dates.That suggests that change pattern analysis is relevant for a better comprehension of the impact of possible solutions.
When comparing the t86-t90's results, the density values for Regions 4 ⃝ remain unchanged between t86 and t87.Subsequently, the values decrease at t87, increase at t89, and reach at the minimum value at t90.Notably, the average value (refer to the figure) is similar to the t87 outcome; however, it can not reflect the temporal changes of density values.This sequence of images illustrates a case in which the assessment of changes over time might be challenging to be done visually as multiple maps would need to be compared simultaneously.Adopting the approach proposed in this paper to encode temporal changes makes it easier to analyze such variations and then make appropriate decisions accordingly, as illustrated next.

3) ENCODING OF TEMPORAL CHANGES
Figure 22 shows the TTDM results for various binary patterns ranging from ''010'' to ''01111110'', where the number of ''1''s represents the consecutive days with increasing access time.Figures 22a and 22b, illustrating short-term temporal changes, exhibit higher elevations in the targeted regions.This suggests that weather conditions lead to more frequent increasing access time for one or two consecutive days when compared to the frequencies observed for the other patterns.As a result, planners may consider allocating additional ambulance resources to speed up response times during 2022Q1 if cruise ships stay for one or two days.On the other hand, Figures 22c, 22d, and 22e depict similar height distributions, indicating fewer temporal changes over time.Notably, Figure 22f shows no elevation, indicating minimal long-term temporal changes.Overall, Region 5 ⃝ is associated with more frequent temporal changes in access time, suggesting a greater need for rapid emergency services in that area.

4) SUPPORT TO COMPLEX DECISION-MAKING
Here, we investigate the potential of using the TTDM in a complex decision-making scenario, involving the analysis of changes over quarters (season-based analysis).The seasonbased analysis is motivated by the fact that snow depth changes considerably across seasons.We focus on changes in quarters 1 (Q1) and 2 (Q2) as are the ones when snow occurs more frequently.We also investigate the impact of selecting different threshold values (60s as the low threshold and 300s the high threshold) for determining if a change occurred.The analysis is based on label maps, i.e., it takes into account how the influence areas of POIs change Figure 23 shows the TTDM results for Regions 4 ⃝ and 5 ⃝, over 2021-2023 period.We can observe that both regions have similar elevation distributions in those quarters over the years.Notably, Figures 23a (2021Q1, T =60s), 23e (2022Q1, T =60s), and 23i (2023Q1, T =60s) reveal more noticeable temporal changes in the first quarter than Figures 23c  (2021Q2, T =60s) and 23g (2022Q2, T =60s) for T = 60s using the data of the second quarter.However, for T = 300s, the results are similar to those shown in Figures 23j (2021Q2, T =300s) and 23d (2023Q1, T =300s).These outcomes demonstrate that weather conditions exert a more significant impact on the results than the regions' topology, while the threshold plays a critical role in the visual outputs.In short, the choice of the best solution depends on the threshold values and the specific time periods considered in the analysis.
Using label maps as the input for the TTDM facilitates the identification of regions with more changes of the closest emergency service point.Figure 24 presents the outcome, which indicates no changes for all quarters.However, the region marked with a dotted line exhibits some temporal changes, particularly in 2023Q1 (as Figure 24e).Given that no significant differences were observed for the two regions concerning these aspects, decision-makers could extend the region of interest to investigate the occurrence of changes in the adjacent areas.

D. DISCUSSION 1) ALTERNATIVE DESIGNS
TTDM incorporates elevation as the visualized dimension, enabling the simultaneous presentation of temporal changes along with the density field and topology.Alternatively, users have the option to compute color mapping based on the density field or the temporal changes, similar to a heatmap.This solution can be seen as a ''trivial'' integration of TDM and CFH. Figure 25 compares this heatmap approach with the TTDM solutions using data related to Case Study 1.The objective is to identify two positions with high-density values (highlighted with orange and light blue circles) that also exhibit significant temporal changes.The red and blue colors with alpha values represent the transition of temporal 110610 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.changes from high to low.With the heatmap solution, users need to switch between the density field data (Figure 25a) and the heatmap representation (Figure 25b) to locate positions with high-density values and observe temporal changes.In contrast, the TTDM approaches (Figures 25c  and 25d) streamline this process through the visual integrator (Section IV-D4), utilizing elevation as a distinct visual aspect.For our case studies, the 3D bar visualization (Figure 25d) is selected as the preferred design due to its enhanced  effectiveness in recognition compared to the original 3D mesh design (Figure 25c).
The also relies on TDM and CFH.In our method, however, their integration is not trivial.For example, in this paper, we elaborate on different aspects related to the kinds of change (label and cost maps) that are expected to be encoded.Also, we discuss strategies for the computation of the representative TDM.Furthermore, the paper includes a comprehensive description of the procedures to generate the final visual representation, including algorithms and running examples.Finally, the proposed method is validated by employing the assessment of the impact of different parameter settings and compelling usage scenarios involving decision-making based on real urban data.

2) LIMITATIONS
Through the evaluation of efficiency and insights from two case studies, the potential of TTDM to provide insights into temporal changes associated with density maps becomes evident.Nevertheless, the TTDM computation seems to be expensive, especially for large topologies.We claim, however, that the computation time is acceptable for realworld usage scenarios as those involving real topologies as the ones considered in the two case studies.Anyway, further research is needed to explore methods for improving the overall performance of the algorithms.
The practical use of the visualization method demands the proper definition of configuration parameters (e.g., threshold value -T , binary pattern -S bp , interpolation scalem, and scaling factor -K ), given the requirements defined in terms of the kind of analysis that will be performed and properties of the target dataset.The proper definition of such parameters might be a time-consuming and laborious task, especially for novice users.
To address this challenge, the utilization of digital twins [12] designed for TTDM may be used to support the analysis of different TTDM parameters' values with regard to the overall quality of the generated maps as well as the efficiency of the algorithms.These capabilities may lead to enhanced efficiency analysis, faster definition of parameter values, and potentially more insightful visualizations, ultimately supporting better-informed decision-making tailored to new datasets and topologies.

VII. CONCLUSION
This paper introduced a new method, Temporal Topology Density Map (TTDM), to intuitively encode and represent temporal variations associated with the topology density map.The proposed solution explores the Change Frequency Heatmap (CFH) that registers the occurrence frequency of change patterns of interest.TTDM provides the possibility for analysing temporal changes associated with density maps along the network, a problem overlooked in the literature.The method starts with a stack of graphs, with each graph representing a network with dynamic edge costs.After computing the representative TDM, the method encodes temporal variations represented in a 3D space.The efficiency of the method was evaluated based on different parameter values and topologies with different sizes.The proposed method was also validated in two compelling case studies related to urban planning activities for assessing walkability and emergency services analysis based on real road networks.
This conducted research opens opportunities for future work in several directions.First, we can explore TTDM in other applications for analysis scenarios involving spatiotemporal data associated with topologies.One promising application is to construct sensor monitoring systems whose spatial distribution relies on a wired or wireless communication network.In this system, the network includes data collection points, such as sensors, which are connected to data centers.The sensors represent non-POI nodes, whereas the data centers serve as POIs.Edge cost is the communication delay.Examples of time-delay-sensitive scenarios include applications involving underwater acoustic sensor networks [46], smart grid-connected power systems [47], nonlinear time delay systems [48].Another research direction involves the investigation of suitable domain-specific CFH configurations.Because the configuration of CFH, including the definition of the metric function and binary patterns of interest, is application-dependent, it would be worthwhile to use packaging them as libraries for future re-use.For instance, the configuration used in the two case studies provided in this research could serve as the first two templates for mobility applications.Users can benefit from them in the assessment of different metric functions and change patterns.

FIGURE 1 .
FIGURE 1. Examples of density maps.The left figure is a density map, while the right one is a Topology Density Map (TDM).

FIGURE 3 .
FIGURE 3. Change Frequency Heatmap algorithm pipeline with a running example.The temporal data S includes a 4 × 4 discretization grid for one limited 2D space.The metric function calculates the changing pattern as Equation 1 for neighbor matrices Ms and outputs into a new stack of matrices D. For example, D 1 is computed on M 1 and M 2 .The pixel value of the final result matrix D CFH encodes how many times a pattern of interest occurs.For instance, D CFH (0) has the number of '0' for each pixel and D CFH (010) counts the pattern change '010.' The figure is adapted from [18].

FIGURE 4 .
FIGURE 4. Temporal topology density map algorithm pipeline in concept.

Figure 4
Figure 4 illustrates the pipeline for the TTDM computation.The input data consist of a stack of graphs, each associated with a different timestamp.Two main branches use this stack of graphs as input.The first one goes through the Compute Representative TDM component and produces two outputs: the network accessibility data related to the 1D network and the density, label, and cost maps associated with the 2D space.Section IV-B describes how this step works.The other branch goes through the component Encode Temporal Changes and creates a change frequency map (matrix).Section IV-C provides the details of this component.The three outputs of

FIGURE 5 .
FIGURE 5.An example of computing TDM.It reflects the temporal changes of the sample graph data when T = 4.

FIGURE 6 .
FIGURE 6.An example of temporal variation mapping.

Algorithm 2
Network Color Mapping 1 Auxiliary Data functions: f norm is a function that returns the linear normalization result of a given node path cost value in the cost set of G cost .The functions f drawNode and f drawEdge are used to draw nodes and edges with the given parameters on the Network Map G c .The variables s v x and s v y are the edge widths at v x and v y .c e is the color of the edge.The edge width scaling factor is K w (default to 1), which controls the size of the edge width.It determines the extent to which the edge width is scaled.Data: Accessibility data G cost Result: Network Map G c includes colored 1D road network 2 foreach vertex v in V do 3 c v ← the color of node v by G cost 4 f drawNode (v, c v ) 5 end 6 foreach edge e = (v x , v y ) in E do 7 F c (POI nr , v x ) ← the path cost of node v x from belonging POI node by G cost 8 s v x ← 1 − f norm (F c (POI nr , v x )) 9 c e ← the color of node v x by G cost 10 F c (POI nr , v y ) ← the path cost of node v y from belonging POI node by G cost 11 s v y ← 1 − f norm (F c (POI nr , v y )) 12 f drawEdge (e, s v x , s v y , c e , K w ) 13 end

Algorithm 3 FIGURE 8 .
FIGURE 8.An example of density field color mapping.The density field color mapping result (left) and the integrated result with the network color mapping (right).

FIGURE 9 .
FIGURE 9.The height mapping result (left) and the integrated result with the network color mapping (right).

FIGURE 10 .
FIGURE 10.The visual integrator example with the multiple layers.

FIGURE 11 .
FIGURE 11.TTDM computation time results.Figure (a) is the result with m = 4 and the number of POIs = 2. Figure (b) reflects the computation time for various numbers of POIs in the same resolution m = 4 and the number of timestamps = 100.

FIGURE 13 .
FIGURE 13.The regions of interest considered in Case Study 1.There are three small regions ( 1 ⃝ -3 ⃝) considered in the analyses.

FIGURE 14 .
FIGURE 14. Overall visual effects for the two solutions.

FIGURE 15 .
FIGURE 15.The results of the sample regions after computing the representative TDM for the two solutions.

FIGURE 16 .
FIGURE16.The TTDM results related to the encoding of temporal changes for two solutions considering binary pattern strings ''010'' and ''0110''.

FIGURE 17 .
FIGURE 17.The comparison of TTDM results for specific types of days considering two solutions.

FIGURE 18 .
FIGURE 18.The TTDM results of the sample regions using label maps as input for the two solutions.Solution 2 leads to significantly higher bars on road N4-N3-N7.
shows the three health services that provide ambulances.These are H 101 (highlighted in blue), H 102 (red), and H 103 (green).H 101 and H 102 are close to the downtown area of interest(latitude from 62.4691749 • to 62.4729876 • and longitude from 6.1496012 • to 6.1605936 • ).From a planning perspective, city planners and emergency coordinators would like to analyze whether they should provide more resources at the time when one or more cruise ships arrive in Ålesund.In this case, the chosen areas of interest include Region 4 ⃝ (N 101, N 102, N 103) and Region 5 ⃝ (N 102, N 104, N 105).Region 4 ⃝ is the port position for cruise ships and Region 5

FIGURE 19 .
FIGURE 19.Regions of interest in Case Study 2. There are two small regions ( 4 ⃝, 5 ⃝) which are connected by node N102.

FIGURE 20 .
FIGURE 20.Overall visual effects for Case Study 2.

FIGURE 21 .
FIGURE 21. results of computing a representative TDM (a) for the average of density maps and the TDM result for interested timestamps (b)-(f).

Figure 21 presents
Figure21presents the representative TDM (Figure21a) as well as the TDM result for the last five days in the period of analysis considered (Figures21b-21f).In this figure, txx represents the xx-th day.No significant differences over time were observed in Regions 4⃝ and 5 ⃝ for all cases.As we can observe, the representative map differs a lot from the ones related to the investigated dates.That suggests that change pattern analysis is relevant for a better comprehension of the impact of possible solutions.When comparing the t86-t90's results, the density values for Regions 4⃝ remain unchanged between t86 and t87.

FIGURE 23 .
FIGURE 23.The comparison of various quarterly time periods with threshold T = 60s and 300s.

FIGURE 24 .
FIGURE 24.The comparison of various quarterly time periods using label maps for TTDM CFH computation.

TABLE 1 .
The road networks considered in the evaluation and case studies.Coordinates are encoded in degrees.

TABLE 2 .
The running time (in seconds) comparison for TTDM by different number of POIs, image size resolution, and topology complexity.

TABLE 3 .
Parameters for simulated temporal changes by snow depth.Topology 3, introduced in Table 1, is selected as the road network for Case Study 1.We retrieve four public service places as the POI nodes: H 1, H 2, H 3, and H 4. For Case Study 2, we choose Topology 4 as the road network and three locations (H 101, H 102, and H 103), associated with emergency services facilities in Ålesund, Norway, as POI nodes.