Mean Seasonal Cycle for a single pixel

using CairoMakie
CairoMakie.activate!()
using Dates
using Statistics

We define the data span. For simplicity, three non-leap years were selected.

t =  Date("2021-01-01"):Day(1):Date("2023-12-31")
NpY = 3

3

and create some seasonal dummy data

x = repeat(range(0, 2π, length=365), NpY)
var = @. sin(x) + 0.1 * randn()

fig, ax, obj = lines(t, var; color = :purple, linewidth=1.25,
    axis=(; xlabel="Time", ylabel="Variable"),
    figure = (; size = (600,400))
    )
ax.xticklabelrotation = π / 4
ax.xticklabelalign = (:right, :center)
fig

Define the cube

julia> using YAXArrays, DimensionalData

julia> using YAXArrays: YAXArrays as YAX

julia> axes = (YAX.Time(t),)

(↓ Time Date("2021-01-01"):Dates.Day(1):Date("2023-12-31"))

julia> c = YAXArray(axes, var)

┌ 1095-element YAXArray{Float64, 1} ┐
├───────────────────────────────────┴──────────────────────────────────── dims ┐
  ↓ Time Sampled{Date} Date("2021-01-01"):Dates.Day(1):Date("2023-12-31") ForwardOrdered Regular Points
├──────────────────────────────────────────────────────────────────── metadata ┤
  Dict{String, Any}()
├──────────────────────────────────────────────────────────── loaded in memory ┤
  data size: 8.55 KB
└──────────────────────────────────────────────────────────────────────────────┘

Let's calculate the mean seasonal cycle of our dummy variable 'var'

function mean_seasonal_cycle(c; ndays = 365)
    ## filterig by month-day
    monthday = map(x->Dates.format(x, "u-d"), collect(c.Time))
    datesid = unique(monthday)
    ## number of years
    NpY = Int(size(monthday,1)/ndays)
    idx = Int.(zeros(ndays, NpY))
    ## get the day-month indices for data subsetting
    for i in 1:ndays
        idx[i,:] = Int.(findall(x-> x == datesid[i], monthday))
    end
    ## compute the mean seasonal cycle
    mscarray = map(x->var[x], idx)
    msc = mapslices(mean, mscarray, dims=2)
    return msc
end

msc = mean_seasonal_cycle(c);

365×1 Matrix{Float64}:
  0.03528277758302477
  0.02345835017256051
  0.039451611552802975
  0.09689360224777122
  0.022312156353890087
  0.11391077060238619
  0.08058305464254185
  0.06159722707791853
  0.08158557886952912
  0.14715175267308206
  ⋮
 -0.10291543325743235
 -0.22990067443344916
 -0.12370988510072528
 -0.0750790675265741
 -0.05611581504766607
 -0.06417594925348342
 -0.062270000476910094
 -0.029247418895843032
 -0.05742042630154354

TODO: Apply the new groupby funtion from DD

Plot results: mean seasonal cycle

fig, ax, obj = lines(1:365, var[1:365]; label="2021", color=:black,
    linewidth=2.0, linestyle=:dot,
    axis = (;  xlabel="Day of Year", ylabel="Variable"),
    figure=(; size = (600,400))
    )
lines!(1:365, var[366:730], label="2022", color=:brown,
    linewidth=1.5, linestyle=:dash
    )
lines!(1:365, msc[:,1]; label="MSC", color=:dodgerblue, linewidth=2.5)
axislegend()
ax.xticklabelrotation = π / 4
ax.xticklabelalign = (:right, :center)
fig
current_figure()