diff --git a/ext/AvizExt/viz/environment/dhw.jl b/ext/AvizExt/viz/environment/dhw.jl
index 85eadb476..70241aaaf 100644
--- a/ext/AvizExt/viz/environment/dhw.jl
+++ b/ext/AvizExt/viz/environment/dhw.jl
@@ -1,8 +1,29 @@
 using Statistics
 
+"""
+    ADRIA.viz.dhw_scenario(
+        dom::Domain,
+        scen_id::Int64;
+        fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
+        axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
+    )
+
+Display the trajectory of a single DHW scenario.
+
+# Arguments
+- `dom` : ADRIA Domain
+- `scen_id` : ID of DHW scenario to display
+- `fig_opts` : Additional options to pass to adjust Figure creation
+- `axis_opts` : Additional options to pass to adjust Axis attributes
+
+# Returns
+Figure
+"""
 function ADRIA.viz.dhw_scenario(
-    dom::Domain, scen_id::Int64;
-    fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(), axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
+    dom::Domain,
+    scen_id::Int64;
+    fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
+    axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
 )
     loc_scens = dom.dhw_scens[:, :, scen_id]
     mean_dhw_scen = dropdims(mean(loc_scens; dims=2); dims=2)
@@ -28,3 +49,75 @@ function ADRIA.viz.dhw_scenario(
 
     return f
 end
+
+"""
+    dhw_scenarios(
+        dom::Domain;
+        ci_level=0.95,
+        fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
+        axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
+    )
+
+Display all DHW scenarios for the currently set RCP/SSP.
+
+The set RCP/SSP can be changed with `switch_RCPs!()`.
+
+# Arguments
+- `dom` : ADRIA Domain
+- `ci_level` : Level of Confidence Interval
+- `fig_opts` : Additional options to pass to adjust Figure creation
+- `axis_opts` : Additional options to pass to adjust Axis attributes
+
+# Returns
+Figure
+"""
+function ADRIA.viz.dhw_scenarios(
+    dom::Domain;
+    ci_level=0.95,
+    fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
+    axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE()
+)
+    data = dom.dhw_scens
+
+    # Extract dimension information
+    years = dom.env_layer_md.timeframe
+
+    fig_size = get(fig_opts, :size, (1200, 600))
+
+    axis_opts[:title] = get(axis_opts, :title, "DHW Scenarios")
+    axis_opts[:xlabel] = get(axis_opts, :xlabel, "Year")
+    axis_opts[:ylabel] = get(axis_opts, :ylabel, "DHW")
+    axis_opts[:xticklabelrotation] = get(axis_opts, :xticklabelrotation, 2 / π)
+    fig = Figure(; size=fig_size)
+    ax = Axis(fig[1, 1]; axis_opts...)
+
+    # Determine CI bounds
+    α = (1 - ci_level) / 2
+
+    means = zeros(length(years))
+    lower_ci = zeros(length(years))
+    upper_ci = zeros(length(years))
+
+    for t in collect(axes(data, 1))
+        # Calculate statistics across scenarios for each timestep
+        means[t] = mean(data[t, :, :])
+        lower_ci[t] = quantile(data[t, :, :][:], α)
+        upper_ci[t] = quantile(data[t, :, :][:], 1 - α)
+    end
+
+    ci_band = band!(ax, years, lower_ci, upper_ci; color=(:red, 0.3))
+    each_scen = series!(
+        ax, years, dropdims(mean(data; dims=:sites); dims=:sites).data';
+        solid_color=(:red, 0.2)
+    )
+    mean_line = lines!(ax, years, means; color=(:black, 0.7))
+
+    # Add legend
+    Legend(
+        fig[1, 2],
+        [ci_band, each_scen, mean_line],
+        ["$(floor(Int64, ci_level*100))% CI", "Individual trajectories", "Mean"]
+    )
+
+    return fig
+end
diff --git a/src/analysis/feature_set.jl b/src/analysis/feature_set.jl
new file mode 100644
index 000000000..cd2117a3c
--- /dev/null
+++ b/src/analysis/feature_set.jl
@@ -0,0 +1,109 @@
+"""
+    _filter_constants(scens::DataFrame)::DataFrame
+
+Filter out features/factors that do not vary.
+"""
+function _filter_constants(scens::DataFrame)::DataFrame
+    varying_cols = [
+        i for (i, col) in enumerate(eachcol(scens))
+        if !all(==(first(col)), col)
+    ]
+
+    return scens[:, varying_cols]
+end
+
+"""
+    _seeded(rs::ResultSet)::Tuple{Vector}
+
+Extract total and average deployment, currently in terms of proportional increase to cover
+relative to the locations' carrying capacity.
+
+# Arguments
+- `rs` : ResultSet holding scenario outcomes
+
+# Returns
+DataFrames of mean and total deployment for each coral group
+"""
+function _seeded(rs::ResultSet)::Tuple{DataFrame,DataFrame}
+    deployed_corals = rs.seed_log.coral_id
+
+    mean_deployment = [
+        mean(sum(rs.seed_log[:, c_id, :, :]; dims=:timesteps); dims=:locations).data[:]
+        for c_id in deployed_corals
+    ]
+
+    total_deployment = [
+        sum(sum(rs.seed_log[:, c_id, :, :]; dims=:timesteps); dims=:locations).data[:]
+        for c_id in deployed_corals
+    ]
+
+    col_names = ["deployed_coral_$(i)" for i in string.(collect(deployed_corals))]
+    μ = DataFrame(hcat(mean_deployment...), "mean_" .* col_names)
+    T = DataFrame(hcat(total_deployment...), "total_" .* col_names)
+
+    return μ, T
+end
+
+"""
+    feature_set(rs::ResultSet)::DataFrame
+
+Extract a feature set from results for analysis purposes.
+"""
+function feature_set(rs::ResultSet)::DataFrame
+    scens = copy(rs.inputs)
+
+    rcp_ids = collect(keys(rs.dhw_stats))
+    rcp_id = first(rcp_ids)
+    if length(rcp_ids) > 1
+        @warn "Multiple RCPs found, assigning stats for first id: $(rcp_id)"
+    end
+
+    # Add DHW statistics
+    dhw_stat = mean(rs.dhw_stats[rcp_id]; dims=:locations)
+    dhw_means = dhw_stat[stat=At("mean")].data[:]
+    dhw_stdevs = dhw_stat[stat=At("std")].data[:]
+
+    insertcols!(scens, 2, :dhw_mean => -99.0, :dhw_stdev => -99.0)
+    for (i, r) in enumerate(eachrow(scens))
+        scens[i, :dhw_mean] = dhw_means[Int64(r.dhw_scenario)]
+        scens[i, :dhw_stdev] = dhw_stdevs[Int64(r.dhw_scenario)]
+    end
+
+    @assert all(scens.dhw_mean .> -99.0) "Unknown DHW scenario found, check rows: $(findall(scens.dhw_mean .== -99.0))"
+    @assert all(scens.dhw_stdev .> -99.0) "Unknown DHW scenario found, check rows: $(findall(scens.dhw_mean .== -99.0))"
+
+    # Add indicators of deployments
+    seed_stats = ADRIA.decision.deployment_summary_stats(rs.ranks, :seed)
+    fog_stats = ADRIA.decision.deployment_summary_stats(rs.ranks, :fog)
+
+    # Only attach mean of deployment effort
+    insertcols!(
+        scens,
+        :n_loc_seed_mean => seed_stats[stats=At(:mean)].data[:],
+        :n_loc_fog_mean => fog_stats[stats=At(:mean)].data[:]
+    )
+
+    # Replace `depth_offset` with maximum depth
+    scens.depth_max = scens.depth_min .+ scens.depth_offset
+    scens = scens[:, Not(:depth_offset)]
+
+    total_seeded, mean_seeded = _seeded(rs)
+    DataFrames.hcat!(scens, mean_seeded)
+    DataFrames.hcat!(scens, total_seeded)
+    scens.total_deployed_coral = sum.(eachrow(total_seeded))
+
+    # Remove `dhw_scenario` as Scenario IDs are not very informative for analyses
+    scens = scens[:, Not(:dhw_scenario)]
+
+    # Remove correlated features
+    # Remove desired seed deployment targets
+    # N_seed_TA, etc, indicate maximum deployment effort, not actual simulated deployment
+    scens = scens[:, Not([:N_seed_TA, :N_seed_CA, :N_seed_SM])]
+
+    # Set missing values to 0
+    for col in eachcol(scens)
+        replace!(col, missing => 0.0)
+    end
+
+    return _filter_constants(scens)
+end
diff --git a/src/viz/viz.jl b/src/viz/viz.jl
index 79ce897e6..412a50b6f 100644
--- a/src/viz/viz.jl
+++ b/src/viz/viz.jl
@@ -47,6 +47,7 @@ function connectivity!() end
 
 # Environment
 function dhw_scenario() end
+function dhw_scenarios() end
 function cyclone_scenario() end
 
 # Coral Dynamics