New tools for visualising and explaining multivariate spatio-temporal data

H. Sherry Zhang, PhD stduent

Monash University, Australia

Apr 17, 2023

Hi!

• A final year PhD student in the Department of Econometrics and Business Statistics, Monash University, Australia

• My research centers on exploring multivariate spatio-temporal data with data wrangling and visualisation tool.

• Find me on

  • Twitter: huizezhangsh,
  • GitHub: huizezhang-sherry, and
  • https://huizezhangsh.netlify.app/

• slides for today: https://sherryzhang-ireland2023.netlify.app/

• Thanks for the invitation to speak

• Today I will be talking about …

• First, a little bit about myself

• I’m Sherry Zhang, …

• …

• Here are the details to find me on Twitter, GitHub, and my website

Multivariate spatio-temporal data can be commonly found in the society: in finance and economics, stock prices and economic indicators are tracked over time; in logistics, supermarkets collect product level data to decide the optimal stock levels for different stores; in meteorology, weather stations record climate variables to monitor climate change and its impact on different sectors.

Spatio-temporal data are recorded with data-time and geographic location information. Multivariate means that multiple variables are recorded. This could include, for example, in weather station data, temperature, precipitation, wind speed and direction

However, all three aspects are ideally considered together to tackle contemporary problems, such as monitoring droughts which requires historical climate data to understand “normal” conditions for any spatial neighborhood, and the interplay of precipitation, temperature and other relevant variable(s), for example, ice melt for high altitude mountain region.

This research addresses the challenge of investigating multivariate spatio-temporal data, by providing new tools for organising, visualizing and explaining relationships.

This illustration here shows how the three research topics are related. The solutions are to provide easy ways to pivot between the three components, to allow focusing on multivariate or spatio-temporal analysis and new pipelines to construct indexes.

When fixing the time, the data are reduced to multivariate and can be analyzed using multivariate methods such as dimension reduction. The particular dimension reduction investigated in this thesis is called projection pursuit, which contains an optimization component. Diagnostic plots can be useful here to track the performance of these optimization algorithms.

When the data are collected at different locations in space, software from geoinformatics can be useful to analyze the spatial aspect of the data. However, existing spatial and temporal data analysis software are built upon different data formats. This creates frictions in the data analysis to constantly rearrange that data format to work with the software. We designed a spatio-temporal data structure to organise the data.

Multivariate spatio-temporal variables are often combined into a single index series for decision making. But index definition and construction are vastly different in different fields and different researchers making it difficult to understand how they might perform slight changes in the formula, or how competing indexes compare.

Visual diagnostics for constrained optimisation with application to guided tours

Projection pursuit for dimension reduction

I will first use a simple example to explain the projection pursuit algorithm, before showing two visualisation on diagnosing the optimisers.

On the left is a density plot with the marginal distribution of the two variables on the top and right. We can construct a projection matrix and multiple it with the data to get a projection.

On the right are 21 projection matrices and their corresponding projections. If we give X1 zero weight and X2 full weight, you will get the distribution of X2 as in the first panel on the right. Similarly, ….

Sometimes, the interesting structure in the high-dimensional data may be unknown and we could use a grand tour to randomly generate projection matrices and view the projections. Once we have some ideas of the structure we’re looking for, for example, clusters or central mass, we will choose an index function, which maps the projection to a single value. The guided tour uses optimisation to maximum the index function given the projection matrices.

Optimisation in projection pursuit

Data: \(\mathbf{X}_{n \times p}\);

Projection matrix: \(\mathbf{A}_{p\times d}\)

Projection: \(\mathbf{Y}_{n \times d} = \mathbf{X} \cdot \mathbf{A}\)

Index function \(f: \mathbb{R}^{n \times d} \mapsto \mathbb{R}\)

Optimisation: \[\arg \max_{\mathbf{A}} f(\mathbf{X} \cdot \mathbf{A}) ~~~ s.t. ~~~ \mathbf{A}^{\prime} \mathbf{A} = I_d\]

We can put these into notations. The data \(X\) is n * p, N being the number of observation and P being the original data dimension. The projection matrix A, has dimension of p * d. The projection that multiply X and A will be the data in the lower dimension d.

The index function maps the projection into a single scale, which will be optimized over the space of A. Also the projection matrices have to subject to the orthonormality condition to be a projection matrix.

