YAXArrays.jl

This tutorial will illustrate how to use SpectralIndices.jl using YAXArrays.jl as input data.

First we need to download the data, like in the previous tutorial. Only this time the data is going to be higher dimensional and slightly more complex, hence the need for YAXArrays.jl. In order to do so we are going to use the load_dataset function:

using YAXArrays, DimensionalData
using SpectralIndices

yaxa = load_dataset("sentinel", YAXArray)

╭─────────────────────────────╮
│ 300×300×4 YAXArray{Int64,3} │
├─────────────────────────────┴───────────────────────────────────── dims ┐
  ↓ x     Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  → y     Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  ↗ bands Categorical{String} ["B02", "B03", "B04", "B08"] ForwardOrdered
├─────────────────────────────────────────────────────────────── metadata ┤
  Dict{String, Any}()
├────────────────────────────────────────────────────────────── file size ┤ 
  file size: 2.75 MB
└─────────────────────────────────────────────────────────────────────────┘

As it is possible to observe we have a YAXArray object with three dimensions: bands, x and y. Each band is one of the 10 m spectral bands of a Sentinel-2 image.

The data is stored as Int64, so let us convert it to Float and rescale it:

yaxa = yaxa./10000

╭───────────────────────────────╮
│ 300×300×4 YAXArray{Float64,3} │
├───────────────────────────────┴─────────────────────────────────── dims ┐
  ↓ x     Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  → y     Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  ↗ bands Categorical{String} ["B02", "B03", "B04", "B08"] ForwardOrdered
├─────────────────────────────────────────────────────────────── metadata ┤
  Dict{String, Any}()
├────────────────────────────────────────────────────────────── file size ┤ 
  file size: 2.75 MB
└─────────────────────────────────────────────────────────────────────────┘

Now let's compute the NDVI for this dataset!

First, let's define the bands to be used:

b8 = yaxa[bands = At("B08")]
b4 = yaxa[bands = At("B04")]

now, let's compute the index

ndvi_compute = compute_index("NDVI"; N=b8, R=b4)

╭─────────────────────────────╮
│ 300×300 YAXArray{Float64,2} │
├─────────────────────────────┴───────────────────── dims ┐
  ↓ x Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  → y Sampled{Int64} 1:300 ForwardOrdered Regular Points
├─────────────────────────────────────────────── metadata ┤
  Dict{String, Any}()
├────────────────────────────────────────────── file size ┤ 
  file size: 703.12 KB
└─────────────────────────────────────────────────────────┘

map

Due to the type design, in order to use map we will use the PartialFunctions package to specify the first initial type on each function, namely

using PartialFunctions
ndvi_p = NDVI.compute $Float64

NDVI_func(Float64, ...)

now, we can compute the index

ndvi_map = map(ndvi_p, b8, b4)

╭─────────────────────────────╮
│ 300×300 YAXArray{Float64,2} │
├─────────────────────────────┴───────────────────── dims ┐
  ↓ x Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  → y Sampled{Int64} 1:300 ForwardOrdered Regular Points
├─────────────────────────────────────────────── metadata ┤
  Dict{String, Any}()
├────────────────────────────────────────────── file size ┤ 
  file size: 703.12 KB
└─────────────────────────────────────────────────────────┘

Let's check that we have the same output:

ndvi_compute.data == ndvi_map.data

true

mapCube

For out of memory calculations then using mapCube is the way to go. This can be implemented as follows. First, we wrap our function of interest into a function compatible with mapCube, namely

function ndvi_out(xout, x1, x2)
    xout .= NDVI.(Float64, x1, x2) # note the .
end

ndvi_out (generic function with 1 method)

next, define the input and output dimensions of your YAXArray's.

in_dims = InDims("x") # the second one will be inferred
out_dims = OutDims("x") # dito

ndvi_cube = mapCube(ndvi_out, (b8, b4), indims=(in_dims, in_dims),
    outdims=OutDims("x", outtype=Float64))

╭─────────────────────────────────────────────╮
│ 300×300 YAXArray{Union{Missing, Float64},2} │
├─────────────────────────────────────────────┴───── dims ┐
  ↓ x Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  → y Sampled{Int64} 1:300 ForwardOrdered Regular Points
├─────────────────────────────────────────────── metadata ┤
  Dict{String, Any}()
├────────────────────────────────────────────── file size ┤ 
  file size: 703.12 KB
└─────────────────────────────────────────────────────────┘

and we check again the data output matches

ndvi_compute.data == ndvi_cube.data

true

Computing index by named dims

As usual we can also just feed a properly constructed YAXArray to the compute_index function. Let's built the array:

index_R = findfirst(yaxa.bands.val .== "B04")
index_N = findfirst(yaxa.bands.val .== "B08")
new_bands_dim = Dim{:Variables}(["R", "N"])

