Group YAXArrays and Datasets
The following examples will use the groupby
function to calculate temporal and spatial averages.
using YAXArrays, DimensionalData
using YAXArrays: YAXArrays as YAX
using NetCDF
using Downloads
using Dates
using Statistics
[ Info: new driver key :netcdf, updating backendlist.
Seasonal Averages from Time Series of Monthly Means
The following reproduces the example in xarray by Joe Hamman.
Where the goal is to calculate the seasonal average. And in order to do this properly, is necessary to calculate the weighted average considering that each month has a different number of days.
Download the data
url_path = "https://github.com/pydata/xarray-data/raw/master/rasm.nc"
filename = Downloads.download(url_path, "rasm.nc")
ds_o = Cube(filename)
┌ 275×205×36 YAXArray{Float64, 3} ┐
├─────────────────────────────────┴────────────────────────────────────── dims ┐
↓ x Sampled{Int64} 1:275 ForwardOrdered Regular Points,
→ y Sampled{Int64} 1:205 ForwardOrdered Regular Points,
↗ time Sampled{CFTime.DateTimeNoLeap} [CFTime.DateTimeNoLeap(1980-09-16T12:00:00), …, CFTime.DateTimeNoLeap(1983-08-17T00:00:00)] ForwardOrdered Irregular Points
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{String, Any} with 7 entries:
"units" => "C"
"coordinates" => "yc xc"
"name" => "Tair"
"long_name" => "Surface air temperature"
"type_preferred" => "double"
"_FillValue" => 9.96921e36
"time_rep" => "instantaneous"
├─────────────────────────────────────────────────────────────── loaded lazily ┤
data size: 15.48 MB
└──────────────────────────────────────────────────────────────────────────────┘
WARNING
The following rebuild should not be necessary in the future, plus is unpractical to use for large data sets. Out of memory groupby currently is work in progress. Related to https://github.com/rafaqz/DimensionalData.jl/issues/642
_FillValue = ds_o.properties["_FillValue"]
ds = replace(ds_o[:,:,:], _FillValue => NaN) # load into memory and replace _FillValue by NaN
┌ 275×205×36 YAXArray{Float64, 3} ┐
├─────────────────────────────────┴────────────────────────────────────── dims ┐
↓ x Sampled{Int64} 1:275 ForwardOrdered Regular Points,
→ y Sampled{Int64} 1:205 ForwardOrdered Regular Points,
↗ time Sampled{CFTime.DateTimeNoLeap} [CFTime.DateTimeNoLeap(1980-09-16T12:00:00), …, CFTime.DateTimeNoLeap(1983-08-17T00:00:00)] ForwardOrdered Irregular Points
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{String, Any} with 7 entries:
"units" => "C"
"coordinates" => "yc xc"
"name" => "Tair"
"long_name" => "Surface air temperature"
"type_preferred" => "double"
"_FillValue" => 9.96921e36
"time_rep" => "instantaneous"
├──────────────────────────────────────────────────────────── loaded in memory ┤
data size: 15.48 MB
└──────────────────────────────────────────────────────────────────────────────┘
GroupBy: seasons
function weighted_seasons(ds) ... end
function weighted_seasons(ds)
# calculate weights
tempo = dims(ds, :time)
month_length = YAXArray((tempo,), daysinmonth.(tempo))
g_tempo = groupby(month_length, YAX.time => seasons(; start=December))
sum_days = sum.(g_tempo, dims=:time)
weights = map(./, g_tempo, sum_days)
# unweighted seasons
g_ds = groupby(ds, YAX.time => seasons(; start=December))
mean_g = mean.(g_ds, dims=:time)
mean_g = dropdims.(mean_g, dims=:time)
# weighted seasons
g_dsW = broadcast_dims.(*, weights, g_ds)
weighted_g = sum.(g_dsW, dims = :time);
weighted_g = dropdims.(weighted_g, dims=:time)
# differences
diff_g = map(.-, weighted_g, mean_g)
seasons_g = lookup(mean_g, :time)
return mean_g, weighted_g, diff_g, seasons_g
end
INFO
In what follows, note how we are referencing the time dimension via YAX.time
. This approach is used to avoid name clashes with time
(Time
) from Base (Dates). For convenience, we have defined the Dimensions time
and Time
in YAXArrays.jl, which are only accessible when explicitly called.
Now, we continue with the groupby
operations as usual
julia> g_ds = groupby(ds, YAX.time => seasons(; start=December))
┌ 4-element DimGroupByArray{YAXArray{Float64,2},1} ┐
├──────────────────────────────────────────────────┴───────────────────── dims ┐
↓ time Categorical{Symbol} [:Dec_Jan_Feb, :Mar_Apr_May, :Jun_Jul_Aug, :Sep_Oct_Nov] Unordered
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{Symbol, Any} with 1 entry:
:groupby => :time=>CyclicBins(month; cycle=12, step=3, start=12)…
├────────────────────────────────────────────────────────────────── group dims ┤
↓ x, → y, ↗ time
└──────────────────────────────────────────────────────────────────────────────┘
:Dec_Jan_Feb 275×205×9 YAXArray
:Mar_Apr_May 275×205×9 YAXArray
:Jun_Jul_Aug 275×205×9 YAXArray
:Sep_Oct_Nov 275×205×9 YAXArray
And the mean per season is calculated as follows
julia> mean_g = mean.(g_ds, dims=:time)
┌ 4-element DimArray{YAXArray{Float64, 3, Array{Float64, 3}, Tuple{Dim{:x, DimensionalData.Dimensions.Lookups.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, Dim{:y, DimensionalData.Dimensions.Lookups.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, YAXArrays.time{DimensionalData.Dimensions.Lookups.Sampled{CFTime.DateTimeNoLeap, Vector{CFTime.DateTimeNoLeap}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Dict{String, Any}}, 1} ┐
├──────────────────────────────────────────────────────────────────────── dims ┤
↓ time Categorical{Symbol} [:Dec_Jan_Feb, :Mar_Apr_May, :Jun_Jul_Aug, :Sep_Oct_Nov] Unordered
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{Symbol, Any} with 1 entry:
:groupby => :time=>CyclicBins(month; cycle=12, step=3, start=12)…
└──────────────────────────────────────────────────────────────────────────────┘
:Dec_Jan_Feb … [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 11.1372 11.3835; NaN NaN … 11.3252 11.5843;;;]
:Mar_Apr_May [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 21.1363 21.018; NaN NaN … 21.4325 21.1762;;;]
:Jun_Jul_Aug [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 28.2818 27.9432; NaN NaN … 28.619 28.0537;;;]
:Sep_Oct_Nov [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 21.7119 21.7158; NaN NaN … 21.9682 21.9404;;;]
dropdims
Note that now the time dimension has length one, we can use dropdims
to remove it
julia> mean_g = dropdims.(mean_g, dims=:time)
┌ 4-element DimArray{YAXArray{Float64, 2, Matrix{Float64}, Tuple{Dim{:x, DimensionalData.Dimensions.Lookups.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, Dim{:y, DimensionalData.Dimensions.Lookups.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Dict{String, Any}}, 1} ┐
├──────────────────────────────────────────────────────────────────────── dims ┤
↓ time Categorical{Symbol} [:Dec_Jan_Feb, :Mar_Apr_May, :Jun_Jul_Aug, :Sep_Oct_Nov] Unordered
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{Symbol, Any} with 1 entry:
:groupby => :time=>CyclicBins(month; cycle=12, step=3, start=12)…
└──────────────────────────────────────────────────────────────────────────────┘
:Dec_Jan_Feb … [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 11.1372 11.3835; NaN NaN … 11.3252 11.5843]
:Mar_Apr_May [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 21.1363 21.018; NaN NaN … 21.4325 21.1762]
:Jun_Jul_Aug [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 28.2818 27.9432; NaN NaN … 28.619 28.0537]
:Sep_Oct_Nov [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 21.7119 21.7158; NaN NaN … 21.9682 21.9404]
seasons
Due to the groupby
function we will obtain new grouping names, in this case in the time dimension:
seasons_g = lookup(mean_g, :time)
Categorical{Symbol} Unordered
wrapping: 4-element Vector{Symbol}:
:Dec_Jan_Feb
:Mar_Apr_May
:Jun_Jul_Aug
:Sep_Oct_Nov
Next, we will weight this grouping by days/month in each group.
GroupBy: weight
Create a YAXArray
for the month length
tempo = dims(ds, :time)
month_length = YAXArray((tempo,), daysinmonth.(tempo))
┌ 36-element YAXArray{Int64, 1} ┐
├───────────────────────────────┴──────────────────────────────────────── dims ┐
↓ time Sampled{CFTime.DateTimeNoLeap} [CFTime.DateTimeNoLeap(1980-09-16T12:00:00), …, CFTime.DateTimeNoLeap(1983-08-17T00:00:00)] ForwardOrdered Irregular Points
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{String, Any}()
├──────────────────────────────────────────────────────────── loaded in memory ┤
data size: 288.0 bytes
└──────────────────────────────────────────────────────────────────────────────┘
Now group it by season
julia> g_tempo = groupby(month_length, YAX.time => seasons(; start=December))
┌ 4-element DimGroupByArray{YAXArray{Int64,0},1} ┐
├────────────────────────────────────────────────┴─────────────────────── dims ┐
↓ time Categorical{Symbol} [:Dec_Jan_Feb, :Mar_Apr_May, :Jun_Jul_Aug, :Sep_Oct_Nov] Unordered
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{Symbol, Any} with 1 entry:
:groupby => :time=>CyclicBins(month; cycle=12, step=3, start=12)…
├────────────────────────────────────────────────────────────────── group dims ┤
↓ time
└──────────────────────────────────────────────────────────────────────────────┘
:Dec_Jan_Feb 9-element YAXArray
:Mar_Apr_May 9-element YAXArray
:Jun_Jul_Aug 9-element YAXArray
:Sep_Oct_Nov 9-element YAXArray
Get the number of days per season
julia> sum_days = sum.(g_tempo, dims=:time)
┌ 4-element DimArray{YAXArray{Int64, 1, DimensionalData.DimVector{Int64, Tuple{YAXArrays.time{DimensionalData.Dimensions.Lookups.Sampled{CFTime.DateTimeNoLeap, Vector{CFTime.DateTimeNoLeap}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Tuple{}, Vector{Int64}, Symbol, DimensionalData.Dimensions.Lookups.NoMetadata}, Tuple{YAXArrays.time{DimensionalData.Dimensions.Lookups.Sampled{CFTime.DateTimeNoLeap, Vector{CFTime.DateTimeNoLeap}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Dict{String, Any}}, 1} ┐
├──────────────────────────────────────────────────────────────────────── dims ┤
↓ time Categorical{Symbol} [:Dec_Jan_Feb, :Mar_Apr_May, :Jun_Jul_Aug, :Sep_Oct_Nov] Unordered
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{Symbol, Any} with 1 entry:
:groupby => :time=>CyclicBins(month; cycle=12, step=3, start=12)…
└──────────────────────────────────────────────────────────────────────────────┘
:Dec_Jan_Feb [270]
:Mar_Apr_May [276]
:Jun_Jul_Aug [276]
:Sep_Oct_Nov [273]
weights
Weight the seasonal groups by sum_days
julia> weights = map(./, g_tempo, sum_days)
┌ 4-element DimArray{YAXArray{Float64, 1, Vector{Float64}, Tuple{YAXArrays.time{DimensionalData.Dimensions.Lookups.Sampled{CFTime.DateTimeNoLeap, Vector{CFTime.DateTimeNoLeap}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Dict{String, Any}}, 1} groupby ┐
├──────────────────────────────────────────────────────────────────────── dims ┤
↓ time Categorical{Symbol} [:Dec_Jan_Feb, :Mar_Apr_May, :Jun_Jul_Aug, :Sep_Oct_Nov] Unordered
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{Symbol, Any} with 1 entry:
:groupby => :time=>CyclicBins(month; cycle=12, step=3, start=12)…
└──────────────────────────────────────────────────────────────────────────────┘
:Dec_Jan_Feb … [0.114815, 0.114815, 0.103704, 0.114815, 0.114815, 0.103704, 0.114815, 0.114815, 0.103704]
:Mar_Apr_May [0.112319, 0.108696, 0.112319, 0.112319, 0.108696, 0.112319, 0.112319, 0.108696, 0.112319]
:Jun_Jul_Aug [0.108696, 0.112319, 0.112319, 0.108696, 0.112319, 0.112319, 0.108696, 0.112319, 0.112319]
:Sep_Oct_Nov [0.10989, 0.113553, 0.10989, 0.10989, 0.113553, 0.10989, 0.10989, 0.113553, 0.10989]
Verify that the sum per season is 1
julia> sum.(weights)
┌ 4-element DimArray{Float64, 1} ┐
├────────────────────────────────┴─────────────────────────────────────── dims ┐
↓ time Categorical{Symbol} [:Dec_Jan_Feb, :Mar_Apr_May, :Jun_Jul_Aug, :Sep_Oct_Nov] Unordered
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{Symbol, Any} with 1 entry:
:groupby => :time=>CyclicBins(month; cycle=12, step=3, start=12)…
└──────────────────────────────────────────────────────────────────────────────┘
:Dec_Jan_Feb 1.0
:Mar_Apr_May 1.0
:Jun_Jul_Aug 1.0
:Sep_Oct_Nov 1.0
weighted seasons
Now, let's weight the seasons
julia> g_dsW = broadcast_dims.(*, weights, g_ds)
┌ 4-element DimArray{YAXArray{Float64, 3, Array{Float64, 3}, Tuple{YAXArrays.time{DimensionalData.Dimensions.Lookups.Sampled{CFTime.DateTimeNoLeap, Vector{CFTime.DateTimeNoLeap}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, Dim{:x, DimensionalData.Dimensions.Lookups.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, Dim{:y, DimensionalData.Dimensions.Lookups.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, DimensionalData.Dimensions.Lookups.NoMetadata}, 1} ┐
├──────────────────────────────────────────────────────────────────────── dims ┤
↓ time Categorical{Symbol} [:Dec_Jan_Feb, :Mar_Apr_May, :Jun_Jul_Aug, :Sep_Oct_Nov] Unordered
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{Symbol, Any} with 1 entry:
:groupby => :time=>CyclicBins(month; cycle=12, step=3, start=12)…
└──────────────────────────────────────────────────────────────────────────────┘
:Dec_Jan_Feb … [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN;;; NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN;;; NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN;;; … ;;; NaN NaN … 1.32149 1.33565; NaN NaN … 1.29564 1.32555; … ; NaN NaN … 1.3188 1.3169; NaN NaN … 1.17863 1.18434;;; NaN NaN … 1.29816 1.34218; NaN NaN … 1.30113 1.35483; … ; NaN NaN … 1.30142 1.31753; NaN NaN … 1.16258 1.17647;;; NaN NaN … 1.34549 1.37878; NaN NaN … 1.36836 1.41634; … ; NaN NaN … 1.34832 1.38364; NaN NaN … 1.17852 1.16713]
:Mar_Apr_May [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN;;; NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN;;; NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN;;; … ;;; NaN NaN … 1.87705 1.90365; NaN NaN … 2.30018 2.35432; … ; NaN NaN … 2.41049 2.43254; NaN NaN … 2.65105 2.69085;;; NaN NaN … 1.86457 1.90712; NaN NaN … 2.2894 2.34818; … ; NaN NaN … 2.3866 2.41241; NaN NaN … 2.61197 2.64976;;; NaN NaN … 1.89237 1.8984; NaN NaN … 2.29473 2.312; … ; NaN NaN … 2.36142 2.36126; NaN NaN … 2.56632 2.59085]
:Jun_Jul_Aug [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN;;; NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN;;; NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN;;; … ;;; NaN NaN … 3.21209 3.25153; NaN NaN … 3.23 3.28008; … ; NaN NaN … 3.12575 3.15532; NaN NaN … 3.2434 3.26274;;; NaN NaN … 3.17434 3.21699; NaN NaN … 3.18892 3.24375; … ; NaN NaN … 3.06755 3.1083; NaN NaN … 3.19241 3.22211;;; NaN NaN … 3.1437 3.15644; NaN NaN … 3.16631 3.18583; … ; NaN NaN … 3.03361 3.05846; NaN NaN … 3.16581 3.16824]
:Sep_Oct_Nov [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN;;; NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN;;; NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … NaN NaN; NaN NaN … NaN NaN;;; … ;;; NaN NaN … 2.97047 3.00388; NaN NaN … 2.77587 2.80759; … ; NaN NaN … 2.60175 2.60918; NaN NaN … 1.4947 1.52419;;; NaN NaN … 2.94534 2.97649; NaN NaN … 2.75891 2.79502; … ; NaN NaN … 2.57695 2.59212; NaN NaN … 1.46506 1.49909;;; NaN NaN … 2.9192 2.93743; NaN NaN … 2.7593 2.77687; … ; NaN NaN … 2.57873 2.63006; NaN NaN … 1.48367 1.50089]
apply a sum
over the time dimension and drop it
julia> weighted_g = sum.(g_dsW, dims = :time);
julia> weighted_g = dropdims.(weighted_g, dims=:time)
┌ 4-element DimArray{YAXArray{Float64, 2, Matrix{Float64}, Tuple{Dim{:x, DimensionalData.Dimensions.Lookups.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, Dim{:y, DimensionalData.Dimensions.Lookups.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, DimensionalData.Dimensions.Lookups.NoMetadata}, 1} ┐
├──────────────────────────────────────────────────────────────────────── dims ┤
↓ time Categorical{Symbol} [:Dec_Jan_Feb, :Mar_Apr_May, :Jun_Jul_Aug, :Sep_Oct_Nov] Unordered
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{Symbol, Any} with 1 entry:
:groupby => :time=>CyclicBins(month; cycle=12, step=3, start=12)…
└──────────────────────────────────────────────────────────────────────────────┘
:Dec_Jan_Feb … [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 11.1181 11.372; NaN NaN … 11.3069 11.5743]
:Mar_Apr_May [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 21.1242 21.0057; NaN NaN … 21.4198 21.1644]
:Jun_Jul_Aug [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 28.2747 27.9362; NaN NaN … 28.6122 28.0465]
:Sep_Oct_Nov [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 21.73 21.7341; NaN NaN … 21.986 21.959]
Calculate the differences
julia> diff_g = map(.-, weighted_g, mean_g)
┌ 4-element DimArray{YAXArray{Float64, 2, Matrix{Float64}, Tuple{Dim{:x, DimensionalData.Dimensions.Lookups.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, Dim{:y, DimensionalData.Dimensions.Lookups.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, DimensionalData.Dimensions.Lookups.NoMetadata}, 1} ┐
├──────────────────────────────────────────────────────────────────────── dims ┤
↓ time Categorical{Symbol} [:Dec_Jan_Feb, :Mar_Apr_May, :Jun_Jul_Aug, :Sep_Oct_Nov] Unordered
├──────────────────────────────────────────────────────────────────── metadata ┤
Dict{Symbol, Any} with 1 entry:
:groupby => :time=>CyclicBins(month; cycle=12, step=3, start=12)…
└──────────────────────────────────────────────────────────────────────────────┘
:Dec_Jan_Feb … [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … -0.019016 -0.0115514; NaN NaN … -0.0183003 -0.00990356]
:Mar_Apr_May [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … -0.0121037 -0.0123091; NaN NaN … -0.0127077 -0.0117519]
:Jun_Jul_Aug [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … -0.00709111 -0.00693713; NaN NaN … -0.00684233 -0.00722034]
:Sep_Oct_Nov [NaN NaN … NaN NaN; NaN NaN … NaN NaN; … ; NaN NaN … 0.0180572 0.0182373; NaN NaN … 0.0178074 0.018571]
All the previous steps are equivalent to calling the function defined at the top:
mean_g, weighted_g, diff_g, seasons_g = weighted_seasons(ds)
Once all calculations are done we can plot the results with Makie.jl
as follows:
using CairoMakie
# define plot arguments/attributes
colorrange = (-30,30)
colormap = Reverse(:Spectral)
highclip = :red
lowclip = :grey15
cb_label = ds_o.properties["long_name"]
"Surface air temperature"
with_theme(theme_ggplot2()) do
hm_o, hm_d, hm_w = nothing, nothing, nothing
# the figure
fig = Figure(; size = (850,500))
axs = [Axis(fig[i,j], aspect=DataAspect()) for i in 1:3, j in 1:4]
for (j, s) in enumerate(seasons_g)
hm_o = heatmap!(axs[1,j], mean_g[time=At(s)]; colorrange, lowclip, highclip, colormap)
hm_w = heatmap!(axs[2,j], weighted_g[time=At(s)]; colorrange, lowclip, highclip, colormap)
hm_d = heatmap!(axs[3,j], diff_g[time=At(s)]; colorrange=(-0.1,0.1), lowclip, highclip,
colormap=:diverging_bwr_20_95_c54_n256)
end
Colorbar(fig[1:2,5], hm_o, label=cb_label)
Colorbar(fig[3,5], hm_d, label="Tair")
hidedecorations!.(axs, grid=false, ticks=false, label=false)
# some labels
[axs[1,j].title = string.(s) for (j,s) in enumerate(seasons_g)]
Label(fig[0,1:5], "Seasonal Surface Air Temperature", fontsize=18, font=:bold)
axs[1,1].ylabel = "Unweighted"
axs[2,1].ylabel = "Weighted"
axs[3,1].ylabel = "Difference"
colgap!(fig.layout, 5)
rowgap!(fig.layout, 5)
fig
end
which shows a good agreement with the results first published by Joe Hamman.