VASA: an exploratory visualization tool for mapping spatio-temporal structure of mobility – a COVID-19 case study

ABSTRACT Exploratory data analysis tools designed to measure global and local spatial autocorrelation (e.g. Moran’s $I$I statistic) have become standard in modern GIS software. However, there has been little development in amending these tools for visualization and analysis of patterns captured in spatio-temporal data. We design and implement an exploratory mapping tool, VASA (Visual Analysis for Spatial Association), that streamlines analytical pipelines in assessing spatio-temporal structure of data and enables enhanced visual display of the patterns captured in data. Specifically, VASA applies a set of cartographic visual variables to map local measures of spatial autocorrelation and helps delineate micro and macro trends in space-time processes. Two visual displays are presented: recency and consistency map and line-scatter plots. The former combines spatial and temporal data view of local clusters, while the latter drills down on the temporal trends of the phenomena. As a case study, we demonstrate the usability of VASA for the investigation of mobility patterns in response to the COVID-19 pandemic throughout 2020 in the United States. Using daily county-level and grid-level mobility metrics obtained from three different sources (SafeGraph, Cuebiq, and Mapbox), we demonstrate cartographic functionality of VASA for a swift exploratory analysis and comparison of mobility trends at different regional scales.


Introduction
The COVID-19 pandemic affected our communities in a variety of ways: from ubiquitous telecommuting, to statewide lockdowns and tier-based contagion mitigation strategies.The shifts in mobility patterns over the course of the pandemic have been unprecedented, and we still observe shifts in daily mobility compared to the pre-pandemic era (Li & Giabbanelli, 2021).At the same time, never before have researchers had access to so much fine-grained mobility data to closely investigate the effect of nonpharmaceutical interventions (NPI) on mobility (Iacus et al., 2020), as well as their effectiveness in containing the virus (Haug et al., 2020;Hsiang et al., 2020).The data provided by Google, Facebook, Apple, SafeGraph, Cuebiq, and others (Noi et al., 2021a) have opened new avenues for mapping and understanding the compound mobility indices and complex structures of human mobility across time, space, and scale.While there are many ways to assess and model the spread of the coronavirus, the Geographic Information Systems and Science (GIScience) are crucial in enabling the analysis (Rosenkrantz et al., 2021) and are frequently used for disease mapping (Lyseen et al., 2014).
Various metrics capturing human mobility and potential exposure have been used to establish the connection between human movement and virus transmission (Kraemer et al., 2020;Nouvellet et al., 2021;Xiong et al., 2020).These data sets describe aggregated mobility and are compiled from a variety of sources: GPS, cellphone units, call-detailed records (CDR), closed-circuit television (CCTV), logs, checkin records, and transaction data.More importantly, these mobility metrics can provide us with a better understanding of the COVID-19 pandemic, when coupled with powerful visual analytics (Leung et al., 2020;Reinert et al., 2020;Sahnan et al., 2021).In particular, only a limited number of tools exist that allow simultaneous spatio-temporal analysis of associations at different scales.
The task of mapping a dynamic spatial process across various levels of aggregation is non-trivial.One must select the appropriate cartographic symbology and techniques for effective mapping (Andrienko et al., 2003;Vasiliev, 1997).Through the course of the pandemic, the use of dashboards generated a lot of interests among scholars, policy-makers, and the

Geovisualization for COVID-19
The role of GIScience in spatio-temporal analysis, epidemic modeling, and mapping has been enormous during the ongoing pandemic (Franch-Pardo et al., 2020).In fact, the geographic maps of Hubei province were the first visual display of the soon-to-be COVID-19 epidemic (Guan et al., 2020).What came next was the avalanche of user-generated visualizations to help understand and curb the spread of SARS-CoV-2 virus, a number of which suffered from major limitations from problems with projections, classification schemes, and scaling issues (Field, 2020).The academic community also contributed immensely 2 , developing maps and visual analytics tools to help understand the spatial structure of the coronavirus (Boulos & Geraghty, 2020).Several (geo)visualizations have been devised and presented to find more intuitive and efficient visual displays for mapping COVID-19 (Lan et al., 2021), including glyph maps (Beecham et al., 2021); bivariate, spike, and space-time cube (Lan et al., 2021), hotspot cluster mapping (Hohl et al., 2020), time geography (Yin et al., 2021), risk mapping (Loo et al., 2021), and flowmaps (Ponce-de Leon et al., 2021).

Exploratory spatial data analysis
Considering an infectious disease as a spatio-temporal phenomenon that spreads in space and time, we can examine the underlying structure of its prevalence by looking at how different factors co-vary at different locations.Various spatial autocorrelation metrics have been proposed to measure the amount of spatial dependence, including Getis and Ord's G (Getis & Ord, 1992), Ripley's K (Ripley, 1976), Geary's C (Geary, 1954), Matheron's Semivariance γðhÞ (Matheron, 1963), and Moran's I (Moran, 1948).While global metrics of spatial association are good at establishing the magnitude of the autocorrelation, local measures can pinpoint locations, where it is strongest or weakest.For instance, Local Indicators of Spatial Association (LISA) (Anselin, 1995) have been adopted from spatial econometrics and applied across multiple domains in exploratory spatial data analysis (ESDA) with numerous case studies related to COVID-19 (Murgante et al., 2020;Noi et al., 2021a;Ramírez-Aldana et al., 2020;Scarpone et al., 2020;Tokey, 2021;Xie et al., 2020).In addition to monitoring trends and hotspot detection, local indicators of spatial association can be used to investigate the effect of governmental non-pharmaceutical interventions (e.g.isolation, contact tracing, physical distancing, school and business closures) on the transmission of the disease (Haug et al., 2020;Hsiang et al., 2020;Kucharski et al., 2020).