nr_data = cat(yaxa[:, :, index_R], yaxa[:, :, index_N], dims=3)

new_yaxa = YAXArray((yaxa.x, yaxa.y, new_bands_dim), nr_data)

╭───────────────────────────────╮
│ 300×300×2 YAXArray{Float64,3} │
├───────────────────────────────┴─────────────────────────── dims ┐
  ↓ x         Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  → y         Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  ↗ Variables Categorical{String} ["R", "N"] ReverseOrdered
├─────────────────────────────────────────────────────── metadata ┤
  Dict{String, Any}()
├────────────────────────────────────────────────────── file size ┤ 
  file size: 1.37 MB
└─────────────────────────────────────────────────────────────────┘

Warning

Please notice how the Dim is called Variables. This is needed for the internal computation to work properly. Also, note that this does not work for out of memory datasets.

Now that we have our YAXArray with the correctly names Dims we can use it direcly into compute_index:

ndvi = compute_index(
    "NDVI", new_yaxa
)

╭─────────────────────────────╮
│ 300×300 YAXArray{Float64,2} │
├─────────────────────────────┴───────────────────── dims ┐
  ↓ x Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  → y Sampled{Int64} 1:300 ForwardOrdered Regular Points
├─────────────────────────────────────────────── metadata ┤
  Dict{String, Any}()
├────────────────────────────────────────────── file size ┤ 
  file size: 703.12 KB
└─────────────────────────────────────────────────────────┘

Computing Kernels for kNDVI

We want to compute the kNDVI for a YAXArray.

kNDVI

kNDVI: Kernel Normalized Difference Vegetation Index
* Application Domain: kernel
* Bands/Parameters: Any["kNN", "kNR"]
* Formula: (kNN-kNR)/(kNN+kNR)
* Reference: https://doi.org/10.1126/sciadv.abc7447

As we see from bands we need the kNN and kNR. In order to compute these values SpectralIndices.jl provides a compute_kernel function. If you are curious about the kNDVI remember that you always have the source of the index in the reference field:

kNDVI.reference

"https://doi.org/10.1126/sciadv.abc7447"

Onto the calculations:

knn = YAXArray((yaxa.x, yaxa.y), fill(1.0, 300, 300));
knr = compute_kernel(
    RBF;
    a = Float64.(yaxa[bands = At("B08")]),
    b = Float64.(yaxa[bands = At("B04")]),
    sigma = yaxa[bands = At("B08")].+yaxa[bands = At("B04")] ./ 2
)

╭─────────────────────────────╮
│ 300×300 YAXArray{Float64,2} │
├─────────────────────────────┴───────────────────── dims ┐
  ↓ x Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  → y Sampled{Int64} 1:300 ForwardOrdered Regular Points
├─────────────────────────────────────────────── metadata ┤
  Dict{String, Any}()
├────────────────────────────────────────────── file size ┤ 
  file size: 703.12 KB
└─────────────────────────────────────────────────────────┘

As always, you can decide to build a YAXArray and feed that to the compute_kernel function if you prefer:

a = Float64.(yaxa[bands = At("B08")])
b = Float64.(yaxa[bands = At("B04")])
sigma = yaxa[bands = At("B08")].+yaxa[bands = At("B04")] ./ 2
kernel_dims = Dim{:Variables}(["a", "b", "sigma"])

params = concatenatecubes([a, b, sigma], kernel_dims)

knr = compute_kernel(RBF, params)

╭─────────────────────────────╮
│ 300×300 YAXArray{Float64,2} │
├─────────────────────────────┴───────────────────── dims ┐
  ↓ x Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  → y Sampled{Int64} 1:300 ForwardOrdered Regular Points
├─────────────────────────────────────────────── metadata ┤
  Dict{String, Any}()
├────────────────────────────────────────────── file size ┤ 
  file size: 703.12 KB
└─────────────────────────────────────────────────────────┘

We can finally compute the kNDVI:

kndvi = compute_index("kNDVI"; kNN = knn, kNR=knr)

╭─────────────────────────────╮
│ 300×300 YAXArray{Float64,2} │
├─────────────────────────────┴───────────────────── dims ┐
  ↓ x Sampled{Int64} 1:300 ForwardOrdered Regular Points,
  → y Sampled{Int64} 1:300 ForwardOrdered Regular Points
├─────────────────────────────────────────────── metadata ┤
  Dict{String, Any}()
├────────────────────────────────────────────── file size ┤ 
  file size: 703.12 KB
└─────────────────────────────────────────────────────────┘

Let's plot it!

using CairoMakie
fig, ax, plt = heatmap(kndvi; colormap=:haline,
    axis = (; aspect=DataAspect()),
    figure = (; size=(600, 400)))
Colorbar(fig[1,2], plt)
colsize!(fig.layout, 1, Aspect(1, 1.0))
fig