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")

YAXArray Dataset
Shared Axes: 
  (↓ X Sampled{Float64} 562980.0:10.0:572970.0 ForwardOrdered Regular Points,
  → Y Sampled{Float64} 1.86002e6:-10.0:1.85003e6 ReverseOrdered Regular Points)

Variables: 
nir, red

Properties: Dict{String, Any}("projection" => "PROJCS[\"WGS 84 / UTM zone 36N\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",33],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32636\"]]")

Convert it into a DGGS dataset:

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)