Mapping spatial association
From a cartographic perspective, the maps generated from LISA typically illustrate statistically significant local clusters (coldspots and hotspots) and local outliers using a choropleth map.While these visual displays work really well for purely spatial data, developing analogous uncluttered visual displays for mapping spatiotemporal data (e.g.mobility data) remains an active area of research.For example, differential Local Moran (Anselin & Rey, 2014), which assesses the standardized differences between the original values of the variable, collapses the observation period to two points in time, while the emerging hot spot analysis (ESRI, 2021) tool breaks down the uncovered clusters into 17 types based on their spatial and temporal proximity, making the maps cluttered and hard to assess.
Inspired by traditional LISA visualizations, our proposed tool, VASA (first introduced in Noi et al., 2021b), seeks to overcome these limitations by providing cartographic representations of spatial clusters over time.These amended spatio-temporal views simplify the assessment of spatial associations, making exploratory spatial data analysis more accessible.VASA further adds flexibility that allows the user to run an analysis at different temporal and spatial levels of aggregation.As an exploratory visual analysis tool, VASA is designed to synthesize information from massive mobility data into easily interpretable insights.

Describing spatio-temporal patterns
Throughout the paper, we characterize both our data and induced spatial clusters in a number of ways.We use the term "structure" to refer specifically to data and how it is structured, as well as describe its attributes.On the other hand, we use the term "patterns" to describe the extracted spatial and/or temporal distributions.Finally, we utilize the word "trend" to characterize patterns of mobility with respect to time and temporal progression of patterns throughout the observation period.

VASA: visual analysis for spatial association
VASA consists of three modules as illustrated in Figure 1: import module, processing module, and visualization module.The import module contains various routines for managing spatio-temporal series and construction of adjacency matrices (Rey & Anselin, 2007).
The processing module enables aggregation of data to various levels of analysis: either spatial (e.g.tract, county, state, grid cells) and/or temporal (weeks, months, etc.) units.The visualization module handles the graphic object representations and provides plotting functions for two multivariate cartographic visualizations: a stacked recency and consistency map, and a linepath scatter plot.Both techniques use LISA (Anselin, 1995) as the base, and utilize the Local Moran's I and permuted p-values.The techniques are best suited for areal data (e.g.county, Census blocks, Census tracts, etc.) capturing aggregated information at different temporal and spatial levels.
For illustration purposes, consider a typical data set of mobility metrics, where the data is provided as spatiotemporal series at the county-level for each day of the 2020 in the USA.The stacked recency and consistency map allows to ascertain the spatio-temporal structure of mobility for the entire data set (i.e.national level) or a sub-sample of locations (e.g.state level).The line-path visualization enables a detailed analysis of patterns at the individual or local level (e.g. as counties within a specific state) over time.For example, the trends in counties identified as hotspots/coldspots for the state of California can be explored.The algorithms used to generate each of these visual displays are described in the initial conference paper Noi et al., 2021b, and the supplemental reproducible code is provided in GitHub.Here, we focus on describing the cartographic details, and provide a number of illustrative examples to demonstrate the applicability of the tool.More technical information on the visualization routines can be found in the Appendix.

Importing and pre-processing data
The user interface of VASA includes a component for data import and processing, which accepts tabular data and spatial data supported via Pandas (McKinney, 2010) and GeoPandas (Jordahl et al., 2020).During the import, the analyst specifies the spatial and temporal units for aggregation.For instance, when importing daily time-series data at the county-level, one might group the data by "state" (spatial aggregation) and/or "week" (temporal aggregation).This step is crucial as it will define the regional scale of analysis (e.g.national/ state level).At the same time, VASA runs a series of checks to identify and handle missing values in the prepossessing step, which could potentially cause issues in further analysis (see Section 5 for discussion of missingness in data).Currently, VASA provides three options for handling missing data.First, the overall average can be used to impute the missing values.For example, for county-level data that is aggregated at the weekly level, the arithmetic average for the entire sample of the county will be used for the week with missing values.Second, a rolling average can be used to infer the missing data for partially missing data points within a specific spatial unit (e.g.county).Finally, there is a method to count the percentage of missing data per each spatial unit and an option to filter the data if it is above a specific user-defined threshold.
After data preprocessing, the adjacency matrix is constructed to run LISA analysis.Following the recommendations by Anselin and Rey (1991); Farber et al. (2009); Florax and Rey (1995); Stakhovych and Bijmolt (2009), VASA uses the first-order Queen's contiguity to infer spatial weights by default, although other spatial weights construction methods are available including distance bands and k-nearest neighbors.Queen's contiguity provides a sufficient approximation of the true weight specification (Smith, 2009).The readers are referred to the Discussion section and Appendix for more information.For the constructed weight matrices, Local Moran's I are calculated for each temporal unit and each of the spatial units is classified as either a local hotspot or a local coldspot.The result of this process is a classification where every spatial unit is denoted as either a hotspot, a coldspot, or not-significant observation for each temporal unit over the observation period.It is important to add that while traditional visualizations of LISA include both the clusters (hotspots/coldspots) and outliers, VASA uses clusters only.The reason is to simplify visualizations and control the number of visual variables in the proposed displays.The classification is used as an input to the visualization techniques outlined below.Mathematically, the Local Moran's I coefficient is calculated for each unit of analysis: where w ij is the spatial-weight matrix, X and � X is the sample mean.Thus, calculated Local Moran's I coefficients will depend on the specification of spatial weights matrix and the sample mean.If our sample is the entire United States, then Local Moran's I will be calculated in relation to the national average.We can as easily calculate Local Moran's I values at the state level, by only including mobility metrics for a specific state.This is further demonstrated in the case studies in Section 4.

