Calculate Sentinel-2 NDVI

Load and subset geographical data:

using DGGS
using YAXArrays
using ArchGDAL
using GLMakie

geo_ds = open_dataset("https://github.com/danlooo/DGGS.jl/raw/refs/heads/main/test/data/s2-ndvi-irrigation.tif")
plot(geo_ds.nir)

Convert it into a DGGS array:

dggs_ds = to_dggs_dataset(geo_ds, 19, geo_ds.nir.properties["projection"])

┌ 1175×937×1 DGGSDataset ┐
├────────────────────────┴─────────────────────────────────────────────── dims ┐
  ↓ dggs_i Sampled{Int64} 351619:352793 ForwardOrdered Regular Points,
  → dggs_j Sampled{Int64} 496491:497427 ForwardOrdered Regular Points,
  ↗ dggs_n Sampled{Int64} 2:2 ForwardOrdered Regular Points
├────────────────────────────────────────────────────────────────────── layers ┤
  :red eltype: Union{Missing, UInt16} dims: dggs_i, dggs_j, dggs_n size: 1175×937×1
  :nir eltype: Union{Missing, UInt16} dims: dggs_i, dggs_j, dggs_n size: 1175×937×1
├──────────────────────────────────────────────────────────────────────── DGGS ┤ 
  DGGSRS:     ISEA4D.Penta
  Resolution: 19 (up to 2.75e+12 cells)
  Geo BBox:   Extent(X = (33.590835539883145, 33.68487376375582), Y = (16.73225838289562, 16.822269591601994))
└──────────────────────────────────────────────────────────────────────────────┘

Calculate the NDVI:

uint16_max = 65535
nir = dggs_ds.nir ./ uint16_max
red = dggs_ds.red ./ uint16_max
ndvi = @. (nir - red) / (nir + red)

┌ 1175×937×1 DGGSArray{Union{Missing, Float64}, 3} ┐
├──────────────────────────────────────────────────┴───────────────────── dims ┐
  ↓ dggs_i Sampled{Int64} 351619:352793 ForwardOrdered Regular Points,
  → dggs_j Sampled{Int64} 496491:497427 ForwardOrdered Regular Points,
  ↗ dggs_n Sampled{Int64} 2:2 ForwardOrdered Regular Points
├──────────────────────────────────────────────────────────────────── metadata ┤
  Dict{String, Any} with 3 entries:
  "missing_value" => 65535.0
  "name"          => "nir"
  "projection"    => "PROJCS[\"WGS 84 / UTM zone 36N\",GEOGCS[\"WGS 84\",DATUM[…
├──────────────────────────────────────────────────────────────────────── DGGS ┤ 
  DGGSRS:     ISEA4D.Penta
  Resolution: 19 (up to 2.75e+12 cells)
  Geo BBox:   Extent(X = (33.590835539883145, 33.68487376375582), Y = (16.73225838289562, 16.822269591601994))
└──────────────────────────────────────────────────────────────────────────────┘

Hereby, we use the dot assignment macro @. to broadcast the function over all cells individually.

Plot the NDVI in DGGS:

plot(ndvi)