We use stochastic optimisers for this problem. These optimisers repetitively generate projection matrices, decide if to accept or reject the new position, to find the maximum. This optimisation problem differs from a common problem in that the input, A, is a matrix that contains p * d numbers. It is not easy to understand where in the matrix space, the optimisers have visited, so we would like some visual tools to help here.

Visualise the projection matrix space

Data: \(\mathbf{X}_{n \times p}\);

Projection matrix: \(\mathbf{A}_{p\times d}\)

Projection: \(\mathbf{Y}_{n \times d} = \mathbf{X} \cdot \mathbf{A}\)

Index function \(f: \mathbb{R}^{n \times d} \mapsto \mathbb{R}\)

Optimisation: \[\arg \max_{\mathbf{A}} f(\mathbf{X} \cdot \mathbf{A}) ~~~ s.t. ~~~ \mathbf{A}^{\prime} \mathbf{A} = I_d\]

Simulation:

• simulated 5D data project to 1D
• two optimisers

In this example, we simulated 5D data and project them into 1D.

Because the data is simulated, we know the projection matrix that gives the most interesting projection, which is the star here.

Viewing a 5D space can be challenging, so here, we perform a PCA on all the bases to obtain a 2D view.

Because of the orthonormality condition, the domain for 1D projection is a p-d unit sphere and it becomes a circle when further reduced by PCA.

Here you can see how the two optimisers progress from the start to find the maximum.

Visualise the projection matrix space

Visualise the projection matrix space

Cubble: An R package for organizing and wrangling multivariate spatio-temporal data in R

Motivation

If you have worked with spatial data in R, you will probably know about the sf package. In an oversimplified term, the package uses nested list to wrap up coordinates into a special column, geom. The package also connects to external libraries that allow us to do geospatial operations.

On the other hand, we have tsibble that arranges time series data into a long form with explicit column for the time variable.

Australian weather station data

stations

# A tibble: 88 × 6
   id            lat  long  elev name              wmo_id
   <chr>       <dbl> <dbl> <dbl> <chr>              <dbl>
 1 ASN00001006 -15.5  128.   3.8 wyndham aero       95214
 2 ASN00002032 -17.0  128. 203   warmun             94213
 3 ASN00003080 -17.6  124.  77.5 curtin aero        94204
 4 ASN00005007 -22.2  114.   5   learmonth airport  94302
 5 ASN00006044 -25.9  114.   9   denham             94402
 6 ASN00007600 -28.1  118. 407   mount magnet aero  94429
 7 ASN00008296 -29.2  116. 271.  morawa airport     94417
 8 ASN00009114 -31.0  115.   4   lancelin           95606
 9 ASN00009240 -32.0  116. 384   bickley            95610
10 ASN00009542 -33.7  122. 142   esperance aero     95638
# … with 78 more rows

ts

# A tibble: 32,208 × 5
   id          date        prcp  tmax  tmin
   <chr>       <date>     <dbl> <dbl> <dbl>
 1 ASN00001006 2020-01-01   164  38.3  25.3
 2 ASN00001006 2020-01-02     0  40.6  30.5
 3 ASN00001006 2020-01-03    16  39.7  27.2
 4 ASN00001006 2020-01-04     0  38.2  27.3
 5 ASN00001006 2020-01-05     2  39.3  26.7
 6 ASN00001006 2020-01-06    60  32.9  25.6
 7 ASN00001006 2020-01-07   146  34.1  25.5
 8 ASN00001006 2020-01-08    40  36.6  26.2
 9 ASN00001006 2020-01-09     0  38.2  27.6
10 ASN00001006 2020-01-10     0  38.9  29.7
# … with 32,198 more rows

Initially I was working with climate data where I have 30 weather stations and its spatial information and time series data to record the daily climate variables in 2020. I was trying to use sf to handle the spatial side and tsibble to handle the temporal side. But the two packages do not work well with each other. This is expected since neither sf or tsibble is designed to handle spatio-temporal data.

This motivates me to think about the spatio-temporal data structure in R.

What’s available for spatio-temporal data? - stars

What’s available at that time is a package called stars, it uses a dense array to structure spatio-temporal data. You can think of it as stacking snapshots of the space along the time axis.

This is great for satellite data, but it may not be the most obvious solution for analysts who prefer to operate on a 2D table format.

Hence, I designed a data structure called cubble to handle saptio-temporal vector data.

Cubble: a spatio-temporal vector data structure