Stacked recency and consistency map (RECO)
The main geographic visual display available within VASA is called recency and consistency map (RECO).For illustration, Figure 2 schematically shows how the cartographic visualization is created for county-level analysis of data in South Dakota.First, the counties are depicted using their centroids (Figure 2a)).After running LISA, the centroids (shown as circle markers) are colored in blue if the associated counties are coldspots at the selected point in time, otherwise if the counties are identified as hotspots, the centroids are colored red (Figure 2b)).The centroids of the counties that are not significant clusters remain invisible.The color value (darkness or lightness) represents the recency (observations toward the end of time series) of hotspots/coldspots in the time series (Figure 2c)).A darker color shows that a statistically significant clustering behavior was observed late in the observation period (more recent), while a lighter color denotes a statistically significant clustering behavior early in the observation period (less recent).The visual variable size is used to represent the consistency of hotspots and coldspots (Figure 2d)).The consistency reflects the frequency at which a county became a significant cluster (hotspot/ coldspot) through the observation period.The radius of the centroid marker shows the total number of times that county was either a hotspot or coldspot in the given period (i.e.larger markers indicate more frequent clusters).

Line-path scatter plot
The second visual display available within VASA is called line-path scatter plot, which provides a temporal view and an aspatial data representation.The visualization process is illustrated schematically in Figure 3 for county-level analysis in South Dakota.Each county is represented as a line, denoting a clustering behavior of a county over time: the cumulative number of hotspot or coldspot classifications is mapped onto y-axis and the time is mapped onto x-axis (Figure 3a)).The line is constructed by connecting the points representing a cumulative number of times the county was classified as either a hotspot or a coldspot at different points in time.For instance, for a county that was a hotspot during weeks 5, 9, 10, and 15, the following points will be generated and connected: Point 1 (5,1), Point 2 (9,2), Point 3 (10,3),and Point 4 (15,4).Once the line segments are connected into a trajectory, the endpoints (the last time the county was a statistically significant hotspots/coldspot) are depicted with a circle marker (Figure 3b)).Next, LISA outputs are used to differentiate hotspots and coldspots: hotspots are colored in red and coldspots are colored in blue (Figure 3c)).Finally, the line endpoints are colored based on the majority rule: for counties that were consistent hotspots -red, otherwise -blue (Figure 3d)).The counties that were classified as hotspots/coldspots once or were not classified at all are not rendered on the plot.
The clustering consistency for different spatial units depends on multiple factors and varies depending on the assessed data set.Three basic types of behavior are common: a consistent hotspot (i.e. a county was classified as a hotspot only), a consistent coldspot (i.e. a county was classified as a coldspot only), and a flipper (i.e. a county was identified as both a coldspot and a hotspot within the observation period).Flips are mapped on the main line-path scatter plot with an empty circle marker and a thick edge endpoint (Figure 3d)).We also added the total tally for these three types of behavior as a text box in the top-right corner.
Additionally, a variation of the line-path scatter plot, which only looks at flipping behavior was developed (see Figure 4).In these plots, the color of the line changes to reflect the clusters' most recent hotspot or coldspot behavior.The week at which the cluster changes behavior is marked by an "X" on the line.The endpoint marker for each cluster can either be a circle denoting a hotspot/coldspot classification or a flip ("X" marker).

Utility functions
We added several utility functions to VASA to complement exploratory analysis: jitter and sampling.In some tested applications, large sample size may result in layering of lines with similar behavior for line-path scatter plots.For counties whose endpoints terminate at the same coordinates, the jitter (i.e.minute variation of coordinate values) was added to differentiate such cases.Sampling functionality is intended to simplify information overload for large sample sizes.The sampling is used to "thin out" the visual display by randomly selecting a subset of spatial units (e.g.counties) and can be adjusted by the user by setting a sampling ratio (i.e. percentage of observations to keep in a sample).Sampling can be applied to both line-path scatter plot and recency and consistency map.

Mapping spatio-temporal information to visual variables
Originally summarized by Bertin (2010), visual variables describe how the information can be encoded into graphical dimensions, describing basic building blocks of any visualization.Bertin believed that visual variables are assessed in an immediate (preconceptual) manner, and not understood cognitively, yielding different types of processing and perception: association perception, selective perception, ordered perception and quantitative perception.We followed general design principles and recommendations for visual semiotics and information visualization, when designing visual displays for VASA (Hegarty, 2011;Slocum et al., 2009).In particular, we used color hue (an associative and selective variable) to distinguish hotspots and coldspots clusters to put an equal weight on the two types of clusters, allowing the user to perceive variations in color hue as two distinct groups and visually isolate these two phenomena.At the same time, color values (a dissociative, selective, ordered variable) were used to make the outlined patterns over time (e.g.clusters of hotspots and coldspots) more discernible, such that the user's eye is attracted to this variation.Finally, the size (selective, ordered, and quantitative variable) allows the user to assess the consistency and frequency of a clustering behavior through the notion of "more" or "less."The recency denoted by the color values is also "ordered" and enables the user to differentiate between "more" or "less" recent clustering behavior.Finally, time is captured in multiple ways: RECO represents time using color values (the darker the more recent cluster).It also provides a dynamic representation of RECO, which can animate the development of the clusters over time.Finally, the line-path plots provide a timeline view, representing the changes in the clusters over time.

