Basic Estimation

Multitapering is a popular technique for estimating the spectral density function of a signal in time[1]. Multitapering can also be extended to spatial processes, including random fields [2], point processes [3], marked point processes [4] and multivariate processes which are a mixture of the three [4].

This tutorial demonstrates the basic usage of the SpatialMultitaper.jl package for spectral estimation from spatial data. We will cover how to prepare your data, perform spectral estimation, and visualize the results.

using SpatialMultitaper, GeoStatsProcesses

import GLMakie as Mke

Preparing Your Data

First, we need to load and prepare our spatial data. Data should be passed into multitapering functions as a either PointSet or GeoTable objects. If the data is multivariate, it should be a Tuple of such objects. If the data is a GeoTable, and the domain of this data is a PointSet, then this is interpreted as a marked point process, where the first column of the references information is the mark. If the domain is a CartesianGrid, then this is interpreted as a random field sampled on that grid. Again we take only the first column, so you should pass multivariate processes as separate GeoTable objects, even if they are recorded at the same location.

As a simple example, say we have three independent point processes observed on the same region:

region = Box(Point(0, 0), Point(100, 100))
X = rand(PoissonProcess(0.01), region)
Y = rand(PoissonProcess(0.01), region)
data = spatial_data((X, Y), region)

Spatial data in 2 dimensions
  processes: 1:2
  region: Box(min: (x: 0.0 m, y: 0.0 m), max: (x: 100.0 m, y: 100.0 m))
  with values of type Tuple{PointSet{𝔼{2}, CoordRefSystems.Cartesian2D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, PointSet{𝔼{2}, CoordRefSystems.Cartesian2D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}

constructing the SpatialData object automatically restricts the data to the specified region.

We can visualise this as follows:

viz(boundary(region), color = :gray, axis = (aspect = 1,))
viz!(X, color = :black)
viz!(Y, color = :red)
Mke.current_figure()

Estimation

We can perform spectral estimation using the spectra function. This function takes the data object we already created as an argument. In addition, we need to specify the tapers to use, and the maximum wavenumber in each dimension kmax. Optionally, we can also specify the number of wavenumbers to evaluate in each dimension nk, or the spacing in wavenumber dk. If neither nk nor dk are specified, then dk=1/L where L is the length of the bounding box region in that dimension, and nk is the number of wavenumbers that would be computed given dk.

tapers = sin_taper_family((4, 4), region)
nk = (100, 100)
kmax = (0.1, 0.1)
spec = spectra(data; tapers = tapers, nk = nk, kmax = kmax)

spectra of a 2 dimensional process
  between processes: 1:2 and 1:2
  evaluated at -0.1:0.002:0.098, -0.1:0.002:0.098
  with values of type Matrix{StaticArraysCore.SMatrix{2, 2, ComplexF64, 4}}

Visualising the output

The spectral estimate is returned as a Spectra object. There are various transformations we can apply to this object. The object is multidimensional, but we can index it to get the estimate between two processes:

spec11 = spec[1, 1]

spectra of a 2 dimensional process
  between processes: 1 and 1
  evaluated at -0.1:0.002:0.098, -0.1:0.002:0.098
  with values of type Matrix{ComplexF64}

this returns the spectrum of the first process with itself, i.e. the power spectrum of the first process. Similarly, spec[1, 2] would return the cross-spectrum between the first and second processes.

Certain marginal transformations are available, such as taking the real part:

r_spec11 = real(spec11)

real(spectra) of a 2 dimensional process
  between processes: 1 and 1
  evaluated at -0.1:0.002:0.098, -0.1:0.002:0.098
  with values of type Matrix{Float64}

Then to plot we can call the appropriate plotting function for the dimensionality of the statistic. In this case, the spectrum is two-dimensional, so we can use a heatmap:

Mke.heatmap(r_spec11)

A more interesting example

Let's consider a more interesting example, where we have a bivariate process with some cross-correlation. We can simulate a very simple process by

region = Box(Point(0, 0), Point(100, 100))
shift = 10.0
big_region = Box(Point(-shift, -shift), Point(100 + shift, 100 + shift))
X = rand(PoissonProcess(0.01), big_region)
Y = Translate(shift, shift)(X)
data = spatial_data((X, Y), region) # automatically restrict to region

tapers = sin_taper_family((4, 4), region)
nk = (100, 100)
kmax = (0.1, 0.1)
spec = spectra(data; tapers = tapers, nk = nk, kmax = kmax)

Mke.heatmap(real(spec[1, 2]))

However, the cross-spectrum is complex-valued, and so it is often more useful to look at the magnitude coherence and phase. These are defined as $|C_{XY}(k)| = {|f_{XY}(k)|}/\sqrt{f_{XX}(k) f_{YY}(k)}$ and $\phi_{XY}(k) = \arg(f_{XY}(k))$ respectively, where $f_{XY}(k)$ is the cross-spectrum between processes $X$ and $Y$, and $f_{XX}(k)$ and $f_{YY}(k)$ are the spectra of processes $X$ and $Y$ respectively.

example_coh = magnitude_coherence(spec)
example_phase = phase(spec)

fig = Mke.Figure()
ax_coh = Mke.Axis(fig[1, 1], title = "Magnitude Coherence", aspect = 1)
Mke.heatmap!(ax_coh, example_coh[1, 2], colorrange = (0, 1))
ax_phase = Mke.Axis(fig[1, 2], title = "Phase", aspect = 1)
Mke.heatmap!(ax_phase, example_phase[1, 2], colorrange = (-π, π))
fig

References

This page was generated using Literate.jl.