From eb405f8cc3fdbf42ded5995d3148471358e04359 Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Wed, 25 Sep 2024 16:24:47 +1000
Subject: [PATCH 01/15] Add metric_feature property to metrics

---
 src/metrics/metrics.jl      | 49 ++++++++++++++++++++++++++++---------
 src/metrics/reef_indices.jl | 15 +++++++-----
 src/metrics/scenario.jl     | 37 ++++++++++++++++++++++------
 src/metrics/utils.jl        |  7 +++---
 4 files changed, 80 insertions(+), 28 deletions(-)

diff --git a/src/metrics/metrics.jl b/src/metrics/metrics.jl
index 6f1b0885a..b326872bd 100644
--- a/src/metrics/metrics.jl
+++ b/src/metrics/metrics.jl
@@ -29,10 +29,12 @@ const IS_NOT_RELATIVE = false
 struct Metric{F<:Function,T<:Tuple,S<:String,B<:Bool} <: Outcome
     func::F
     dims::T     # output dimension axes ?
+    feature::S
     is_relative::B
     unit::S
 end
-Metric(f, d, r) = Metric(f, d, r, "")
+
+Metric(func, dims, feature, is_relative) = Metric(func, dims, feature, is_relative, "")
 
 """
     (f::Metric)(raw, args...; kwargs...)
@@ -70,7 +72,9 @@ end
 function _relative_cover(rs::ResultSet)::YAXArray{<:Real}
     return rs.outcomes[:relative_cover]
 end
-relative_cover = Metric(_relative_cover, (:timesteps, :locations, :scenarios), IS_RELATIVE)
+relative_cover = Metric(
+    _relative_cover, (:timesteps, :locations, :scenarios), "Relative Cover", IS_RELATIVE
+)
 
 """
     total_absolute_cover(X::AbstractArray{<:Real}, k_area::Vector{<:Real})::AbstractArray{<:Real}
@@ -96,7 +100,11 @@ function _total_absolute_cover(rs::ResultSet)::AbstractArray{<:Real}
     return _total_absolute_cover(rs.outcomes[:relative_cover], site_k_area(rs))
 end
 total_absolute_cover = Metric(
-    _total_absolute_cover, (:timesteps, :locations, :scenarios), IS_NOT_RELATIVE, UNIT_AREA
+    _total_absolute_cover,
+    (:timesteps, :locations, :scenarios),
+    "Cover",
+    IS_NOT_RELATIVE,
+    UNIT_AREA
 )
 
 """
@@ -144,7 +152,7 @@ function _relative_taxa_cover(rs::ResultSet)::AbstractArray{<:Real,3}
     return rs.outcomes[:relative_taxa_cover]
 end
 relative_taxa_cover = Metric(
-    _relative_taxa_cover, (:timesteps, :species, :scenarios), IS_RELATIVE
+    _relative_taxa_cover, (:timesteps, :species, :scenarios), "Cover", IS_RELATIVE
 )
 
 """