Case study
This section describes two case studies to demonstrate the capacities of VASA in assessing mobility patterns in response to the COVID-19 pandemic through exploratory data analysis and visualization at multiple geographic scales (e.g.county-level, state-level, and 100 m grid-level).The first case study applies VASA to investigate and compare county mobility patterns at the state and national levels using data structured as polygonal layers, while the second case study shows how VASA visual displays can be adapted to assess data that are structured as grid or other forms of tessellation.Leveraging data obtained from multiple sources (SafeGraph and Cuebiq, Mapbox), the case studies aim to 1) investigate and compare mobility trends at different regional levels; 2) explore the intra-urban spatial heterogeneity of human mobility.Since LISA is essentially a cross-product statistics between a spatially lagged variable and a variable itself expressed in terms of deviation from the mean, by changing the scope of the analysis (national versus state), we can explore the structure of mobility in relation to a different point of reference.For case one, two regional scales of analysis are used: national and state.For analysis at the national level, all of the available data from 3,200 counties is used to calculate the Local Moran's I, while for the state-level analysis only the data for counties within a specific state is used to generate the Local Moran's values.Finally, in case study 2 using Mapbox grids, only the data for the greater Los Angeles area are used.

Mobility data sets and mobility metrics
The case studies utilize data from three different sources: SafeGraph, Cuebiq, and Mapbox for 2020.SafeGraph and Mapbox provide mobility data openly through an initiative Data for Good, while Cuebiq provides the data on a commercial basis and requires license fees.While SafeGraph and Cuebiq publish several metrics of aggregate mobility at different spatial and temporal granularity, we utilize the following two metrics as detailed in (Noi et al., 2021a): Cuebiq % sheltered in place quantifies the percentage of users staying at home in any given state or county.It is calculated daily by measuring the number of users moved less than 330 feet from home.SafeGraph % sheltered in place is calculated as: completely home device count=device count, where completely home device count represents the number of devices that did not leave the geohash-7 (153 m � 153 m) area in which their home is located during the aggregation time period.These data are structured in aggregate form at the county level for each of the 366 days in 2020, representing the percentage of population staying at home across all geographic units.The differences in how the two sources define and estimate "home," as well as measurement biases inherent from smartphone data, are likely to influence the analytical outcomes and inferences as discussed in (Noi et al., 2021a).The Mapbox activity index reflects the level of movement activity based on cell-phone records on a given day within a spatial unit (OpenStreetMap tiles: 100 m by 100 m).The index value ranges from 0 (no activity) to 1 (high activity) and is normalized by a scaling factor calculated from the "99.9th percentile of anonymized counts across all geographic units and days in the baseline time span" (Mapbox, 2022).

Case study 1: comparing pandemic mobility responses in Georgia and South Dakota using SafeGraph and cuebiq data
This case study focuses on the states of Georgia and South Dakota because these states exhibited distinct geographic differences (see Figure 5) stemming from regional mobility patterns.They also exhibited high levels of spatial heterogeneity and pronounced differences between the mobility metrics across the two data sets.We use social distancing metrics representing the percentage of population staying at home obtained from SafeGraph and Cuebiq.Therefore, the inferred hotspots (described in the following sections) can be interpreted as counties with higher percentage of people sheltering at home, located next to other counties with higher percentage of people sheltering at home.In contrast, the coldspots are interpreted as counties with lower percentage of people sheltering at home, located next to other counties with lower percentage of people staying in.

Mobility patterns in Georgia
For Georgia, the recency and consistency map is provided in two views: Figure 6a) captures induced clusters comparing the county data to the national average, while Figure 6b) captures clusters deduced by comparing the county data to the state average.Atlanta (denoted as a set of contiguous beige county polygons) exhibits a hotspot behavior (more people staying at home), while the rest of the state exhibits a coldspot behavior (fewer people staying at home).This makes sense, since large metropolitan areas are more densely populated and there is an increased risk of transmission with every trip (Hamidi et al., 2020;Souch et al., 2021).Another confounding factor could be an increasingly democratic political affiliation (Grossman et al., 2020) for the majority of the population in Atlanta metro area, which are generally more supportive of non-pharmaceutical interventions.While both SafeGraph and Cuebiq agree on this divide, the "amount of agreement" varies: SafeGraph classifies fie counties around Atlanta as recent and consistent hotspots (i.e.darker and bigger markers), while for Cuebiq, only two counties are identified as hotspots, with much lower recency and consistency (smaller lighter red markers).There seems to be a lot more agreement between the two data sources with respect to coldspots recency and consistency, particularly in the Northeast and the Southeast of the state (as indicated by similar color values and marker sizes).
Overall, assessing data in relation to the national mobility average helps to uncover the underlying gradient in mobility responses between urban and rural areas in Georgia.This effect is not necessarily discernible for every state (for instance in South Dakota as described in the next section), and may be attributed to the Modifiable Areal Unit Problem (Fotheringham & Wong, 1991) and how well the urban-rural gradient is fitted into the administrative boundaries where aggregation occurs.
On the other hand, the spatial configuration of clusters produced with data from Georgia counties (using the state average) is more uniform across the two data sources (Figure 6b)), albeit with slight differences in recency and consistency levels.For SafeGraph, the coldspots are less recent and consistent (smaller markers with lower color intensity).More generally, we observe that the state is segmented into two parts: the Northwest, where the counties with high percentage of sheltered people are clustered with other counties with high share of sheltering people; and the Southeast, where the counties with low percentage of sheltered people are clustered with other counties with low share of sheltering people.Once again, this is most likely attributed to the urban-rural gradient, but the effect of this gradient is less attenuated because we are only looking at 159 counties in Georgia.These outcomes suggest a more homogeneous mobility pattern across the counties within the state, which may be indicative of the state-wide NPI policies.
All RECO maps point to the differences in socialdistancing behavior across the urban-rural divide: with urbanized areas exhibiting higher levels of social distancing.The hotspots identified in relation to the national average are not as consistent as other hotspots in large dense metro areas in different states (e.g.Chicago, New York).This is likely caused by the sprawling type of urbanization in Atlanta, which might attenuate social distancing values within urbanized areas.
Figure 7 provides the line-path scatter plots of mobility behavior for individual counties in GA: Figure 7a) illustrates county trajectories during the 2020 with the Local Moran's I values calculated in relation to the national average, and Figure 7b) illustrates county trajectories with the Local Moran's I values calculated in relation to the state averages.The blue line-paths inferred from mobility patterns dominate the visual display (Figure 7a)), pointing to the fact that compared to counties nation-wide, the majority of counties in Georgia were classified as coldspots: counties with lower percentage of people clustered together with other counties with lower percentage of people sheltering at home.We can see that the state stay-at-home order (dated 4 April 2020) was more effective in Atlanta-Sandy Springs-Alpharetta metro area, captured by red lines that start to exhibit a hotspot behavior from week 12 across both data sources.Furthermore, there are roughly two groups of coldspots: those that remain consistent (i.e.frequently became a cluster) throughout the pandemic (lines with slope close to 1), and those whose consistency rose in the first 20 weeks of 2020 and then stabilized until week 44 started increasing again as a result of the Presidential Election (3 November 2020) and the Thanksgiving holiday (26 November 2020).Not-surprisingly, the only significant hotspots for both data sources stopped being significant shortly before the Presidential Election (Cuebiq) or sometime between these two dates (SafeGraph).This trend is attenuated at the state-level (Figure 7b).When using the state average, there are fewer counties classified as statistically significant hotspots/coldspots (as quantified by the number of lines and a counter of clusters).Looking at the time aspect of clustering, two groups of counties can be identified: before and after week 25.The former contains the counties that stopped exhibiting clustering behavior around week 25 (June 2020), which could be indicative of the changing mobility patterns (less sheltering, more movement).The other group of counties remained a hotspot/ coldspot at least once during the weeks 25-52.The number of flipping counties in Georgia is relatively low compared to the overall number of counties.