Cubble is a nested object built on tibble that allow easy pivoting between the spatial and temporal form.

The nested form is similar to the sf data frame you seen before, with an additional list column called ts that nests all the temporal variables

The long form mimics the long table in tsibble where each row is cross identified by the site and date in a long table

Cubble: a spatio-temporal vector data structure

Cubble is a nested object built on tibble that allow easy pivoting between spatial and temporal form.

What cubble do is to provide a linking between the two forms through the pair face_temporal() and face_spatial().

With face_temporal(), the focus of the data is now on the temporal face of the spatio-temporal cube and this corresponds to switch the data to the long form.

With face_spatial(), the long cubble is switched back to the nested form, the spatial face of the datacube.

Cast your data into a cubble

(weather <- as_cubble(
  list(spatial = stations, temporal = ts),
  key = id, index = date, coords = c(long, lat)
))

# cubble:   id [88]: nested form
# bbox:     [113.53, -43.49, 153.64, -10.58]
# temporal: date [date], prcp [dbl], tmax [dbl], tmin [dbl]
   id            lat  long  elev name              wmo_id ts                
   <chr>       <dbl> <dbl> <dbl> <chr>              <dbl> <list>            
 1 ASN00001006 -15.5  128.   3.8 wyndham aero       95214 <tibble [366 × 4]>
 2 ASN00002032 -17.0  128. 203   warmun             94213 <tibble [366 × 4]>
 3 ASN00003080 -17.6  124.  77.5 curtin aero        94204 <tibble [366 × 4]>
 4 ASN00005007 -22.2  114.   5   learmonth airport  94302 <tibble [366 × 4]>
 5 ASN00006044 -25.9  114.   9   denham             94402 <tibble [366 × 4]>
 6 ASN00007600 -28.1  118. 407   mount magnet aero  94429 <tibble [366 × 4]>
 7 ASN00008296 -29.2  116. 271.  morawa airport     94417 <tibble [366 × 4]>
 8 ASN00009114 -31.0  115.   4   lancelin           95606 <tibble [366 × 4]>
 9 ASN00009240 -32.0  116. 384   bickley            95610 <tibble [366 × 4]>
10 ASN00009542 -33.7  122. 142   esperance aero     95638 <tibble [366 × 4]>
# … with 78 more rows

• the spatial data (stations) can be an sf object and temporal data (ts) can be a tsibble object.

• To cast the two separate tables into a cubble, you can supply them in a named list.

• You also need to tell cubble some identifiers it looks for

• The key argument is the spatial identifier that connects the two tables.

• The index argument is the temporal identifier that prescribes the timestamp.

• The coords argument is to used to specify the coordinate

[breath]

• From the cubble header, you can read that the key is id, there are 30 stations and it is in the nested form.

• The third line here shows you the available temporal variables and their types.

• Also, if the spatial and temporal data is an sf or tsibble object, they will be indicated in the header as well.

Switch between the two forms

long form

(weather_long <- weather %>% 
  face_temporal())

# cubble:  date, id [88]: long form
# bbox:    [113.53, -43.49, 153.64, -10.58]
# spatial: lat [dbl], long [dbl], elev [dbl],
#   name [chr], wmo_id [dbl]
   id          date        prcp  tmax  tmin
   <chr>       <date>     <dbl> <dbl> <dbl>
 1 ASN00001006 2020-01-01   164  38.3  25.3
 2 ASN00001006 2020-01-02     0  40.6  30.5
 3 ASN00001006 2020-01-03    16  39.7  27.2
 4 ASN00001006 2020-01-04     0  38.2  27.3
 5 ASN00001006 2020-01-05     2  39.3  26.7
 6 ASN00001006 2020-01-06    60  32.9  25.6
 7 ASN00001006 2020-01-07   146  34.1  25.5
 8 ASN00001006 2020-01-08    40  36.6  26.2
 9 ASN00001006 2020-01-09     0  38.2  27.6
10 ASN00001006 2020-01-10     0  38.9  29.7
# … with 32,198 more rows

back to the nested form:

(weather_back <- weather_long %>% 
   face_spatial())