@@ -183,7 +191,10 @@ function _relative_loc_taxa_cover(
     return replace!(taxa_cover, NaN => 0.0)
 end
 relative_loc_taxa_cover = Metric(
-    _relative_loc_taxa_cover, (:timesteps, :species, :locations, :scenarios), IS_RELATIVE
+    _relative_loc_taxa_cover,
+    (:timesteps, :species, :locations, :scenarios),
+    "Relative Cover",
+    IS_RELATIVE
 )
 
 """
@@ -210,7 +221,7 @@ function _relative_juveniles(rs::ResultSet)::AbstractArray{<:Real,3}
     return rs.outcomes[:relative_juveniles]
 end
 relative_juveniles = Metric(
-    _relative_juveniles, (:timesteps, :locations, :scenarios), IS_RELATIVE
+    _relative_juveniles, (:timesteps, :locations, :scenarios), "Relative Cover", IS_RELATIVE
 )
 
 """
@@ -234,7 +245,11 @@ function _absolute_juveniles(rs::ResultSet)::AbstractArray{<:Real,3}
     return rs.outcomes[:relative_juveniles] .* site_k_area(rs)'
 end
 absolute_juveniles = Metric(
-    _absolute_juveniles, (:timesteps, :locations, :scenarios), IS_NOT_RELATIVE, UNIT_AREA
+    _absolute_juveniles,
+    (:timesteps, :locations, :scenarios),
+    "Cover",
+    IS_NOT_RELATIVE,
+    UNIT_AREA
 )
 
 """
@@ -282,7 +297,10 @@ function _juvenile_indicator(rs::ResultSet)::AbstractArray{<:Real,3}
     return rs.outcomes[:juvenile_indicator]
 end
 juvenile_indicator = Metric(
-    _juvenile_indicator, (:timesteps, :locations, :scenarios), IS_NOT_RELATIVE,
+    _juvenile_indicator,
+    (:timesteps, :locations, :scenarios),
+    "Density Indicator",
+    IS_NOT_RELATIVE,
     UNIT_AREA_INVERSE)
 
 """
@@ -323,7 +341,10 @@ function _coral_evenness(rs::ResultSet)::AbstractArray{<:Real,3}
     return rs.outcomes[:coral_evenness]
 end
 coral_evenness = Metric(
-    _coral_evenness, (:timesteps, :locations, :scenarios), IS_NOT_RELATIVE
+    _coral_evenness,
+    (:timesteps, :locations, :scenarios),
+    "Evenness Indicator",
+    IS_NOT_RELATIVE
 )
 
 """
@@ -567,7 +588,10 @@ function _absolute_shelter_volume(rs::ResultSet)::AbstractArray
     return rs.outcomes[:absolute_shelter_volume]
 end
 absolute_shelter_volume = Metric(
-    _absolute_shelter_volume, (:timesteps, :locations, :scenarios), IS_NOT_RELATIVE,
+    _absolute_shelter_volume,
+    (:timesteps, :locations, :scenarios),
+    "Volume",
+    IS_NOT_RELATIVE,
     UNIT_VOLUME
 )
 
@@ -702,7 +726,10 @@ function _relative_shelter_volume(rs::ResultSet)::YAXArray
     return rs.outcomes[:relative_shelter_volume]
 end
 relative_shelter_volume = Metric(
-    _relative_shelter_volume, (:timesteps, :locations, :scenarios), IS_RELATIVE
+    _relative_shelter_volume,
+    (:timesteps, :locations, :scenarios),
+    "Relative Volume",
+    IS_RELATIVE
 )
 
 include("pareto.jl")
diff --git a/src/metrics/reef_indices.jl b/src/metrics/reef_indices.jl
index 5e9657b32..99b508b16 100644
--- a/src/metrics/reef_indices.jl
+++ b/src/metrics/reef_indices.jl
@@ -104,7 +104,10 @@ function _reef_condition_index(rs::ResultSet)::AbstractArray{<:Real}
     return _reef_condition_index(rc, evenness, sv, juves)
 end
 reef_condition_index = Metric(
-    _reef_condition_index, (:timesteps, :locations, :scenarios), IS_NOT_RELATIVE
+    _reef_condition_index,
+    (:timesteps, :locations, :scenarios),
+    "RCI",
+    IS_NOT_RELATIVE
 )
 
 """
@@ -127,7 +130,7 @@ function _scenario_rci(rs::ResultSet; kwargs...)
 
     return _scenario_rci(rci, tac; kwargs...)
 end
-scenario_rci = Metric(_scenario_rci, (:timesteps, :scenarios), IS_NOT_RELATIVE)
+scenario_rci = Metric(_scenario_rci, (:timesteps, :scenarios), "RCI", IS_NOT_RELATIVE)
 
 """
     reef_tourism_index(rc::AbstractArray, evenness::AbstractArray, sv::AbstractArray, juves::AbstractArray, intcp_u::Vector)::AbstractArray
@@ -192,7 +195,7 @@ function _reef_tourism_index(rs::ResultSet; intcp_u::Bool=false)::AbstractArray
     return _reef_tourism_index(rc, evenness, sv, juves, intcp)
 end
 reef_tourism_index = Metric(
-    _reef_tourism_index, (:timesteps, :locations, :scenarios), IS_NOT_RELATIVE
+    _reef_tourism_index, (:timesteps, :locations, :scenarios), "RTI", IS_NOT_RELATIVE
 )
 
 """
@@ -208,7 +211,7 @@ end
 function _scenario_rti(rs::ResultSet; kwargs...)
     return _scenario_rti(reef_tourism_index(rs); kwargs...)
 end
-scenario_rti = Metric(_scenario_rti, (:timesteps, :scenarios), IS_NOT_RELATIVE)
+scenario_rti = Metric(_scenario_rti, (:timesteps, :scenarios), "RTI", IS_NOT_RELATIVE)
 
 """
     reef_fish_index(rc::AbstractArray)
@@ -276,7 +279,7 @@ function _reef_fish_index(rs::ResultSet; intcp_u1::Bool=false, intcp_u2::Bool=fa
     return _reef_fish_index(rc, icp1, icp2)
 end
 reef_fish_index = Metric(
-    _reef_fish_index, (:timesteps, :locations, :scenarios), IS_NOT_RELATIVE
+    _reef_fish_index, (:timesteps, :locations, :scenarios), "RFI", IS_NOT_RELATIVE
 )
 
 """
@@ -292,4 +295,4 @@ end
 function _scenario_rfi(rs::ResultSet; kwargs...)
     return _scenario_rfi(reef_fish_index(rs); kwargs...)
 end
-scenario_rfi = Metric(_scenario_rfi, (:timesteps, :scenarios), IS_NOT_RELATIVE)
+scenario_rfi = Metric(_scenario_rfi, (:timesteps, :scenarios), "RFI", IS_NOT_RELATIVE)
diff --git a/src/metrics/scenario.jl b/src/metrics/scenario.jl
index 70e311e05..868fd4906 100644
--- a/src/metrics/scenario.jl
+++ b/src/metrics/scenario.jl
@@ -49,7 +49,7 @@ function _scenario_total_cover(rs::ResultSet; kwargs...)::AbstractArray{<:Real}
     return _scenario_total_cover(tac::AbstractArray; kwargs...)
 end
 scenario_total_cover = Metric(
-    _scenario_total_cover, (:timesteps, :scenarios), IS_NOT_RELATIVE, UNIT_AREA
+    _scenario_total_cover, (:timesteps, :scenarios), "Cover", IS_NOT_RELATIVE, UNIT_AREA
 )
 
 """
@@ -64,7 +64,7 @@ function _scenario_relative_cover(rs::ResultSet; kwargs...)::AbstractArray{<:Rea
     return _scenario_total_cover(rs; kwargs...) ./ target_area
 end
 scenario_relative_cover = Metric(
-    _scenario_relative_cover, (:timesteps, :scenarios), IS_RELATIVE
+    _scenario_relative_cover, (:timesteps, :scenarios), "Relative Cover", IS_RELATIVE
 )
 
 """
@@ -107,7 +107,10 @@ function _scenario_relative_juveniles(rs::ResultSet; kwargs...)::YAXArray
     return dropdims(sum(aj; dims=:locations); dims=:locations) ./ sum(site_k_area(rs))
 end
 scenario_relative_juveniles = Metric(
-    _scenario_relative_juveniles, (:timesteps, :scenarios), IS_RELATIVE
+    _scenario_relative_juveniles,
+    (:timesteps, :scenarios),
+    "Relative Juveniles",
+    IS_RELATIVE
 )
 
 """
@@ -136,7 +139,11 @@ function _scenario_absolute_juveniles(rs::ResultSet; kwargs...)::AbstractArray{<
     return dropdims(sum(absolute_juveniles(rs); dims=:locations); dims=:locations)
 end
 scenario_absolute_juveniles = Metric(
-    _scenario_absolute_juveniles, (:timesteps, :scenarios), IS_NOT_RELATIVE, UNIT_AREA
+    _scenario_absolute_juveniles,
+    (:timesteps, :scenarios),
+    "Number of Juveniles",
+    IS_NOT_RELATIVE,
+    UNIT_AREA
 )
 
 """
@@ -164,7 +171,10 @@ function _scenario_juvenile_indicator(rs::ResultSet; kwargs...)::AbstractArray{<
     return dropdims(mean(juvenile_indicator(rs); dims=:locations); dims=:locations)
 end
 scenario_juvenile_indicator = Metric(
-    _scenario_juvenile_indicator, (:timesteps, :scenarios), IS_RELATIVE
+    _scenario_juvenile_indicator,
+    (:timesteps, :scenarios),
+    "Juvenile Indicator",
+    IS_RELATIVE
 )
 
 """
@@ -185,7 +195,11 @@ function _scenario_asv(rs::ResultSet; kwargs...)::AbstractArray{<:Real}
     return _scenario_asv(rs.outcomes[:absolute_shelter_volume]; kwargs...)
 end
 scenario_asv = Metric(
-    _scenario_asv, (:timesteps, :scenarios), IS_NOT_RELATIVE, "$UNIT_VOLUME/$UNIT_AREA"
+    _scenario_asv,
+    (:timesteps, :scenarios),
+    "Volume",
+    IS_NOT_RELATIVE,
+    "$UNIT_VOLUME/$UNIT_AREA"
 )
 
 """
@@ -201,7 +215,9 @@ end
 function _scenario_rsv(rs::ResultSet; kwargs...)::AbstractArray{<:Real}
     return _scenario_rsv(rs.outcomes[:relative_shelter_volume]; kwargs...)
 end
-scenario_rsv = Metric(_scenario_rsv, (:timesteps, :scenarios), IS_RELATIVE)
+scenario_rsv = Metric(
+    _scenario_rsv, (:timesteps, :scenarios), "Relative Volume", IS_RELATIVE
+)
 
 """
     scenario_evenness(ev::YAXArray; kwargs...)::AbstractArray{<:Real}
@@ -217,7 +233,12 @@ end
 function _scenario_evenness(rs::ResultSet; kwargs...)::AbstractArray{<:Real}
     return _scenario_evenness(rs.outcomes[:coral_evenness]; kwargs...)
 end
-scenario_evenness = Metric(_scenario_evenness, (:timesteps, :scenarios), IS_NOT_RELATIVE)
+scenario_evenness = Metric(
+    _scenario_evenness,
+    (:timesteps, :scenarios),
+    "Evenness Indicator",
+    IS_NOT_RELATIVE
+)
 
 """
     scenario_outcomes(rs::ResultSet, metrics::Vector{Metric})::YAXArray
diff --git a/src/metrics/utils.jl b/src/metrics/utils.jl
index 4ab23b03d..688fb84c2 100644
--- a/src/metrics/utils.jl
+++ b/src/metrics/utils.jl
@@ -50,7 +50,7 @@ end
 Units for each metric axis.
 """
 function axes_units(axes_names::Union{Vector{Symbol},Tuple})::Tuple
-    return values((timesteps="year", species="", locations="", scenarios="")[axes_names])
+    return values((timesteps="years", species="", locations="", scenarios="")[axes_names])
 end
 
 """
@@ -101,9 +101,10 @@ Fill `:axes_names` and `:axes_units` properties of the datacube.
 - `datacube` : YAXArray datacube
 """
 function fill_axes_properties(metric::Metric, metric_result::YAXArray)::YAXArray
-    metric_result.properties[:metric_name] = metric_label(metric)
-    metric_result.properties[:metric_unit] = metric.unit
+    metric_result.properties[:metric_name] = to_string(metric; is_titlecase=true)
+    metric_result.properties[:metric_feature] = metric.feature
     metric_result.properties[:is_relative] = metric.is_relative
+    metric_result.properties[:metric_unit] = metric.unit
 
     _axes_names::Tuple = axes_names(metric_result)
     metric_result.properties[:axes_names] = collect(

From 7776ff59c0cedaf09d25965c953f5dad1aacc166 Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Wed, 25 Sep 2024 16:25:35 +1000
Subject: [PATCH 02/15] Use metrics metadata to fill scenarios and clustering
 viz functions

---
 ext/AvizExt/outcome_metadata.jl     | 14 ++++++++++++++
 ext/AvizExt/viz/clustering.jl       |  4 ++++
 ext/AvizExt/viz/scenarios.jl        | 12 ++++++++----
 ext/AvizExt/viz/viz.jl              |  1 +
 test/metrics/test_metrics_helper.jl | 13 ++++++-------
 5 files changed, 33 insertions(+), 11 deletions(-)
 create mode 100644 ext/AvizExt/outcome_metadata.jl

diff --git a/ext/AvizExt/outcome_metadata.jl b/ext/AvizExt/outcome_metadata.jl
new file mode 100644
index 000000000..21be00383
--- /dev/null
+++ b/ext/AvizExt/outcome_metadata.jl
@@ -0,0 +1,14 @@
+function outcome_title(outcomes::YAXArray)::String
+    return get(outcomes.properties, :metric_name, "")
+end
+
+function outcome_label(outcomes::YAXArray)::String
+    outcome_metadata = outcomes.properties
+
+    return if haskey(outcome_metadata, :metric_feature) &&
+        haskey(outcome_metadata, :metric_unit)
+        "$(outcome_metadata[:metric_feature]) [$(outcome_metadata[:metric_unit])]"
+    else
+        ""
+    end
+end
diff --git a/ext/AvizExt/viz/clustering.jl b/ext/AvizExt/viz/clustering.jl
index 558a97434..9ae471da8 100644
--- a/ext/AvizExt/viz/clustering.jl
+++ b/ext/AvizExt/viz/clustering.jl
@@ -79,6 +79,10 @@ function ADRIA.viz.clustered_scenarios(
     fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
     axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
 )::Figure
+    if !haskey(axis_opts, :title)
+        axis_opts[:title] = "$(outcome_title(outcomes)) Clusters"
+    end
+
     return ADRIA.viz.scenarios(
         outcomes, clusters; opts=opts, fig_opts=fig_opts, axis_opts=axis_opts
     )
diff --git a/ext/AvizExt/viz/scenarios.jl b/ext/AvizExt/viz/scenarios.jl
index 9ae3024bf..c35fc3372 100644
--- a/ext/AvizExt/viz/scenarios.jl
+++ b/ext/AvizExt/viz/scenarios.jl
@@ -129,6 +129,11 @@ function ADRIA.viz.scenarios!(
     # Ensure last year is always shown in x-axis
     xtick_vals = get(axis_opts, :xticks, _time_labels(timesteps(outcomes)))
     xtick_rot = get(axis_opts, :xticklabelrotation, 2 / π)
+
+    if !haskey(axis_opts, :title)
+        axis_opts[:title] = outcome_title(outcomes)
+    end
+
     ax = Axis(g[1, 1]; xticks=xtick_vals, xticklabelrotation=xtick_rot, axis_opts...)
 
     _scenarios = copy(scenarios[1:end .∈ [outcomes.scenarios], :])
@@ -179,9 +184,8 @@ function ADRIA.viz.scenarios!(
         _render_legend(g, scen_groups, legend_position, legend_labels)
     end
 
-    ax.xlabel = "Year"
-    # ax.ylabel = metric_label(metric)
-
+    ax.xlabel = "Time [years]"
+    ax.ylabel = outcome_label(outcomes)
     return g
 end
 
@@ -198,7 +202,7 @@ function _confints(
     agg_dim = symdiff(axes_names(outcomes), [:timesteps])[1]
     for (idx, group) in enumerate(group_names)
         confints[:, idx, :] = series_confint(
-            outcomes[:, scen_groups[group]]; agg_dim=agg_dim
+            outcomes.data[:, scen_groups[group]]; agg_dim=agg_dim
         )
     end
 
diff --git a/ext/AvizExt/viz/viz.jl b/ext/AvizExt/viz/viz.jl
index 140539677..ea09b61bb 100644
--- a/ext/AvizExt/viz/viz.jl
+++ b/ext/AvizExt/viz/viz.jl
@@ -85,6 +85,7 @@ function _calc_gridsize(n_factors::Int64; max_cols::Int64=4)::Tuple{Int64,Int64}
     return n_rows, n_cols
 end
 
+include("../outcome_metadata.jl")
 include("scenarios.jl")
 include("sensitivity.jl")
 include("clustering.jl")
diff --git a/test/metrics/test_metrics_helper.jl b/test/metrics/test_metrics_helper.jl
index 63c2cc036..3abc14d7e 100644
--- a/test/metrics/test_metrics_helper.jl
+++ b/test/metrics/test_metrics_helper.jl
@@ -42,16 +42,15 @@ function test_metric(metric::metrics.Metric, params::Tuple)::Nothing
     return nothing
 end
 
-function _test_properties(metric::metrics.Metric, metric_result::YAXArray)
+function _test_properties(metric::metrics.Metric, outcomes::YAXArray)
     prop_labels::NTuple{4,Symbol} = (:metric_name, :metric_unit, :axes_names, :axes_units)
-    all_properties_present = all(haskey.([metric_result.properties], prop_labels))
+    outcome_metadata = outcomes.properties
+    all_properties_present = all(haskey.([outcome_metadata], prop_labels))
     @test all_properties_present
     if all_properties_present
-        result_prop_axes_names::Vector{String} = metric_result.properties[:axes_names]
-        result_prop_axes_units::Vector{String} = metric_result.properties[:axes_units]
-        result_axes_names = string.(axes_names(metric_result))
-
-        @test metric_result.properties[:metric_name] == metrics.metric_label(metric)
+        result_prop_axes_names::Vector{String} = outcome_metadata[:axes_names]
+        result_prop_axes_units::Vector{String} = outcome_metadata[:axes_units]
+        result_axes_names = string.(axes_names(outcomes))
 
         # All keys in result properties match YAXArray
         @test all(result_prop_axes_names .∈ [result_axes_names]) ||

From e9909241f6d7306b443ed5694879f7ad35307390 Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Fri, 27 Sep 2024 18:11:09 +1000
Subject: [PATCH 03/15] Extract metadata related code to metrics/metadata.jl
 and refactor

---
 src/metrics/metadata.jl                   | 72 +++++++++++++++++++++++
 src/metrics/metrics.jl                    |  8 +--
 src/metrics/{site_level.jl => spatial.jl} |  6 +-
 src/metrics/utils.jl                      | 33 -----------
 4 files changed, 80 insertions(+), 39 deletions(-)
 create mode 100644 src/metrics/metadata.jl
 rename src/metrics/{site_level.jl => spatial.jl} (95%)

diff --git a/src/metrics/metadata.jl b/src/metrics/metadata.jl
new file mode 100644
index 000000000..85e24ad16
--- /dev/null
+++ b/src/metrics/metadata.jl
@@ -0,0 +1,72 @@
+const METRIC_METADATA = (:metric_name, :metric_feature, :is_relative, :metric_unit)
+const AXIS_METADATA = (:axes_names, :axes_units)
+
+function metadata(outcomes::YAXArray)::Dict{Symbol,Any}
+    return outcomes.properties
+end
+
+"""
+    fill_metadata!(outcomes::YAXArray, metric::Metric)::Nothing
+    fill_metadata!(outcomes::YAXArray; kwargs...)::Nothing
+
+Fill `:axes_names` and `:axes_units` properties of the datacube.
+
+# Arguments
+- `outcomes` : YAXArray datacube of metric outcomes
+- `metric` : ADRIA.metrics.Metric object
+"""
+function fill_metadata!(
+    outcomes::YAXArray{T,N,A}, metric::Metric
+)::YAXArray{T,N,A} where {T,N,A}
+    fill_axes_metadata!(@views(outcomes))
+
+    outcomes.properties[:metric_name] = to_string(metric; is_titlecase=true)
+    outcomes.properties[:metric_feature] = metric.feature
+    outcomes.properties[:is_relative] = metric.is_relative
+    outcomes.properties[:metric_unit] = metric.unit
+
+    return outcomes
+end
+function fill_metadata!(
+    outcomes::YAXArray{T,N,A}, metadata::Dict{Symbol,Any}
+)::YAXArray{T,N,A} where {T,N,A}
+    # If `outcomes.properties`` is not a Dict{Symbol,Any} we need to build a new cube
+    # because `properties`` type can't be updated
+    if !(typeof(outcomes.properties) == Dict{Symbol,Any})
+        _axis_names = axes_names(outcomes)
+        _axis_labels = axis_labels.([outcomes], _axis_names)
+        outcomes = DataCube(outcomes.data; NamedTuple{_axis_names}(_axis_labels)...)
+    end
+
+    fill_axes_metadata!(@views(outcomes))
+
+    for metric_metadata in METRIC_METADATA
+        if haskey(metadata, metric_metadata)
+            outcomes.properties[metric_metadata] = metadata[metric_metadata]
+        else
+            @warn "Metric metadata \"$metric_metadata\" not found and won't be filled."
+        end
+    end
+
+    return outcomes
+end
+
+function fill_axes_metadata!(outcomes::YAXArray)::Nothing
+    _axes_names::Tuple = axes_names(outcomes)
+    outcomes.properties[:axes_names] = collect(
+        parentmodule(metrics).human_readable_name.(
+            _axes_names
+        )
+    )
+    outcomes.properties[:axes_units] = collect(axes_units(_axes_names))
+    return nothing
+end
+
+"""
+    axes_units(axes_names::Tuple)::Tuple
+
+Units for each metric axis.
+"""
+function axes_units(axes_names::Union{Vector{Symbol},Tuple})::Tuple
+    return values((timesteps="years", species="", locations="", scenarios="")[axes_names])
+end
diff --git a/src/metrics/metrics.jl b/src/metrics/metrics.jl
index b326872bd..ee4375187 100644
--- a/src/metrics/metrics.jl
+++ b/src/metrics/metrics.jl
@@ -44,11 +44,10 @@ Makes Metric types callable with arbitary arguments that are passed to associate
 """
 function (f::Metric)(raw, args...; kwargs...)::YAXArray
     axes::Tuple = (:timesteps, :species, :locations, :scenarios)[1:ndims(raw)]
-
-    return fill_axes_properties(f, f.func(DataCube(raw, axes), args...; kwargs...))
+    return fill_metadata!(@views(f.func(DataCube(raw, axes), args...; kwargs...)), f)
 end
 function (f::Metric)(rs::ResultSet, args...; kwargs...)::YAXArray
-    return fill_axes_properties(f, (f.func(rs, args...; kwargs...)))
+    return fill_metadata!(@views(f.func(rs, args...; kwargs...)), f)
 end
 
 """
@@ -732,11 +731,12 @@ relative_shelter_volume = Metric(
     IS_RELATIVE
 )
 
+include("metadata.jl")
 include("pareto.jl")
 include("ranks.jl")
 include("reef_indices.jl")
 include("scenario.jl")
-include("site_level.jl")
+include("spatial.jl")
 include("temporal.jl")
 include("utils.jl")
 
diff --git a/src/metrics/site_level.jl b/src/metrics/spatial.jl
similarity index 95%
rename from src/metrics/site_level.jl
rename to src/metrics/spatial.jl
index d2edaf1a6..47710070f 100644
--- a/src/metrics/site_level.jl
+++ b/src/metrics/spatial.jl
@@ -96,7 +96,7 @@ function summarize(
     # Only use this approach when data occupies more than 70% of available RAM
     if data_size > available_ram
         # `D.` is ensuring the returned YAXArray has the same type as the input `data`
-        return D.(mapslices(metric, data; dims=alongs_axis))
+        return fill_metadata!(D.(mapslices(metric, data; dims=alongs_axis)), metadata(data))
     end
 
     alongs = sort([axis_index(data, axis) for axis in alongs_axis])
@@ -110,7 +110,9 @@ function summarize(
     new_dims = setdiff(axes_names(data), alongs_axis)
     new_axis = [axis_labels(data, ax) for ax in new_dims]
 
-    return DataCube(summarized_data; NamedTuple{Tuple(new_dims)}(new_axis)...)
+    return fill_metadata!(
+        DataCube(summarized_data; NamedTuple{Tuple(new_dims)}(new_axis)...), metadata(data)
+    )
 end
 function summarize(
     data::YAXArray{D,T,N,A},
diff --git a/src/metrics/utils.jl b/src/metrics/utils.jl
index 688fb84c2..981ebc0fa 100644
--- a/src/metrics/utils.jl
+++ b/src/metrics/utils.jl
@@ -44,15 +44,6 @@ function metric_label(f::Function, unit::String)::String
     return n
 end
 
-"""
-    axes_units(axes_names::Tuple)::Tuple
-
-Units for each metric axis.
-"""
-function axes_units(axes_names::Union{Vector{Symbol},Tuple})::Tuple
-    return values((timesteps="years", species="", locations="", scenarios="")[axes_names])
-end
-
 """
     dims(m::Metric)::Tuple
 
@@ -92,30 +83,6 @@ function call_metric(metric::Union{Function,Metric}, data::YAXArray, args...; kw
     end
 end
 
-"""
-    fill_axes_properties(metric::Metric, metric_result::YAXArray)::YAXArray
-
-Fill `:axes_names` and `:axes_units` properties of the datacube.
-
-# Arguments
-- `datacube` : YAXArray datacube
-"""
-function fill_axes_properties(metric::Metric, metric_result::YAXArray)::YAXArray
-    metric_result.properties[:metric_name] = to_string(metric; is_titlecase=true)
-    metric_result.properties[:metric_feature] = metric.feature
-    metric_result.properties[:is_relative] = metric.is_relative
-    metric_result.properties[:metric_unit] = metric.unit
-
-    _axes_names::Tuple = axes_names(metric_result)
-    metric_result.properties[:axes_names] = collect(
-        parentmodule(metrics).human_readable_name.(
-            _axes_names
-        )
-    )
-    metric_result.properties[:axes_units] = collect(axes_units(_axes_names))
-    return metric_result
-end
-
 """
     slice_results(data::YAXArray; timesteps=(:), species=(:), locations=(:), scenarios=(:))
 

From e24c3056255fd42646172b32b19c02f3e2314264 Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Fri, 27 Sep 2024 18:11:26 +1000
Subject: [PATCH 04/15] Use metadata for spatial cluster plots

---
 ext/AvizExt/outcome_metadata.jl |  8 ++++----
 ext/AvizExt/viz/clustering.jl   | 21 ++++++++++++++-------
 2 files changed, 18 insertions(+), 11 deletions(-)

diff --git a/ext/AvizExt/outcome_metadata.jl b/ext/AvizExt/outcome_metadata.jl
index 21be00383..72b98ce75 100644
--- a/ext/AvizExt/outcome_metadata.jl
+++ b/ext/AvizExt/outcome_metadata.jl
@@ -2,12 +2,12 @@ function outcome_title(outcomes::YAXArray)::String
     return get(outcomes.properties, :metric_name, "")
 end
 
-function outcome_label(outcomes::YAXArray)::String
+function outcome_label(outcomes::YAXArray; label_case=titlecase)::String
     outcome_metadata = outcomes.properties
 
-    return if haskey(outcome_metadata, :metric_feature) &&
-        haskey(outcome_metadata, :metric_unit)
-        "$(outcome_metadata[:metric_feature]) [$(outcome_metadata[:metric_unit])]"
+    return if all(haskey.([outcome_metadata], [:metric_feature, :metric_unit]))
+        _metric_feature = label_case(outcome_metadata[:metric_feature])
+        "$(_metric_feature) [$(outcome_metadata[:metric_unit])]"
     else
         ""
     end
diff --git a/ext/AvizExt/viz/clustering.jl b/ext/AvizExt/viz/clustering.jl
index 9ae471da8..2b0c2b00e 100644
--- a/ext/AvizExt/viz/clustering.jl
+++ b/ext/AvizExt/viz/clustering.jl
@@ -105,7 +105,7 @@ Visualize clustered time series for each location and map.
 
 # Arguments
 - `rs` : ResultSet
-- `data` : Vector of summary statistics data for each location
+- `loc_outcomes` : Vector of summary statistics data for each location
 - `clusters` : Vector of numbers corresponding to clusters
 - `opts` : Options specific to this plotting method
     - `highlight` : Vector of colors indicating cluster membership for each location.
@@ -117,7 +117,7 @@ Figure
 """
 function ADRIA.viz.map(
     rs::Union{Domain,ResultSet},
-    data::AbstractArray{<:Real},
+    loc_outcomes::AbstractVector{<:Real},
     clusters::Union{BitVector,AbstractVector{Int64}};
     opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
     fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
@@ -125,14 +125,21 @@ function ADRIA.viz.map(
 )::Figure
     f = Figure(; fig_opts...)
     g = f[1, 1] = GridLayout()
-    ADRIA.viz.map!(g, rs, data, clusters; opts=opts, axis_opts=axis_opts)
+
+    if !haskey(axis_opts, :title)
+        axis_opts[:title] = "$(outcome_title(loc_outcomes)) Clusters"
+    end
+
+    opts[:colorbar_label] = get(opts, :colorbar_label, outcome_label(loc_outcomes))
+
+    ADRIA.viz.map!(g, rs, loc_outcomes, clusters; opts=opts, axis_opts=axis_opts)
 
     return f
 end
 function ADRIA.viz.map!(
     g::Union{GridLayout,GridPosition},
     rs::Union{Domain,ResultSet},
-    data::AbstractVector{<:Real},
+    loc_outcomes::AbstractVector{<:Real},
     clusters::Union{BitVector,Vector{Int64}};
     opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
     axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
@@ -141,7 +148,7 @@ function ADRIA.viz.map!(
     loc_groups::Dict{Symbol,BitVector} = ADRIA.analysis.scenario_clusters(clusters)
     group_colors::Dict{Symbol,Union{Symbol,RGBA{Float32}}} = colors(loc_groups)
 
-    legend_params::Tuple = _cluster_legend_params(data, loc_groups, group_colors)
+    legend_params::Tuple = _cluster_legend_params(loc_outcomes, loc_groups, group_colors)
 
     _colors::Vector{Union{Symbol,RGBA{Float32}}} = Vector{Union{Symbol,RGBA{Float32}}}(
         undef, length(clusters)
@@ -154,7 +161,7 @@ function ADRIA.viz.map!(
     opts[:highlight] = get(opts, :highlight, _colors)
     opts[:legend_params] = get(opts, :legend_params, legend_params)
 
-    ADRIA.viz.map!(g, rs, data; opts=opts, axis_opts=axis_opts)
+    ADRIA.viz.map!(g, rs, loc_outcomes; opts=opts, axis_opts=axis_opts)
 
     return g
 end
@@ -191,7 +198,7 @@ function _cluster_legend_params(
 
     legend_labels =
         labels(group_keys) .* ": " .* ADRIA.to_scientific.(label_means, digits=2)
-    legend_title = "Cluster mean"
+    legend_title = "Cluster mean $(outcome_label(data; label_case=lowercase))"
 
     return (legend_entries, legend_labels, legend_title)
 end

From 0cb52288ec6482c494f75249bfecb719eedd1df4 Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Fri, 27 Sep 2024 19:19:42 +1000
Subject: [PATCH 05/15] Fix outcome_label

---
 ext/AvizExt/outcome_metadata.jl |  7 ++++++-
 src/metrics/metadata.jl         | 25 ++++++++++++++++++-------
 2 files changed, 24 insertions(+), 8 deletions(-)

diff --git a/ext/AvizExt/outcome_metadata.jl b/ext/AvizExt/outcome_metadata.jl
index 72b98ce75..6b661923f 100644
--- a/ext/AvizExt/outcome_metadata.jl
+++ b/ext/AvizExt/outcome_metadata.jl
@@ -7,7 +7,12 @@ function outcome_label(outcomes::YAXArray; label_case=titlecase)::String
 
     return if all(haskey.([outcome_metadata], [:metric_feature, :metric_unit]))
         _metric_feature = label_case(outcome_metadata[:metric_feature])
-        "$(_metric_feature) [$(outcome_metadata[:metric_unit])]"
+        _metric_label = if !isempty(outcome_metadata[:metric_unit])
+            "[$(outcome_metadata[:metric_unit])]"
+        else
+            ""
+        end
+        "$(_metric_feature) $(_metric_label)"
     else
         ""
     end
diff --git a/src/metrics/metadata.jl b/src/metrics/metadata.jl
index 85e24ad16..94a9c06d7 100644
--- a/src/metrics/metadata.jl
+++ b/src/metrics/metadata.jl
@@ -1,19 +1,25 @@
 const METRIC_METADATA = (:metric_name, :metric_feature, :is_relative, :metric_unit)
 const AXIS_METADATA = (:axes_names, :axes_units)
 
+"""
+    metadata(outcomes::YAXArray)::Dict{Symbol,Any}
+
+Helper function to extract metadata from YAXArrays.
+"""
 function metadata(outcomes::YAXArray)::Dict{Symbol,Any}
     return outcomes.properties
 end
 
 """
-    fill_metadata!(outcomes::YAXArray, metric::Metric)::Nothing
-    fill_metadata!(outcomes::YAXArray; kwargs...)::Nothing
+    fill_metadata!(outcomes::YAXArray{T,N,A}, metric::Metric)::YAXArray{T,N,A} where {T,N,A}
+    fill_metadata!(outcomes::YAXArray{T,N,A}, metadata::Dict{Symbol,Any})::YAXArray{T,N,A} where {T,N,A}
 
-Fill `:axes_names` and `:axes_units` properties of the datacube.
+Fill outcomes YAXArray metadata (`properties` attribute).
 
 # Arguments
-- `outcomes` : YAXArray datacube of metric outcomes
-- `metric` : ADRIA.metrics.Metric object
+- `outcomes` : YAXArray datacube of metric outcomes.
+- `metric` : ADRIA.metrics.Metric object.
+- `metadata` : Dict to be used to fill outcomes metrics metadata.
 """
 function fill_metadata!(
     outcomes::YAXArray{T,N,A}, metric::Metric
@@ -51,6 +57,11 @@ function fill_metadata!(
     return outcomes
 end
 
+"""
+    fill_axes_metadata!(outcomes::YAXArray)::Nothing
+
+Fill outcomes axes metadata.
+"""
 function fill_axes_metadata!(outcomes::YAXArray)::Nothing
     _axes_names::Tuple = axes_names(outcomes)
     outcomes.properties[:axes_names] = collect(
@@ -63,9 +74,9 @@ function fill_axes_metadata!(outcomes::YAXArray)::Nothing
 end
 
 """
-    axes_units(axes_names::Tuple)::Tuple
+    axes_units(axes_names::Union{Vector{Symbol},Tuple})::Tuple
 
-Units for each metric axis.
+Get units for each metric axis.
 """
 function axes_units(axes_names::Union{Vector{Symbol},Tuple})::Tuple
     return values((timesteps="years", species="", locations="", scenarios="")[axes_names])

From 69a5e5bdb937a3479559cd9de0082b573af760fa Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Mon, 30 Sep 2024 13:34:49 +1000
Subject: [PATCH 06/15] Replace taxa by species in viz/taxa_dynamics

---
 ext/AvizExt/viz/taxa_dynamics.jl | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/ext/AvizExt/viz/taxa_dynamics.jl b/ext/AvizExt/viz/taxa_dynamics.jl
index 98b00cc55..d347590a8 100644
--- a/ext/AvizExt/viz/taxa_dynamics.jl
+++ b/ext/AvizExt/viz/taxa_dynamics.jl
@@ -20,12 +20,12 @@ ADRIA.viz.taxonomy(scenarios, relative_taxa_cover)
 # Arguments
 - `rs` : ADRIA result set
 - `scenarios` : Scenario specification
-- `relative_taxa_cover` : YAXArray of dimensions [timesteps ⋅ taxa ⋅ scenarios]
+- `relative_taxa_cover` : YAXArray of dimensions [timesteps ⋅ species ⋅ scenarios]
 - `opts` : Aviz options
     - `by_RCP` : Split plots by RCP otherwise split by scenario type. Defaults to false.
     - `by_functional_groups` : If true, split plots by scenario types, otherwise split by taxonomy. Defaults to true.
     - `show_confints` : Show confidence intervals around series. Defaults to true.
-    - `colors` : Colormap for each taxonomy or scenario type. Defaults to Set1_5 for taxa and ADRIA defaults for scenario type.
+    - `colors` : Colormap for each taxonomy or scenario type. Defaults to Set1_5 for species and ADRIA defaults for scenario type.
 - `axis_opts` : Additional options to pass to adjust Axis attributes.
   See: https://docs.makie.org/v0.19/api/index.html#Axis
 - `series_opts` : Additional options to pass to adjust Series attributes
@@ -103,7 +103,7 @@ function ADRIA.viz.taxonomy!(
     by_functional_groups::Bool = get(opts, :by_functional_groups, true)
     if by_functional_groups
         # Create colors
-        n_functional_groups::Int64 = length(relative_taxa_cover.taxa)
+        n_functional_groups::Int64 = length(relative_taxa_cover.species)
         default_color = Symbol("Set1_" * string(n_functional_groups))
         color = get(opts, :colors, default_color)
         _colors = categorical_colors(color, n_functional_groups)
@@ -193,12 +193,13 @@ function taxonomy_by_intervention!(
     series_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
 )::Nothing where {T<:Float32}
     n_timesteps::Int64 = length(relative_taxa_cover.timesteps)
-    n_functional_groups::Int64 = length(relative_taxa_cover.taxa)
+    functional_groups::Vector{Int64} = ADRIA.axis_labels(relative_taxa_cover, :species)
+    n_functional_groups::Int64 = length(functional_groups)
 
     # Plot and calculate confidence intervals
     confints = zeros(n_timesteps, n_functional_groups, 3)
-    for (idx, taxa) in enumerate(relative_taxa_cover.taxa)
-        confints[:, idx, :] = series_confint(relative_taxa_cover[taxa=At(taxa)])
+    for (idx, group) in enumerate(functional_groups)
+        confints[:, idx, :] = series_confint(relative_taxa_cover[species=At(group)])
         show_confints ?
         band!(
             ax, 1:n_timesteps, confints[:, idx, 1], confints[:, idx, 3];
@@ -247,7 +248,7 @@ function intervention_by_taxonomy!(
 
         intervention_by_taxonomy!(
             ax,
-            relative_taxa_cover[taxa=idx],
+            relative_taxa_cover[species=idx],
             colors,
             scen_groups;
             show_confints=show_confints,

From e0e159940dcf8c4b1d034ef9cf6c7627239a3806 Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Tue, 1 Oct 2024 13:41:45 +1000
Subject: [PATCH 07/15] Add basic tests for viz/taxa_dynamics

---
 test/runtests.jl          |  1 +
 test/viz/taxa_dynamics.jl | 36 ++++++++++++++++++++++++++++++++++++
 2 files changed, 37 insertions(+)
 create mode 100644 test/viz/taxa_dynamics.jl

diff --git a/test/runtests.jl b/test/runtests.jl
index 49df00260..d1713dc82 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,6 +24,7 @@ include("metrics/scenario.jl")
 include("metrics/metrics.jl")
 include("utils/scale.jl")
 include("utils/text_display.jl")
+include("viz/taxa_dynamics.jl")
 
 # TODO Fix spatial_clustering and site_selection tests
 # include("site_selection.jl")
diff --git a/test/viz/taxa_dynamics.jl b/test/viz/taxa_dynamics.jl
new file mode 100644
index 000000000..11fab41c3
--- /dev/null
+++ b/test/viz/taxa_dynamics.jl
@@ -0,0 +1,36 @@
+using WGLMakie, GeoMakie, GraphMakie
+
+using ADRIA
+
+if !@isdefined(TEST_RS)
+    const TEST_DOM, TEST_N_SAMPLES, TEST_SCENS, TEST_RS = test_rs()
+end
+
+@testset "taxonomy" begin
+    @testset "Default opts" begin
+        @test ADRIA.viz.taxonomy(TEST_RS) isa Figure
+    end
+
+    # Default opts values
+    opts::Dict{Symbol,Any} = Dict(
+        :by_functional_groups => true, :by_RCP => false, :show_confints => true
+    )
+
+    @testset "Split by taxonomy" begin
+        opts[:by_functional_groups] = false
+        @test ADRIA.viz.taxonomy(TEST_RS; opts=opts) isa Figure
+    end
+
+    @testset "Split by RCP" begin
+        opts[:by_functional_groups] = true
+        opts[:by_RCP] = true
+        @test ADRIA.viz.taxonomy(TEST_RS; opts=opts) isa Figure
+    end
+
+    @testset "Don't show confidence intervals" begin
+        opts[:by_functional_groups] = true
+        opts[:by_RCP] = true
+        opts[:show_confints] = false
+        @test ADRIA.viz.taxonomy(TEST_RS; opts=opts) isa Figure
+    end
+end

From 938ee0c4b273bbe438915149e1e430621c3b3d5e Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Tue, 1 Oct 2024 13:42:07 +1000
Subject: [PATCH 08/15] Add title for taxa dynamics plot

---
 ext/AvizExt/viz/taxa_dynamics.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/ext/AvizExt/viz/taxa_dynamics.jl b/ext/AvizExt/viz/taxa_dynamics.jl
index d347590a8..3825308d9 100644
--- a/ext/AvizExt/viz/taxa_dynamics.jl
+++ b/ext/AvizExt/viz/taxa_dynamics.jl
@@ -139,6 +139,10 @@ function ADRIA.viz.taxonomy!(
         )
     end
 
+    Label(
+        g[1, :, Top()], "Taxa dynamics"; padding=(0, 0, 30, 0), font=:bold, valign=:bottom
+    )
+
     return g
 end
 

From 3764b1f7d0cabc970f0ef2e34ebb164ba437dfeb Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Tue, 1 Oct 2024 13:55:28 +1000
Subject: [PATCH 09/15] Update x-axis label for scenarios plot

---
 ext/AvizExt/viz/scenarios.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ext/AvizExt/viz/scenarios.jl b/ext/AvizExt/viz/scenarios.jl
index c35fc3372..f83acbcd7 100644
--- a/ext/AvizExt/viz/scenarios.jl
+++ b/ext/AvizExt/viz/scenarios.jl
@@ -184,7 +184,7 @@ function ADRIA.viz.scenarios!(
         _render_legend(g, scen_groups, legend_position, legend_labels)
     end
 
-    ax.xlabel = "Time [years]"
+    ax.xlabel = "Year"
     ax.ylabel = outcome_label(outcomes)
     return g
 end

From 68782e33fcb785430ea9d01f02101e474bb83b74 Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Tue, 1 Oct 2024 16:52:59 +1000
Subject: [PATCH 10/15] Add colorbar label for map plots

Also use set_figure_defaults() for cluster maps to follow the same standard as other map plots.
---
 ext/AvizExt/outcome_metadata.jl | 22 +++++++++++++++++++---
 ext/AvizExt/viz/clustering.jl   | 10 +++-------
 ext/AvizExt/viz/spatial.jl      |  2 ++
 3 files changed, 24 insertions(+), 10 deletions(-)

diff --git a/ext/AvizExt/outcome_metadata.jl b/ext/AvizExt/outcome_metadata.jl
index 6b661923f..672cc158d 100644
--- a/ext/AvizExt/outcome_metadata.jl
+++ b/ext/AvizExt/outcome_metadata.jl
@@ -2,11 +2,27 @@ function outcome_title(outcomes::YAXArray)::String
     return get(outcomes.properties, :metric_name, "")
 end
 
-function outcome_label(outcomes::YAXArray; label_case=titlecase)::String
+function set_plot_opts!(
+    outcomes::AbstractArray,
+    opts::OPT_TYPE,
+    opts_key::Symbol;
+    metadata_key::Symbol=:metric_feature,
+    label_case=titlecase
+)
+    if !haskey(opts, opts_key) && (outcomes isa YAXArray)
+        opts[opts_key] = outcome_label(
+            outcomes; metadata_key=metadata_key, label_case=label_case
+        )
+    end
+end
+
+function outcome_label(
+    outcomes::YAXArray; metadata_key::Symbol=:metric_feature, label_case=titlecase
+)::String
     outcome_metadata = outcomes.properties
 
-    return if all(haskey.([outcome_metadata], [:metric_feature, :metric_unit]))
-        _metric_feature = label_case(outcome_metadata[:metric_feature])
+    return if all(haskey.([outcome_metadata], [metadata_key, :metric_unit]))
+        _metric_feature = label_case(outcome_metadata[metadata_key])
         _metric_label = if !isempty(outcome_metadata[:metric_unit])
             "[$(outcome_metadata[:metric_unit])]"
         else
diff --git a/ext/AvizExt/viz/clustering.jl b/ext/AvizExt/viz/clustering.jl
index 2b0c2b00e..e320feef0 100644
--- a/ext/AvizExt/viz/clustering.jl
+++ b/ext/AvizExt/viz/clustering.jl
@@ -120,17 +120,13 @@ function ADRIA.viz.map(
     loc_outcomes::AbstractVector{<:Real},
     clusters::Union{BitVector,AbstractVector{Int64}};
     opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
-    fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
-    axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
+    fig_opts::OPT_TYPE=set_figure_defaults(DEFAULT_OPT_TYPE()),
+    axis_opts::OPT_TYPE=set_axis_defaults(DEFAULT_OPT_TYPE())
 )::Figure
     f = Figure(; fig_opts...)
     g = f[1, 1] = GridLayout()
 
-    if !haskey(axis_opts, :title)
-        axis_opts[:title] = "$(outcome_title(loc_outcomes)) Clusters"
-    end
-
-    opts[:colorbar_label] = get(opts, :colorbar_label, outcome_label(loc_outcomes))
+    set_plot_opts!(loc_outcomes, opts, :colorbar_label)
 
     ADRIA.viz.map!(g, rs, loc_outcomes, clusters; opts=opts, axis_opts=axis_opts)
 
diff --git a/ext/AvizExt/viz/spatial.jl b/ext/AvizExt/viz/spatial.jl
index 6722e3e9a..cb8ae7840 100644
--- a/ext/AvizExt/viz/spatial.jl
+++ b/ext/AvizExt/viz/spatial.jl
@@ -182,6 +182,8 @@ function ADRIA.viz.map(
     fig_opts::OPT_TYPE=set_figure_defaults(DEFAULT_OPT_TYPE()),
     axis_opts::OPT_TYPE=set_axis_defaults(DEFAULT_OPT_TYPE())
 )
+    set_plot_opts!(y, opts, :colorbar_label; metadata_key=:metric_name)
+
     f = Figure(; fig_opts...)
     g = f[1, 1] = GridLayout()
 

From 184b9ff826c9e1e0c8743341c3f9dee68ab690c1 Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Tue, 1 Oct 2024 17:26:33 +1000
Subject: [PATCH 11/15] Add title to viz/rule_extraction

---
 ext/AvizExt/viz/rule_extraction.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/ext/AvizExt/viz/rule_extraction.jl b/ext/AvizExt/viz/rule_extraction.jl
index aee06d257..acf9d5390 100644
--- a/ext/AvizExt/viz/rule_extraction.jl
+++ b/ext/AvizExt/viz/rule_extraction.jl
@@ -141,6 +141,11 @@ function ADRIA.viz.rules_scatter!(
         margin=(5, 5, 5, 5)
     )
 
+    Label(
+        g[1, :, Top()], "SIRUS extracted rules"; padding=(0, 0, 50, 0), font=:bold,
+        valign=:bottom
+    )
+
     return g
 end
 

From b08615ca29d078ed5a0e5e4b96efb3a43b311ccd Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Tue, 1 Oct 2024 17:27:19 +1000
Subject: [PATCH 12/15] Set default values for viz/rule_extraction figure size

---
 ext/AvizExt/viz/rule_extraction.jl | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/ext/AvizExt/viz/rule_extraction.jl b/ext/AvizExt/viz/rule_extraction.jl
index acf9d5390..ce59fc73e 100644
--- a/ext/AvizExt/viz/rule_extraction.jl
+++ b/ext/AvizExt/viz/rule_extraction.jl
@@ -34,12 +34,18 @@ function ADRIA.viz.rules_scatter(
     fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
     axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
 )::Figure
-    f = Figure(; fig_opts...)
-    g = f[1, 1] = GridLayout()
-
     # For now we only plot conditions with two clauses
     rules = filter(r -> length(r.condition) == 2, rules)
 
+    # Setup figure size
+    n_rows = ceil(10 / 4)
+    base_width = 1200
+    base_height = n_rows * base_width / 5
+    fig_opts[:size] = get(fig_opts, :size, (base_width, base_height))
+
+    f = Figure(; fig_opts...)
+    g = f[1, 1] = GridLayout()
+
     ADRIA.viz.rules_scatter!(
         g, rs, scenarios, clusters, rules; opts=opts, axis_opts=axis_opts
     )

From 7d16dd20effdb7ddf7b1214d4253d6dc89678f33 Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Tue, 1 Oct 2024 17:51:50 +1000
Subject: [PATCH 13/15] Set cluster map default size

---
 ext/AvizExt/viz/clustering.jl | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/ext/AvizExt/viz/clustering.jl b/ext/AvizExt/viz/clustering.jl
index e320feef0..da0833191 100644
--- a/ext/AvizExt/viz/clustering.jl
+++ b/ext/AvizExt/viz/clustering.jl
@@ -120,9 +120,13 @@ function ADRIA.viz.map(
     loc_outcomes::AbstractVector{<:Real},
     clusters::Union{BitVector,AbstractVector{Int64}};
     opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
-    fig_opts::OPT_TYPE=set_figure_defaults(DEFAULT_OPT_TYPE()),
+    fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
     axis_opts::OPT_TYPE=set_axis_defaults(DEFAULT_OPT_TYPE())
 )::Figure
+    # Setup fig_opts size before the other default fig_opts
+    fig_opts[:size] = get(fig_opts, :size, (800, 700))
+    set_figure_defaults(fig_opts)
+
     f = Figure(; fig_opts...)
     g = f[1, 1] = GridLayout()
 

From f5a3005ebfcf486522f0b11380f512e4687edee3 Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Thu, 3 Oct 2024 18:15:04 +1000
Subject: [PATCH 14/15] Update rule extraction function to ignore constant
 factors

---
 src/analysis/rule_extraction.jl | 42 ++++++++++++++++++++++++++-------
 1 file changed, 34 insertions(+), 8 deletions(-)

diff --git a/src/analysis/rule_extraction.jl b/src/analysis/rule_extraction.jl
index 712aed9ba..9cd2d6b33 100644
--- a/src/analysis/rule_extraction.jl
+++ b/src/analysis/rule_extraction.jl
@@ -103,16 +103,18 @@ function print_rules(rules::Vector{Rule{Vector{Vector},Vector{Float64}}})::Nothi
 end
 
 """
-    cluster_rules(clusters::Vector{T}, X::DataFrame, max_rules::T; seed::Int64=123, kwargs...) where {T<:Integer,F<:Real}
-    cluster_rules(clusters::Union{BitVector,Vector{Bool}}, X::DataFrame, max_rules::T; kwargs...) where {T<:Int64}
+    cluster_rules(domain::Domain, clusters::Vector{T}, scenarios::DataFrame, factors::Vector{Symbol}, max_rules::T; seed::Int64=123, kwargs...) where {T<:Integer,F<:Real}
+    cluster_rules(domain::Domain, clusters::Union{BitVector,Vector{Bool}}, scenarios::DataFrame, factors::Vector{Symbol}, max_rules::T; kwargs...) where {T<:Int64}
 
 Use SIRUS package to extract rules from time series clusters based on some summary metric
 (default is median). More information about the keyword arguments accepeted can be found in
 MLJ's doc (https://juliaai.github.io/MLJ.jl/dev/models/StableRulesClassifier_SIRUS/).
 
 # Arguments
+- `domain` : Domain
 - `clusters` : Vector of cluster indexes for each scenario outcome
-- `X` : Features to be used as input by SIRUS
+- `scenarios` : Scenarios DataFrame
+- `factors` : Vector of factors of interest
 - `max_rules` : Maximum number of rules, to be used as input by SIRUS
 - `seed` : Seed to be used by RGN
 
@@ -129,8 +131,24 @@ A StableRules object (implemented by SIRUS).
    Electron. J. Statist. 15 (1) 427 - 505.
    https://doi.org//10.1214/20-EJS1792
 """
-function cluster_rules(clusters::Vector{T}, X::DataFrame, max_rules::T;
-    seed::Int64=123, kwargs...) where {T<:Int64}
+function cluster_rules(
+    domain::ADRIA.Domain,
+    clusters::Vector{T},
+    scenarios::DataFrame,
+    factors::Vector{Symbol},
+    max_rules::T;
+    seed::Int64=123,
+    kwargs...
+) where {T<:Int64}
+    ms = ADRIA.model_spec(domain)
+    variable_factors_filter::BitVector = .!ms[ms.fieldname .∈ [factors], :is_constant]
+    variable_factors::Vector{Symbol} = factors[variable_factors_filter]
+
+    if isempty(variable_factors)
+        throw(ArgumentError("Factors of interest have to be non constant."))
+    end
+
+    X = scenarios[:, variable_factors]
     # Set seed and Random Number Generator
     rng = StableRNG(seed)
 
@@ -150,9 +168,17 @@ function cluster_rules(clusters::Vector{T}, X::DataFrame, max_rules::T;
 
     return rules(mach.fitresult)
 end
-function cluster_rules(clusters::Union{BitVector,Vector{Bool}}, X::DataFrame, max_rules::T;
-    kwargs...) where {T<:Int64}
-    return cluster_rules(convert.(Int64, clusters), X, max_rules; kwargs...)
+function cluster_rules(
+    domain::ADRIA.Domain,
+    clusters::Union{BitVector,Vector{Bool}},
+    scenarios::DataFrame,
+    factors::Vector{Symbol},
+    max_rules::T;
+    kwargs...
+) where {T<:Int64}
+    return cluster_rules(
+        domain, convert.(Int64, clusters), scenarios, factors, max_rules; kwargs...
+    )
 end
 
 """

From 263a125c90d5ea6af40b214243332d8a9ed14824 Mon Sep 17 00:00:00 2001
From: Pedro Ribeiro de Almeida <pedroribalmeida@gmail.com>
Date: Fri, 4 Oct 2024 12:22:17 +1000
Subject: [PATCH 15/15] Update rule extraction documentation

---
 docs/src/usage/analysis.md | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/src/usage/analysis.md b/docs/src/usage/analysis.md
index 9d47aa63e..308908f55 100644
--- a/docs/src/usage/analysis.md
+++ b/docs/src/usage/analysis.md
@@ -405,18 +405,18 @@ scenario, it is possible to use a Rule Induction algorithm (SIRUS) and plot each
 rule as a scatter graph:
 
 ```julia
-# Select only desired features
-fields_iv = ADRIA.component_params(rs, [Intervention, CriteriaWeights]).fieldname
-scenarios_iv = scens[:, fields_iv]
+# Select features of interest
+foi = ADRIA.component_params(rs, [Intervention, SeedCriteriaWeights]).fieldname
 
 # Use SIRUS algorithm to extract rules
 max_rules = 10
-rules_iv = ADRIA.analysis.cluster_rules(target_clusters, scenarios_iv, max_rules)
+rules_iv = ADRIA.analysis.cluster_rules(dom, target_clusters, scens, foi, max_rules)
+
 
 # Plot scatters for each rule highlighting the area selected them
 rules_scatter_fig = ADRIA.viz.rules_scatter(
     rs,
-    scenarios_iv,
+    scens,
     target_clusters,
     rules_iv;
     fig_opts=fig_opts,