Mobility patterns in South Dakota
The recency and consistency map of mobility in South Dakota is provided in two views: Figure 8a) shows clusters deduced by comparing the county data to the national average, while Figure 8b) captures clusters deduced by comparing the county data to the state average.Compared to Georgia, the influence of the rural-urban gradient is not observed in these mobility patterns at the national level.None of the most populated counties exhibit uniform behavior: Pennington County (home of the Rapid City) exhibits a very inconsistent hotspot behavior, while Minnehaha County (home of Sioux Falls) is a very recent and consistent coldspot.Nonetheless, some geographic and temporal differences are captured: coldspots on the east of the state (east of the Missouri River) and hotspots on the west of the state (west of the Missouri River).The coldspots drawn on the Cuebiq data appear more recent and consistent (markers of larger size with higher color intensity); more counties with lower percentage of people clustering next to other counties with lower percentage of sheltering people through the course of the pandemic.For SafeGraph data, the clusters are less recent with lower degrees of consistency.South Dakota is often partitioned into West and East River regions, with the Missouri River splitting the territory of the state into roughly two halves.The West River is predominantly ranch and dryland farming oriented, while the East river is mostly corn and wheat agricultural production.Furthermore, most of SD mining operations are in the West River.In South Dakota, there are eight counties, where the majority of the population is Native American (Bennett, Buffalo, Corson, Dewey, Mellette, Oglala Lakota, Todd, and Ziebach), seven of which are located in the West River.Both Cuebiq and SafeGraph agree on classification for two counties  There is much more disagreement between the sources for LISA values calculated based on the state average (Figure 8b)): for Cuebiq, two hotspots are very recent and centered in the southwestern corner of the state, with less recent coldspots located along the Missouri River.For SafeGraph, the less recent hotspots are clustered along the western state border and the coldspot counties are located in the middle of the state.The consistency of clusters is noticeably lower (noted in smaller marker size) across local clusters when compared to the clusters deduced from the national averages.The Missouri divide cannot be clearly observed here, and the overall number of counties that were not statistically significant (counties with no circle marker over them) is larger.Looking at both figures and comparing them across two regional scales of analysis, we notice that mobility patterns uncovered in comparison to the national averages do not persist through comparison to the state averages.We also observe a fair amount of inconsistencies between data sources that could be attributed to data quality issues, resulting in representative bias.
The line-path scatter plot provides additional insights into how the counties behaved during the pandemic (Figure 9).First, for mobility patterns calculated at the national level (Figure 9a)), Cuebiq captures consistent coldspot behavior across counties starting from week 20.South Dakota is an example of the state that never implemented a strict state-wide shelter-in-place order (Blackwell, 2021) for the entire population (only the residents of Minnehaha and Lincoln counties who were at least 65 years old or with underlying medical conditions had to stay home).Generally, with a statewide stay-at-home order, we typically expect a hotspot behavior (more people staying home across counties concentrated in a specific location).The comparison to the national average mobility patterns points to the fact that multiple coldspots (fewer people staying at home than on average across the United States).This signal is attenuated for coldspots generated from SafeGraph data: the slope of lines remains flat until week 35 when the lines start rising.For LISA clusters generated from state averages (Figure 9b)), only a few counties exhibit a consistent behavior.Overall, there are no counties that were statistically significant clusters for more than 10 weeks (no values with Count [y-axis] higher than 10).This also points to higher uniformity in individual mobility patterns across the counties and very little predisposition to clustering behavior.Another important observation revealed in Figure 9 is the high number of flipping counties.In fact, 23/66 counties for both Cuebiq and SafeGraph change their clustering behavior during the observation period.Looking closely at the flipper scatter-path (Figure 14), most of the flips are observed for counties that were classified as hotspots within the first 10-20 weeks of 2020 (between January 1 and June 1), followed by some coldspots toward the second half of 2020.

