SpectralIndices.jl

Appearance

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:

julia
using YAXArrays, DimensionalData
using SpectralIndices
julia
yaxa = load_dataset("sentinel", YAXArray)
┌ 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", …, "B08"] ForwardOrdered
├───────────────────────────────────────────── loaded in memory ┤
  data 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.

We now rescale the data:

julia
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", …, "B08"] ForwardOrdered
├──────────────────────────────────────────────── loaded lazily ┤
  data size: 2.75 MB
└───────────────────────────────────────────────────────────────┘

Now let's compute the NDVI for this dataset!

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

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

now, let's compute the index

julia
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
├─────────────────────────────────────────── loaded lazily ┤
  data 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

julia
using PartialFunctions
ndvi_p = NDVI.compute $Float64
NDVI_func(Float64, ...)

now, we can compute the index

julia
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
├─────────────────────────────────────────── loaded lazily ┤
  data size: 703.12 KB
└──────────────────────────────────────────────────────────┘

Let's check that we have the same output:

julia
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

julia
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.

julia
in_dims = InDims("x") # the second one will be inferred
out_dims = OutDims("x") # dito
julia
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
├──────────────────────────────────────── loaded in memory ┤
  data size: 703.12 KB
└──────────────────────────────────────────────────────────┘

and we check again the data output matches

julia
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:

julia
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)
julia
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
├────────────────────────────────────────── loaded in memory ┤
  data 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:

julia
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
├─────────────────────────────────────────── loaded lazily ┤
  data size: 703.12 KB
└──────────────────────────────────────────────────────────┘

Computing Kernels for kNDVI ​

We want to compute the kNDVI for a YAXArray.

julia
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:

julia
kNDVI.reference
"https://doi.org/10.1126/sciadv.abc7447"

Onto the calculations:

julia
knn = YAXArray((yaxa.x, yaxa.y), fill(1.0, 300, 300));
knr = compute_kernel(
    RBF;
    a = yaxa[bands = At("B08")],
    b = 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
├─────────────────────────────────────────── loaded lazily ┤
  data size: 703.12 KB
└──────────────────────────────────────────────────────────┘

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

julia
a = yaxa[bands = At("B08")]
b = 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)
julia
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
├─────────────────────────────────────────── loaded lazily ┤
  data size: 703.12 KB
└──────────────────────────────────────────────────────────┘

We can finally compute the kNDVI:

julia
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
├─────────────────────────────────────────── loaded lazily ┤
  data size: 703.12 KB
└──────────────────────────────────────────────────────────┘