SpatialMultitaper

SpatialMultitaper.jl provides tools for spectral analysis of spatial point processes and random fields using multitaper methods.

Installation

using Pkg
Pkg.add("SpatialMultitaper")

Quick Start

using SpatialMultitaper, GeoStatsProcesses

import GLMakie as Mke

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)

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}}

We can subset this object

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

An then plot transformations

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

We have other quantities that can be computed, either directly or from the spectra:

coh = coherence(spec)
Mke.heatmap(abs(coh[1,2]))
Example block output
rot_spec = rotational_estimate(spec)
Mke.lines(real(rot_spec[1,1]))
Example block output

Quick Start K function

using SpatialMultitaper, GeoStatsProcesses

import GLMakie as Mke

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)

tapers = sin_taper_family((4, 4), region)
nk = (100, 100)
kmax = (0.1, 0.1)
radii = range(0, 30, 100)
kfun = k_function(data; tapers = tapers, nk = nk, kmax = kmax, radii = radii)
K function of a 2 dimensional process
  between processes: 1:2 and 1:2
  evaluated at 0.0:0.30303030303030304:30.0
  with values of type Vector{StaticArraysCore.SMatrix{2, 2, Float64, 4}}

We can then plot the estimated cross K function as follows

Mke.lines(kfun[1,2])
Example block output

We can transform this to an L function

lfun = l_function(kfun)
L function of a 2 dimensional process
  between processes: 1:2 and 1:2
  evaluated at 0.0:0.30303030303030304:30.0
  with values of type Vector{StaticArraysCore.SMatrix{2, 2, Float64, 4}}

which can be easier to visualise

Mke.lines(lfun[1,2])
Example block output

Or we can compute this from the indexed version

lfun[1,2] == l_function(kfun[1,2])
true

Contributors