Case study 2: investigating intra-urban mobility in greater Los Angeles in 2020
Los Angeles is a major metropolitan area and is also the second-largest city in the United States with a population of 3.9 million.Furthermore, this area provides a good mix of urban and semi-urban landuse patterns that can influence mobility in different ways, including commercial areas (Glendale, Downtown Los Angeles), tourist areas (Huntington Beach, Santa Monica), logistic centers (Long Beach Port), and others.Since the Mapbox data are highly detailed, only grid cells within the US Census Populated Place Areas for Los Angeles County were considered, thus excluding areas with low population density, forested areas, and over-the-sea grid cells.To simplify the analysis and computations, we generate the recency and consistency map and the line-path scatter plot (Figure 10) by aggregating a 100 m grid (OpenStreetMap zoom level 18) to the "neighborhood" level of 2000 m grid (OpenStreetMap zoom level 14).One major difference between SafeGraph/Cuebiq and Mapbox data, is the tessellation type.Mapbox provides regular gridded data that is continuously distributed over the entire study area including forested areas and ocean tiles, whereas SafeGraph/Cuebiq aggregates mobility at Census designated areal units drawn over land only.This necessitates extra data processing and filtering for gridded Mapbox data due to unequal volume and distribution of mobility over different types of surface (e.g.land/ocean).
The RECO map (Figure 10a)) shows the amount of recency and consistency of local clusters, corresponding to the underlying aggregate mobility values of the entire LA region.The hotspots (shown in red) highlight the tiles with higher activity levels located next to other tiles with high values, whereas coldspots (in blue) correspond to tiles with lower activity levels located next to other tiles with low values.
Figure 10a) highlights the spatial distribution of clusters in Los Angeles.First, we notice two large recent hotspots: in central Los Angeles (Hollywood, Fashion District, Skid Row, Koreatown, Beverly Groves) and in coastal Santa Monica and Venice.More localized hotspots are also noticeable in Van Nuys, Costa Mesa, and Newport Beach, where many local points of interest, such as shopping and recreational activities are located.The coldspots are situated in Long Beach (home to the second largest port in the United States).Since COVID-19 disrupted many world supply chains, the signal we observe here is mostly due to a drastic decrease in cargo shipping.Other coldspots are also localized around Cypress and Santa Ana areas.A series of coldspots to the South-East of the Newport Beach are observed in areas that have a large number of coastal rental properties.Looking at Figure 10b) we see that the hotspots (areas of high mobility) remain consistent throughout 2020, while the coldspots consistency remains low (below 5-10 weeks during 2020).The consistency of hotspots peaks around 30 weeks.The consistency of coldspots is significantly lower, with maximum around 10 weeks.Ideally, we would observe the effect of the stay-at-home order on mobility in the increased number of coldspots (areas of low mobility situated next to other areas of low mobility), but the effect seems to be delayed by 5 weeks.We observe a sharp increase in the amount of statistically significant coldspots around 15 April 2020, when the new social distancing guidelines for the city of Los Angeles were established.Both the RECO map and the line-path scatter plot uncover mobility patterns in relation to the overall average of the study area (i.e. the greater Los Angeles metro area).Therefore, it makes sense that the areas with most activity are concentrated around downtown Los Angeles, where there are many businesses and commercial areas.The concentration of high activity tiles in Newport Beach and in Santa Monica/Venice Beach can be attributed to clustering of hotels, resorts, and other tourist attractions.At this level of analysis, we are unable to differentiate any major transportation arteries: the aggregated spatial grid size is too large to capture the outline of the highways and freeways.The number of flipper tiles is pretty low and is not clustered within a specific geographical area.

Discussion: analytical challenges and future development
In this section, we summarize VASA features and limitations in light of existing literature, and outline some possible avenues for further development.We also provide a few general recommendations on VASA's application with regard to data type and size, geography of features, interactivity, and missingness in the data.

Effective mapping of spatial and spatio-temporal autocorrelation indices
Overall, there are multiple challenges to mapping COVID-19 and its mobility responses effectively (Juergens, 2020), which VASA tries to address.Building on the well-established LISA, VASA is designed for swift exploratory analysis with a particular focus on the spatial and temporal patterns across various levels of analysis.VASA allows assessing the spatial clusters over time across multiple geographic scales, by adding recency and consistency metrics to the available hotspot/coldspot choropleth displays.It further enables the user to track the behavior of individual counties, states, and other political units via the line-path scatter plot.In contrast to differential Local Moran (Anselin, 2019), our visual displays contain more information on the temporal aspects and can capture the change points in time at each spatial unit (for instance, change in county-level responses to imposing/lifting nationwide or statewide stay-at-home orders).In comparison with the space-time scan statistic (Hohl et al., 2020), which allows to effectively monitor and identify COVID-19 clusters, the proposed visualization technique is more appropriate for exploratory data analysis and does not require the user to specify the size of the spatial and temporal windows.While space-time cube may be used as a viable alternative for mapping raw or gridded dynamic series (Mo et al., 2020;Purwanto et al., 2021), the visual display becomes easily cluttered when the areal units (e.g.county, state) are rendered in 3D over time.
While there are multiple implementations of combined spatio-temporal autocorrelation (Gao et al., 2019;Griffith & Heuvelink, 2012;Lee & Li, 2017), we opted to use purely spatial Local Moran's I for three reasons.First, Local Moran's I is easier to interpret and is more commonly used across various disciplines, including regional science, spatial econometrics, economics, spatial statistics, and geography.Second, it provides an easy mapping from I to four clustering types (which are more easily interpreted than continuous unit-less I values).Third, for combined spatiotemporal indices, it might be harder to separate the individual effect of spatial and temporal variations.By coupling the spatial Local Moran's I with temporal representation, we were able to demonstrate temporal mobility trends in relation to non-pharmaceutical interventions.