# cubble:   id [88]: nested form
# bbox:     [113.53, -43.49, 153.64, -10.58]
# temporal: date [date], prcp [dbl], tmax [dbl],
#   tmin [dbl]
   id        lat  long  elev name  wmo_id ts      
   <chr>   <dbl> <dbl> <dbl> <chr>  <dbl> <list>  
 1 ASN000… -15.5  128.   3.8 wynd…  95214 <tibble>
 2 ASN000… -17.0  128. 203   warm…  94213 <tibble>
 3 ASN000… -17.6  124.  77.5 curt…  94204 <tibble>
 4 ASN000… -22.2  114.   5   lear…  94302 <tibble>
 5 ASN000… -25.9  114.   9   denh…  94402 <tibble>
 6 ASN000… -28.1  118. 407   moun…  94429 <tibble>
 7 ASN000… -29.2  116. 271.  mora…  94417 <tibble>
 8 ASN000… -31.0  115.   4   lanc…  95606 <tibble>
 9 ASN000… -32.0  116. 384   bick…  95610 <tibble>
10 ASN000… -33.7  122. 142   espe…  95638 <tibble>
# … with 78 more rows

identical(weather_back, weather)

[1] TRUE

• Here is what a cubble look like when being switched between the long and the nested form.

  • With the weather object we just created, we turn it into the long form with the function face_temporal()

• Notice that the third line in the header now changes to see the available spatial variables

[breath]

• On the right, weather_long is switched back the nested form with the function face_spatial()

• As you can see from the last line of code, face_temporal() and face_spatial() are the exact inverse.

• Hence weather_back and weather are identical

Access variables in the other form

Reference temporal variables with $

weather %>% 
  mutate(avg_tmax = mean(ts$tmax, na.rm = TRUE))

# cubble:   id [88]: nested form
# bbox:     [113.53, -43.49, 153.64, -10.58]
# temporal: date [date], prcp [dbl], tmax [dbl], tmin [dbl]
   id            lat  long  elev name              wmo_id ts                 avg_tmax
   <chr>       <dbl> <dbl> <dbl> <chr>              <dbl> <list>                <dbl>
 1 ASN00001006 -15.5  128.   3.8 wyndham aero       95214 <tibble [366 × 4]>     36.7
 2 ASN00002032 -17.0  128. 203   warmun             94213 <tibble [366 × 4]>     35.8
 3 ASN00003080 -17.6  124.  77.5 curtin aero        94204 <tibble [366 × 4]>     35.9
 4 ASN00005007 -22.2  114.   5   learmonth airport  94302 <tibble [366 × 4]>     33.2
 5 ASN00006044 -25.9  114.   9   denham             94402 <tibble [366 × 4]>     27.1
 6 ASN00007600 -28.1  118. 407   mount magnet aero  94429 <tibble [366 × 4]>     30.3
 7 ASN00008296 -29.2  116. 271.  morawa airport     94417 <tibble [366 × 4]>     28.8
 8 ASN00009114 -31.0  115.   4   lancelin           95606 <tibble [366 × 4]>     24.8
 9 ASN00009240 -32.0  116. 384   bickley            95610 <tibble [366 × 4]>     22.8
10 ASN00009542 -33.7  122. 142   esperance aero     95638 <tibble [366 × 4]>     22.8
# … with 78 more rows

Move spatial variables into the long form

weather_long %>% unfold(long, lat)

# cubble:  date, id [88]: long form
# bbox:    [113.53, -43.49, 153.64, -10.58]
# spatial: lat [dbl], long [dbl], elev [dbl], name [chr], wmo_id [dbl]
   id          date        prcp  tmax  tmin  long   lat
   <chr>       <date>     <dbl> <dbl> <dbl> <dbl> <dbl>
 1 ASN00001006 2020-01-01   164  38.3  25.3  128. -15.5
 2 ASN00001006 2020-01-02     0  40.6  30.5  128. -15.5
 3 ASN00001006 2020-01-03    16  39.7  27.2  128. -15.5
 4 ASN00001006 2020-01-04     0  38.2  27.3  128. -15.5
 5 ASN00001006 2020-01-05     2  39.3  26.7  128. -15.5
 6 ASN00001006 2020-01-06    60  32.9  25.6  128. -15.5
 7 ASN00001006 2020-01-07   146  34.1  25.5  128. -15.5
 8 ASN00001006 2020-01-08    40  36.6  26.2  128. -15.5
 9 ASN00001006 2020-01-09     0  38.2  27.6  128. -15.5
10 ASN00001006 2020-01-10     0  38.9  29.7  128. -15.5
# … with 32,198 more rows

• Sometimes, you may need to access variables from the other form for your analysis.

• For example, we may want to calculate some per station summary of the time series data.

