Mean Seasonal Cycle for a sigle 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()
lines(1:length(t), var; color = :purple, linewidth=1.25,
    axis=(; xlabel="Time", ylabel="Variable"),
    figure = (; resolution = (600,400))
    )
Example block output

Currently makie doesn't support time axis natively, but the following function can do the work for now.

function time_ticks(dates; frac=8)
    tempo = string.(dates)
    lentime = length(tempo)
    slice_dates = range(1, lentime, step=lentime ÷ frac)
    return slice_dates, tempo[slice_dates]
end
xpos, ticks = time_ticks(t; frac=8)

In order to apply the previous output, we split the plotting function into his 3 components, figure, axis and plotted object, namely

fig, ax, obj = lines(1:length(t), var; color = :purple, linewidth=1.25,
    axis=(; xlabel="Time", ylabel="Variable"),
    figure = (; resolution = (600,400))
    )
ax.xticks = (xpos, ticks)
ax.xticklabelrotation = π / 4
ax.xticklabelalign = (:right, :center)
fig
Example block output

Define the cube

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.026828628171312697
  0.04343008524971354
 -0.013963376578375446
  0.09632602646594195
  0.06785911352648893
  0.13640560842186403
 -0.019247143570796727
  0.15365962258524676
  0.2042693941687408
  0.16278960247700497
  ⋮
 -0.21315062702210044
 -0.083148807147485
 -0.07288337556769314
 -0.09166918551190778
 -0.06616128638997433
 -0.16221375342858915
 -0.020939015931497735
 -0.05280213695003191
 -0.06690457207246205
TODO: Apply the new groupby funtion from DD
Plot results: mean seasonal cycle
xpos, ticks = time_ticks(t[1:365]; frac=8)

fig, ax, obj = lines(1:365, var[1:365]; label="2021", color=:black,
    linewidth=2.0, linestyle=:dot,
    axis = (;  xlabel="Time", 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.xticks = (xpos, ticks)
ax.xticklabelrotation = π / 4
ax.xticklabelalign = (:right, :center)
fig
current_figure()
Example block output