Specification of spatial weight matrices
The configuration of local spatial clusters depends on the specification of spatial weights matrices.VASA provides most commonly used weight construction methods, including kNN, Queens, and Rooks (and their combination).We generated a series of visualizations to investigate how different specifications affect the overall spatial patterns in the data using recency and consistency maps (see Figure 5, 12, 13 in Appendix).Overall, for national-level county data, we observe very similar cluster configurations, with the highest number of statistically significant clusters for Queens weights, second highest for Rooks, and third largest for kNN.The observed spatial distribution of clusters appears very similar across all three figures.

Dealing with big data
While the case studies outlined in this paper provide a good combination of possible analytical scenarios, where VASA can be used, the size of the problem and, specifically, the number of units under analysis should be carefully managed to avoid information overload.With larger data sets, the computational time is likely to increase and make the visual displays more cluttered.Consider a national map of the United States generated for 3,200 counties (Figure 5).The map provides an easyto-read visual display to ascertain major mobility patterns: Pacific-West and the Northeast are hotspots (more people sheltering at home), Illinois is divided along the urban-rural gradient, with Chicago exhibiting a hotspot behavior and the south of the state is a coldspot.At the same time, looking specifically at the Northeast is troublesome in terms of graphic representation, due to high concentration of counties and high consistency coldspot behavior, which results in the clumping of larger markers.The RECO map can be used to chart a few thousand features.However, the density of features matters, particularly for Census units that are attached to the population (i.e.Census block groups).Therefore, depending on the context, it may be required to aggregate the data one level up (i.e.Census tracts) and tweak the size of centroid markers to remove occlusion from the map (see Figure 5).
For the generated line-path scatter plots above, the number of visualized features was 159 in Georgia and 66 in South Dakota.As the number of features (lines) increases, the visual display may become more cluttered.For cases with over 1000 observations on the plot, we recommend utilizing VASA's "sampling" functionality that makes such visual displays easier to read and differentiate between clusters.Thus, sampling sufficient amount of observations could help assess the temporal structure of clustering behavior, which is the primary purpose of the line-path scatter plot.At the same time, it is also important to keep the duration of observation of the study period relatively high (e.g.above 10).For example, with only 5 days/weeks of observations the amount of variation within a ,5000 geographic features might be very low, rendering the line-path scatter plot useless.For such cases, we recommend the users to use "jitter" and "sampling" functionality.In future developments, we plan to add bundling to the lines, to highlight the trend better in such situations.

Handling missing data
Missing data can have major impacts on the assessment.For the proposed analytical pipeline, the data is assessed sequentially, with Local Moran's I calculated for each time step.Therefore, a missing value may result in unconnected/broken spatial matrix.VASA provides functionality to assess missingness or gaps in the provided data and different options for constructing spatial weights matrices (e.g.contiguity, kNN, distance-bands, and others).Normally, there are two ways of dealing with such scenarios: imputing missing values and adjusting the type of neighborhood (i.e. using distance threshold instead of Queen neighborhood).Depending on the type of missingness (missing at random, missing completely at random, missing not at random), the imputation methods should be tuned appropriately for each data set (Rubin, 1976).Since the gaps in spatio-temporal data are rarely random, the spatial and temporal neighborhood structure should be considered.At this beta-stage of development, VASA utilizes commonly used time-series tools for such imputation: the global average and the rolling average.We recommend that the spatial units that are missing values for the substantial amount of time (i.e.several weeks) are dropped from the analysis.Otherwise, the imputation techniques cannot provide a good enough approximation.In the future, more advanced spatiotemporal interpolation methods can be added to the package.

Handling areal data and tessellations
Evidently, the analysis on the regular Mapbox grid posed multiple challenges as compared to the polygonal data used here.These challenges are mostly related to the format and specifications of the Mapbox data and not VASA's visual displays.First, because the square grid does not extend naturally to the spatial features of the urban forms (as in Census Blocks), the tool may yield coldspots at locations, where there is no or very little movement, for instance, in forested areas and open waters.Additionally, filtering out the "no-movement" tiles generates additional difficulties for constructing the adjacency matrices (especially for contiguity defined spatial weights matrices -e.g.Rooks, Queens) as many tiles may be removed rendering a disjoint study area.One possible solution would be to filter the tiles that are specifically located within city limits, where there is more movement.Unfortunately, this would skew the focus of the studies toward investigating business areas and commercial districts due to the way the data is collected and aggregated by Mapbox, which effectively limits the signal in residential areas via privacy algorithms.
The fine spatial resolution for Mapbox poses two major challenges for visual analytics: 1) RAM limitations for in-memory processing, and 2) occluded visual display due to the large size of the study area.To make computations more manageable, aggregation is required.At this stage of development, VASA does not contain automatic procedures to identify an ideal sampling scale.As a rule of thumb, we recommend using OpenStreetMap tile zoom levels: town (12), village (13), small road (15), street (16), block/park (17), some buildings (18).For analysis in this paper, we selected the zoom level 14, which is situated between "village" and "small road" and is indicative of the neighborhoods.Additionally, aggregation based on administrative units can be considered (e.g.Census block groups).Future extensions may consider applying a multi-scale rendering where the aggregation level changes on the fly based on the selected zoom level.

Adding interactivity
The developed visualizations can be made interactive and web-based, enhancing user experience and VASA exploration capabilities.Examples of such visualizations have been developed and are available on the Lab website 3 .In future extensions of VASA, we also plan to add an interactive version of the RECO and the linepath scatter plot.Specifically, by utilizing Altair API, we plan to make the following adjustments: 1) enable the zooming and panning over the map canvas; 2) interactive selection of hotspots/coldspots from the legend entries; 3) adding a dropdown selection widget to look at states (or other units of analysis) individually; 4) using a slider widget to investigate temporal trends; 5) adding linking and brushing to couple RECO and linepath scatter plots; 6) Adding pipeline to switch between spatial and temporal units of aggregation.