• We can refer to the temporal variables from the nested form with the $ sign.

• Here I’m calculating the average maximum temperature across the whole year for each station and I need to get access to tmax from the list-column ts.

• In the long form, you need the cubble verb unfold() to move the spatial variables into the long form.

• Here I move the two coordinate columns into the long form and later we will see how it can help us to create a glyph map.

Explore temporal pattern across space with a glyph map

Sometimes, although we technically have spatio-temporal data, we may choose to fix a few stations to explore their temporal patterns, or select a few timestamps to explore their spatial distributions.

A more holistic approach is to explore the space and time simultaneously and now we will see an example of how to use cubble to explore temporal pattern across space in a glyph map.

Why do you need a glyph map?

Here is a typical plot you may have seen when someone tries to visualise their spatio-temporal data. The x and y axes are the coordinates, here I simplify it with only two points, but in reality you may see a collection of points in space or a raster image. Each facet here shows the space in different timestamp and the values are mapped into color.

The problem of this type of visualisation is that it becomes difficult to comprehend the temporal structure of the data since you have to compare points at the same location across panels to digest the pattern.

Why do you need a glyph map?

Instead the temporal pattern is much easier to observe if shown in a time series plot.

What a glyph map do is to put the time series glyph in the place of the location, so you can see the temporal trend in the space.

Glyph map transformation

DATA %>%
  ggplot() +
  geom_glyph(
    aes(x_major = X_MAJOR, x_minor = X_MINOR,
        y_major = Y_MAJOR, y_minor = Y_MINOR)) +
  ...

• I have a short illustration to show you how the transformation works

• Here (1) shows a single station on the map with its long and lat coordinate and (2) is its associated time series.

• Here you know the range of your x and y axis and you can use linear algebra to transform them into a different scale.

• In step (3), the time series in still the same but its scale has been transformed to a width of 1 and heights of 0.3 and the center in this scale is where the original point lays.

• Once we have the time series in the transformed axes, they can be placed onto the map as in (4)

• To make a glyph map, you can use the geom_glyph function from the cubble package.

• It requires a pair of major and a pair of minor variable as required aesthetics

• The major variable are the spatial coordinates, long and lat here and the minor variable are the temporal coordinates, date and tmax here.

Avg. max. temperature on the map

cb <- as_cubble(
  list(spatial = stations, temporal = ts),
  key = id, index = date, coords = c(long, lat)
)

cb_glyph <- cb %>%
  face_temporal() %>%
  group_by(month = lubridate::month(date)) %>% 
  summarise(tmax = mean(tmax, na.rm = TRUE)) %>% 
  unfold(long, lat)

cb_glyph %>% 
  ggplot(aes(x_major = long, x_minor = month,
             y_major = lat, y_minor = tmax)) +
  geom_sf(data = oz_simp, fill = "grey90", 
          color = "white", inherit.aes = FALSE) +
  geom_glyph_box(width = 1.3, height = 0.5) + 
  geom_glyph(width = 1.3, height = 0.5) + 
  ggthemes::theme_map()

Here is an example of using cubble

First you create a cubble object with the stations and ts data using as_cubble. This by default creates the cubble in the nested form, easy for spatial operations.

The second block involves wrangling the data using the nested and long form. Sampling 20 stations is a spatial operations, so it is performed in the nested form. Then we need to do a summary of average maximum temperature by month. It is a temporal operation, so the cubble is then switched to the long form with face_temporal().

unfold is another cubble verb that moves the spatial variable into the long form temporarily. We need to do this because ggplot requires all the variables to be in the same table

geom_glyph() is also implemented in cubble, it is intended to make a glyph map.

[breath]

• On the map, you can see that the temperature curve in the north and south (the Tasmania Island) are relative constant throughout the year.

• Those inland stations, for example in the eastern Australia, have a much visible variation in the year, as compared to the coastline ones.

Additional Information

Slides created via quarto available at

https://sherryzhang-ireland2023.netlify.app/

All the materials used to prepare the slides are available at

https://github.com/huizezhang-sherry/ireland2023

cubble package: https://CRAN.R-project.org/package=cubble

ferrn package: https://CRAN.R-project.org/package=ferrn

ferrn paper: https://doi.org/10.32614/RJ-2021-105

H. Sherry Zhang

Supervised by Dianne Cook, Patricia Menéndez, Ursula Laa, and Nicolas Langrené