Animated movement
Animation plays an important role for cartographic visualization of spatial patterns in time (DiBiase et al., 1992).When prototyping VASA, we implemented various animated views to demonstrate how the spatio-temporal patterns of mobility changed over 2020.These visualizations are provided on the https://loquacious-frangipane-b89dda.netlify.appwebsite.The animated choropleth maps are able to show weekly county behavior over time.The weekly LISA classification of counties is shown by hotspot/coldspots being colored red/blue and non-significant counties not being colored.The animation cycles through each week, updating the maps as it plays.The user has the ability to pause the animation, skip forward or backward a week, or drag the timeline cursor to a specific week.There is an additional option to make the animation cumulative, layering classifications of previous weeks together with the current week, while also fading away distant classifications as the animation plays forward.There is some interactivity to these maps as well.Clicking on a state will center a view on that state.This also opens a table that lists all of the counties in the state and displays the counties' raw mobility indicators values used and the derived LISA classifications.

Conclusion
This paper introduced an analytical framework for exploratory analysis on aggregate movement data along with the corresponding package VASA.We demonstrated how mobility data can be effectively mapped to infer spatial and temporal structure of patterns in response to the COVID-19 pandemic and relevant non-pharmaceutical interventions.Using examples with regular and irregular tessellation data, we showed how VASA can be adapted for ESDA across scales: for national-, state-, and regional-level analysis using multiple data sets such as SafeGraph, Cuebiq, and Mapbox mobility indices.VASA can be an important and effective tool in a spatial analyst toolbox, helping guide the knowledge creation and hypothesis generation process at the start of data analysis pipelines.By looking into three different cases, we were able to describe the evolution of spatial clusters of mobility during the COVID-19 pandemic in different areas.The preliminary findings from these research cases can be used to guide further research: 1) the differences between urban and rural counties in Georgia; 2) the influence of political affiliation on movement patterns; 3) the mobility patterns in counties with predominantly American Indian populations.Overall, VASA demonstrated to be more effective for largescale areal (county-level) data.The assessment of the regular finer resolution gridded data poses a number of challenges for spatial analysts: the influence of scale and modifiable areal unit problem might affect the configuration of spatial clusters in an unpredictable manner.Further research is needed to investigate these effects.Another major thread of future research will focus on assessing the effectiveness of the proposed visual displays and visual encodings in typical user tasks for exploratory data analysis.
• SafeGraph -free of charge under the "Data for Good" initiative for academic purposes.• Cuebiq -proprietary data, which requires licensing fee for access.
Restrictions apply to the availability of these data, which were used under license for this study.

Figure 1 .
Figure 1.VASA analytical framework.The inputs are denoted in yellow quadrilaterals, package objects and classes are denoted in purple hexagons, intermediate processing methods are illustrated using blue rectangles, and visual displays are represented using arrow-like green boxes.The modules are denoted in dashed rectangles.

Figure 2 .
Figure 2. a) Centroids are generated and mapped as circle markers; b) LISA is run to generate hotspots/coldspots, which are mapped to red/blue colors; c) Recency is mapped to color value; d) Consistency (number of times a county (or other object-level unit) is a cluster) is mapped to a marker size; e) Recency and consistency are finalized in a RECO map.

Figure 3 .
Figure 3. a) Cumulative hotspot/coldspot line segments are generated and connected; b) Line endpoints for the most recent hotspot/ coldspot behavior are added to lines; c) Line-paths for hotspot/coldspots are colored in red/blue; d) Line endpoints for hotspot/ coldspots are colored red/blue, counties with both hotspots and coldspots (i.e.flippers) are not filled.The tally of hotspots, coldspots and flips is added to the plot.

Figure 4 .
Figure 4. a) Cumulative hotspot/coldspot line segments are generated and connected; b) Line endpoints for the most recent hotspot/ coldspot behavior are added to lines; c) Line segments are classified as hotspot/coldspots depending on the county's most recent significant value and colored red/blue, endpoints are colored by the county's last week behavior d) the week the flipped behavior occurred is marked by an X.

Figure 5 .
Figure 5. Recency and consistency map for 3,200 US counties over the 2020 COVID-19 pandemic calculated for cuebiq and SafeGraph %-sheltered in place metric.

Figure 6 .
Figure 6.Recency and consistency mapping for aggregate mobility in Georgia.a) Comparing mobility patterns in GA counties to national average; b) Comparing mobility patterns in GA counties to the state average.The counties that intersect with the Atlanta metro area are denoted in beige.
only: Buffalo is a coldspot, while Bennett is a hotspot.Corson County is classified as a recent hotspot by Cuebiq and as a non-recent coldspot by SafeGraph.Further research is necessary to establish what caused such inconsistencies in the classification across the two data sources.It is possible that the coverage and representativeness of the sample is at fault.

Figure 7 .
Figure 7. Line-path scatter plot for % sheltered in place in Georgia.a) Comparing mobility patterns in GA to the national average; b) Comparing mobility patterns in GA counties to the state average.

Figure 8 .
Figure 8. Recency and consistency mapping for aggregate mobility in South Dakota.a) Comparing mobility patterns in SD counties to the national average; b) Comparing mobility patterns in SD counties to the state average.The Missouri River is denoted in blue, the counties where majority of population is Native American is denoted in beige.

Figure 9 .
Figure 9. Line-path scatter plot for % sheltered in place in South Dakota.a) Comparing mobility patterns in SD counties to the national average; b) Comparing mobility patterns in SD counties to the state average.

Figure 10 .
Figure 10.Mobility patterns captured in mapbox activity index data: a) the recency and consistency map; and b) the line-path scatter plot.