From 2d5efd44e89a711817132bc92a21b49fafb4a26c Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Sat, 4 Jan 2025 23:03:07 +1100 Subject: [PATCH 01/13] Initial makie migration --- Project.toml | 6 + ext/MakieExt/MakieExt.jl | 380 +++++++++++++++++++++++++++++++++++++++ src/Streamfall.jl | 1 + src/viz/viz.jl | 7 + 4 files changed, 394 insertions(+) create mode 100644 ext/MakieExt/MakieExt.jl create mode 100644 src/viz/viz.jl diff --git a/Project.toml b/Project.toml index 6555a2d5..8ed4a3e3 100644 --- a/Project.toml +++ b/Project.toml @@ -42,6 +42,12 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" +[weakdeps] +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" + +[extensions] +MakieExt = "Makie" + [compat] BlackBoxOptim = "0.5, 0.6" Bootstrap = "2" diff --git a/ext/MakieExt/MakieExt.jl b/ext/MakieExt/MakieExt.jl new file mode 100644 index 00000000..d027b4e1 --- /dev/null +++ b/ext/MakieExt/MakieExt.jl @@ -0,0 +1,380 @@ +module MakieExt + +using Statistics +using DataFrames +using Dates + +using Streamfall +import Streamfall: TemporalCrossSection + +using Makie + + +function Streamfall.viz.quickplot(node::NetworkNode) + fig = lines(node.outflow) + return fig +end + +function Streamfall.viz.quickplot(node::NetworkNode, climate::Climate) + date = timesteps(climate) + + @assert length(date) == length(node.outflow) || "Date length and result lengths do not match!" + + f, ax, sp = lines(date, node.outflow) + + return f +end + +function Streamfall.viz.quickplot( + obs::Vector, node::NetworkNode, climate::Climate; + label="Modeled", log=false, burn_in=1, limit=nothing, metric=Streamfall.mKGE +) + return quickplot(obs, node.outflow, climate; label, log, burn_in=burn_in, limit=limit, metric=metric) +end + +function Streamfall.viz.quickplot( + obs::DataFrame, node::NetworkNode, climate::Climate; + label="", log=false, burn_in=1, limit=nothing, metric=Streamfall.mKGE +) + return Streamfall.viz.quickplot(obs[:, node.name], node.outflow, climate; label, log, burn_in, limit, metric) +end + +function Streamfall.viz.quickplot( + obs::Vector, sim::Vector, climate::Climate; + label="Modeled", log=false, burn_in=1, limit=nothing, metric=Streamfall.mKGE +) + date = timesteps(climate) + last_e = !isnothing(limit) ? limit : lastindex(obs) + show_range = burn_in:last_e + return Streamfall.viz.quickplot(obs[show_range], sim[show_range], date[show_range], label, log; metric=metric) +end +function Streamfall.viz.quickplot(obs::Vector, sim::Vector, xticklabels::Vector, label="Modeled", log=false; metric=Streamfall.mKGE) + @assert length(xticklabels) == length(obs) || "x-axis tick label length and observed lengths do not match!" + @assert length(xticklabels) == length(sim) || "x-axis tick label length and simulated lengths do not match!" + + score = round(metric(obs, sim), digits=4) + metric_name = String(Symbol(metric)) + + if log + # Add small constant in case of 0-flow + obs = copy(obs) + sim = copy(sim) + obs[obs .== 0.0] .+= 1e-4 + sim[sim .== 0.0] .+= 1e-4 + end + + scale = log == false ? identity : log10 + + f = Figure(size=(850,400)) + flow_ax = Axis(f[1,1]; xlabel="Date", ylabel="Streamflow", yscale=scale) + qq_ax = Axis(f[1,2]; xlabel="Observed", ylabel="Modeled", xscale=scale, yscale=scale) + + label = "$(label) ($(metric_name): $(score))" + lines!(flow_ax, xticklabels, obs, label="Observed") + lines!(flow_ax, xticklabels, sim; label=label, alpha=0.5) + + qqplot!( + qq_ax, obs, sim; + qqline=:identity, strokewidth=0.01, strokecolor=(:blue, 0.01), # legend=false, strokewidth=0.03, strokealpha=0.1, + markercolor=(:blue, 0.02) + ) + leg = Legend(f[2, 1:2], flow_ax; orientation=:horizontal) + + return f +end + + +using LaTeXStrings + +""" + temporal_cross_section( + dates, obs; + title="", + ylabel="Mean", + label=nothing, + period=monthday, + kwargs... + ) + +Provides indication of temporal variation and uncertainty across time, grouped by `period`. + +Notes: +Assumes daily data. +Filters out leap days. + +# Arguments +- `dates` : Date of each observation +- `obs` : observed data +- `title` : Plot title +- `ylabel` : Optional replacement ylabel. Uses name of `func` if not provided. +- `period` : Method from `Dates` package to group (defaults to `monthday`) +- `kwargs` : Additional plotting keyword arguments +""" +function Streamfall.viz.temporal_cross_section( + dates, obs; + title="", ylabel="Mean", label=nothing, + period::Function=monthday, + kwargs... +) + if isnothing(label) + label = ylabel + end + + arg_keys = keys(kwargs) + format_func = y -> y + logscale = [:log, :log10] + tmp = nothing + + xsect_res = TemporalCrossSection(dates, obs, period) + + # TODO: This is currently broken + if :yscale in arg_keys || :yaxis in arg_keys + tmp = (:yscale in arg_keys) ? kwargs[:yscale] : kwargs[:yaxis] + + if tmp in logscale + log_obs = symlog(copy(obs)) + format_func = y -> (y != 0) ? L"%$(Int(round(sign(y)) * 10))^{%$(round(abs(y), digits=1))}" : L"0" + log_xsect_res = TemporalCrossSection(dates, log_obs, period) + target = log_xsect_res.cross_section + else + target = xsect_res.cross_section + end + else + target = xsect_res.cross_section + end + + x_section = target.mean + lower_75 = target.lower_75 + upper_75 = target.upper_75 + lower_95 = target.lower_95 + upper_95 = target.upper_95 + + sp = target.subperiod + xlabels = join.(sp, "-") + + # Calculate statistics for legend + xs_mean = xsect_res.cross_section.mean + m_ind = round(mean(xs_mean), digits=2) + sd_ind = round(std(xs_mean), digits=2) + wr75_m_ind = round(xsect_res.mean_75, digits=2) + wr75_sd_ind = round(xsect_res.std_75, digits=2) + wr95_m_ind = round(xsect_res.mean_95, digits=2) + wr95_sd_ind = round(xsect_res.std_95, digits=2) + + # Create figure + fig = Figure(size=(800, 600)) + ax = Axis( + fig[1, 1], + xlabel=string(nameof(period)), + ylabel=ylabel, + title=title + ) + + # Plot confidence intervals and median line + x = 1:length(xlabels) + + # 95% CI + b95 = band!( + ax, x, lower_95, upper_95, + color=(:lightblue, 0.3) + ) + + # 75% CI + b75 = band!( + ax, x, lower_75, upper_75, + color=(:lightblue, 0.5) + ) + + # Median line + med_line = lines!( + ax, x, x_section, + color=:black + ) + + # Set x-axis labels + if period == monthday + ax.xticks = (1:14:length(xlabels), xlabels[1:14:end]) + elseif period == yearmonth + ax.xticks = (1:12:length(xlabels), xlabels[1:12:end]) + end + + ax.xticklabelrotation = π/4 + + # Apply log scale if specified + if !isnothing(tmp) && (tmp in logscale) + ax.yscale = log10 + ax.formatter = format_func + end + + # Add legend at bottom + Legend( + fig[2, 1], + [med_line, b75, b95], + [ + "Mean of $(label) μ: $(m_ind), σ: $(sd_ind)", + "CI₇₅ μ: $(wr75_m_ind), σ: $(wr75_sd_ind)", + "CI₉₅ μ: $(wr95_m_ind), σ: $(wr95_sd_ind)" + ], + orientation=:horizontal, + tellwidth=false, + tellheight=true + ) + + # Adjust layout + rowsize!(fig.layout, 1, Relative(0.9)) + rowsize!(fig.layout, 2, Relative(0.1)) + + return fig +end + +""" + temporal_cross_section( + dates, obs, sim; + title="", ylabel="Median Error", label=nothing, period::Function=monthday, kwargs... + ) + +Provides indication of temporal variation and uncertainty across time, grouped by `period`. + +Notes: +Assumes daily data. +Filters out leap days. + +# Arguments +- `dates` : Date of each observation +- `obs` : Observed data +- `sim` : Simulated data +- `title` : Optional plot title. Blank if not provided. +- `ylabel` : Optional replacement ylabel. Uses name of `period` if not provided. +- `label` : Optional legend label. Uses `ylabel` if not provided. +- `period` : Method from `Dates` package to group (defaults to `month`) +""" +function Streamfall.viz.temporal_cross_section( + dates, obs, sim; + title="", ylabel="Median Error", label=nothing, period::Function=monthday, kwargs... +) + if isnothing(label) + label = ylabel + end + + arg_keys = keys(kwargs) + format_func = y -> y + logscale = [:log, :log10] + tmp = nothing + + xsect_res = TemporalCrossSection(dates, obs, sim, period) + target = xsect_res.cross_section + + # TODO: This is currently broken + if :yscale in arg_keys || :yaxis in arg_keys + tmp = (:yscale in arg_keys) ? kwargs[:yscale] : kwargs[:yaxis] + + if tmp in logscale + log_obs = symlog(copy(obs)) + log_sim = symlog(copy(sim)) + + # Format function for y-axis tick labels (e.g., 10^x) + format_func = y -> (y != 0) ? L"%$(Int(round(sign(y)) * 10))^{%$(round(abs(y), digits=1))}" : L"0" + + log_xsect_res = TemporalCrossSection(dates, log_obs, log_sim, period) + target = log_xsect_res.cross_section + end + end + + x_section = target.median + lower_75 = target.lower_75 + upper_75 = target.upper_75 + lower_95 = target.lower_95 + upper_95 = target.upper_95 + + sp = target.subperiod + xlabels = join.(sp, "-") + + xsect_df = xsect_res.cross_section + + # Calculate statistics for legend + med = xsect_df.median + m_ind = round(mean(med), digits=2) + sd_ind = round(std(med), digits=2) + wr75_m_ind = round(xsect_res.mean_75, digits=2) + wr75_sd_ind = round(xsect_res.std_75, digits=2) + wr95_m_ind = round(xsect_res.mean_95, digits=2) + wr95_sd_ind = round(xsect_res.std_95, digits=2) + + x = 1:length(xlabels) + + # Create figure + fig = Figure(size=(800, 600)) + ax = Axis( + fig[1, 1], + xlabel=string(nameof(period)), + ylabel=ylabel, + title=title + ) + + # Plot confidence intervals and median line + # 95% CI + b95 = band!( + ax, x, lower_95, upper_95, + color=(:lightblue, 0.3), + label="CI₉₅ μ: $(wr95_m_ind), σ: $(wr95_sd_ind)" + ) + + # 75% CI + b75 = band!( + ax, x, lower_75, upper_75, + color=(:lightblue, 0.5), + label="CI₇₅ μ: $(wr75_m_ind), σ: $(wr75_sd_ind)" + ) + + # Median line + med_line = lines!( + ax, x, x_section, + color=:black, + label="$(label) μ: $(m_ind), σ: $(sd_ind)" + ) + + # Set x-axis labels + if period == monthday + ax.xticks = (1:14:length(xlabels), xlabels[1:14:end]) + elseif period == yearmonth + ax.xticks = (1:12:length(xlabels), xlabels[1:12:end]) + end + + ax.xticklabelrotation = π/4 + + # Apply log scale if specified + if !isnothing(tmp) && (tmp in logscale) + ax.yscale = log10 + ax.formatter = format_func + end + + # Add legend at bottom + Legend( + fig[2, 1], + [med_line, b75, b95], + [ + "$(label) μ: $(m_ind), σ: $(sd_ind)", + "CI₇₅ μ: $(wr75_m_ind), σ: $(wr75_sd_ind)", + "CI₉₅ μ: $(wr95_m_ind), σ: $(wr95_sd_ind)" + ], + orientation=:horizontal, + tellwidth=false, + tellheight=true, + framevisible=false + ) # Equivalent to transparent background + + # Adjust layout + rowsize!(fig.layout, 1, Relative(0.9)) + rowsize!(fig.layout, 2, Relative(0.1)) + + # Process any additional keyword arguments + for (key, value) in kwargs + if key ∉ [:yscale, :yaxis] # Skip already processed arguments + setproperty!(ax, key, value) + end + end + + return fig +end + +end \ No newline at end of file diff --git a/src/Streamfall.jl b/src/Streamfall.jl index f4aea096..182dc1a5 100644 --- a/src/Streamfall.jl +++ b/src/Streamfall.jl @@ -275,6 +275,7 @@ end include("Analysis/Analysis.jl") include("plotting.jl") +include("viz/viz.jl") # Nodes export NetworkNode, GenericNode, GenericDirectNode diff --git a/src/viz/viz.jl b/src/viz/viz.jl new file mode 100644 index 00000000..748e5533 --- /dev/null +++ b/src/viz/viz.jl @@ -0,0 +1,7 @@ +module viz + +function quickplot end + +function temporal_cross_section end + +end \ No newline at end of file From 7b88b162d4f25fe159a259edefc36dc634f0d758 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Sun, 5 Jan 2025 11:00:16 +1100 Subject: [PATCH 02/13] Change case of module name `viz` -> `Viz` Also export all methods --- ext/MakieExt/MakieExt.jl | 20 ++++++++++---------- src/viz/{viz.jl => Viz.jl} | 4 +++- 2 files changed, 13 insertions(+), 11 deletions(-) rename src/viz/{viz.jl => Viz.jl} (55%) diff --git a/ext/MakieExt/MakieExt.jl b/ext/MakieExt/MakieExt.jl index d027b4e1..cdea68ee 100644 --- a/ext/MakieExt/MakieExt.jl +++ b/ext/MakieExt/MakieExt.jl @@ -10,12 +10,12 @@ import Streamfall: TemporalCrossSection using Makie -function Streamfall.viz.quickplot(node::NetworkNode) +function Streamfall.Viz.quickplot(node::NetworkNode) fig = lines(node.outflow) return fig end -function Streamfall.viz.quickplot(node::NetworkNode, climate::Climate) +function Streamfall.Viz.quickplot(node::NetworkNode, climate::Climate) date = timesteps(climate) @assert length(date) == length(node.outflow) || "Date length and result lengths do not match!" @@ -25,30 +25,30 @@ function Streamfall.viz.quickplot(node::NetworkNode, climate::Climate) return f end -function Streamfall.viz.quickplot( +function Streamfall.Viz.quickplot( obs::Vector, node::NetworkNode, climate::Climate; label="Modeled", log=false, burn_in=1, limit=nothing, metric=Streamfall.mKGE ) return quickplot(obs, node.outflow, climate; label, log, burn_in=burn_in, limit=limit, metric=metric) end -function Streamfall.viz.quickplot( +function Streamfall.Viz.quickplot( obs::DataFrame, node::NetworkNode, climate::Climate; label="", log=false, burn_in=1, limit=nothing, metric=Streamfall.mKGE ) - return Streamfall.viz.quickplot(obs[:, node.name], node.outflow, climate; label, log, burn_in, limit, metric) + return Streamfall.Viz.quickplot(obs[:, node.name], node.outflow, climate; label, log, burn_in, limit, metric) end -function Streamfall.viz.quickplot( +function Streamfall.Viz.quickplot( obs::Vector, sim::Vector, climate::Climate; label="Modeled", log=false, burn_in=1, limit=nothing, metric=Streamfall.mKGE ) date = timesteps(climate) last_e = !isnothing(limit) ? limit : lastindex(obs) show_range = burn_in:last_e - return Streamfall.viz.quickplot(obs[show_range], sim[show_range], date[show_range], label, log; metric=metric) + return Streamfall.Viz.quickplot(obs[show_range], sim[show_range], date[show_range], label, log; metric=metric) end -function Streamfall.viz.quickplot(obs::Vector, sim::Vector, xticklabels::Vector, label="Modeled", log=false; metric=Streamfall.mKGE) +function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, xticklabels::Vector, label="Modeled", log=false; metric=Streamfall.mKGE) @assert length(xticklabels) == length(obs) || "x-axis tick label length and observed lengths do not match!" @assert length(xticklabels) == length(sim) || "x-axis tick label length and simulated lengths do not match!" @@ -110,7 +110,7 @@ Filters out leap days. - `period` : Method from `Dates` package to group (defaults to `monthday`) - `kwargs` : Additional plotting keyword arguments """ -function Streamfall.viz.temporal_cross_section( +function Streamfall.Viz.temporal_cross_section( dates, obs; title="", ylabel="Mean", label=nothing, period::Function=monthday, @@ -248,7 +248,7 @@ Filters out leap days. - `label` : Optional legend label. Uses `ylabel` if not provided. - `period` : Method from `Dates` package to group (defaults to `month`) """ -function Streamfall.viz.temporal_cross_section( +function Streamfall.Viz.temporal_cross_section( dates, obs, sim; title="", ylabel="Median Error", label=nothing, period::Function=monthday, kwargs... ) diff --git a/src/viz/viz.jl b/src/viz/Viz.jl similarity index 55% rename from src/viz/viz.jl rename to src/viz/Viz.jl index 748e5533..78175195 100644 --- a/src/viz/viz.jl +++ b/src/viz/Viz.jl @@ -1,7 +1,9 @@ -module viz +module Viz function quickplot end function temporal_cross_section end +export quickplot, temporal_cross_section + end \ No newline at end of file From d23c997a1807d33d7a961ccfb0f103a20edc6bcc Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Sun, 5 Jan 2025 11:00:37 +1100 Subject: [PATCH 03/13] Initial switchover to GraphMakie --- Project.toml | 2 ++ ext/GraphMakieExt.jl/GraphMakieExt.jl | 16 ++++++++++++++++ src/Streamfall.jl | 3 ++- src/viz/NetworkViz.jl | 10 ++++++++++ 4 files changed, 30 insertions(+), 1 deletion(-) create mode 100644 ext/GraphMakieExt.jl/GraphMakieExt.jl create mode 100644 src/viz/NetworkViz.jl diff --git a/Project.toml b/Project.toml index 8ed4a3e3..4ce6a4bd 100644 --- a/Project.toml +++ b/Project.toml @@ -44,9 +44,11 @@ ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [weakdeps] Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" +GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2" [extensions] MakieExt = "Makie" +GraphMakieExt = ["Makie", "GraphMakie"] [compat] BlackBoxOptim = "0.5, 0.6" diff --git a/ext/GraphMakieExt.jl/GraphMakieExt.jl b/ext/GraphMakieExt.jl/GraphMakieExt.jl new file mode 100644 index 00000000..13e5aa91 --- /dev/null +++ b/ext/GraphMakieExt.jl/GraphMakieExt.jl @@ -0,0 +1,16 @@ +module MakieExt + +using Statistics +using DataFrames +using Dates + +using Streamfall +import Streamfall: TemporalCrossSection + +using Makie + +function Streamfall.NetworkViz.plot(sn::StreamfallNetwork) + return graphplot(mg) +end + +end \ No newline at end of file diff --git a/src/Streamfall.jl b/src/Streamfall.jl index 182dc1a5..e471bd1d 100644 --- a/src/Streamfall.jl +++ b/src/Streamfall.jl @@ -275,7 +275,8 @@ end include("Analysis/Analysis.jl") include("plotting.jl") -include("viz/viz.jl") +include("viz/Viz.jl") +include("viz/NetworkViz.jl") # Nodes export NetworkNode, GenericNode, GenericDirectNode diff --git a/src/viz/NetworkViz.jl b/src/viz/NetworkViz.jl new file mode 100644 index 00000000..e478fd5f --- /dev/null +++ b/src/viz/NetworkViz.jl @@ -0,0 +1,10 @@ +module NetworkViz + +using GraphMakie +using Streamfall + +function plot end + +export plot + +end From c618a15711d0375208976a8d768920eff78ae718 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Sun, 5 Jan 2025 12:58:31 +1100 Subject: [PATCH 04/13] Complete initial switch over to Makie --- examples/run_nodes.jl | 1 + ext/GraphMakieExt.jl/GraphMakieExt.jl | 16 ----------- ext/GraphMakieExt/GraphMakieExt.jl | 40 ++++++++++++++++++++++++++ ext/MakieExt/MakieExt.jl | 41 ++++++++++++++++++++++++--- src/Streamfall.jl | 7 +++-- src/viz/NetworkViz.jl | 7 ++--- src/viz/Viz.jl | 5 +++- 7 files changed, 89 insertions(+), 28 deletions(-) delete mode 100644 ext/GraphMakieExt.jl/GraphMakieExt.jl create mode 100644 ext/GraphMakieExt/GraphMakieExt.jl diff --git a/examples/run_nodes.jl b/examples/run_nodes.jl index eb23ff3f..7cfd7206 100644 --- a/examples/run_nodes.jl +++ b/examples/run_nodes.jl @@ -1,3 +1,4 @@ +using GLMakie, GraphMakie using YAML, DataFrames, CSV, Plots using Statistics using Streamfall diff --git a/ext/GraphMakieExt.jl/GraphMakieExt.jl b/ext/GraphMakieExt.jl/GraphMakieExt.jl deleted file mode 100644 index 13e5aa91..00000000 --- a/ext/GraphMakieExt.jl/GraphMakieExt.jl +++ /dev/null @@ -1,16 +0,0 @@ -module MakieExt - -using Statistics -using DataFrames -using Dates - -using Streamfall -import Streamfall: TemporalCrossSection - -using Makie - -function Streamfall.NetworkViz.plot(sn::StreamfallNetwork) - return graphplot(mg) -end - -end \ No newline at end of file diff --git a/ext/GraphMakieExt/GraphMakieExt.jl b/ext/GraphMakieExt/GraphMakieExt.jl new file mode 100644 index 00000000..42b0d113 --- /dev/null +++ b/ext/GraphMakieExt/GraphMakieExt.jl @@ -0,0 +1,40 @@ +module GraphMakieExt + +using Streamfall +import Streamfall: StreamfallNetwork + +using Makie + +using Graphs +using GraphMakie + + +""" + plot_network(sn::StreamfallNetwork) + +Simple plot of stream network. +""" +function Streamfall.NetworkViz.plot_network(sn::StreamfallNetwork) + node_labels = ["$(sn[i].name)\n"*string(nameof(typeof(sn[i]))) for i in vertices(sn.mg)] + f, ax, sp = graphplot( + sn.mg; + nlabels=node_labels, + nlabels_align=(:center, :bottom), + node_size=16, + node_color=:blue + ) + # f.scene.padding = (50, 50, 50, 50) # (left, right, bottom, top) + + # Turn off grid lines + ax.xgridvisible = false + ax.ygridvisible = false + + # Hide ticks and labels + hidedecorations!(ax) + + reset_limits!(ax) + + return f +end + +end \ No newline at end of file diff --git a/ext/MakieExt/MakieExt.jl b/ext/MakieExt/MakieExt.jl index cdea68ee..652cc885 100644 --- a/ext/MakieExt/MakieExt.jl +++ b/ext/MakieExt/MakieExt.jl @@ -3,9 +3,11 @@ module MakieExt using Statistics using DataFrames using Dates +using LaTeXStrings using Streamfall -import Streamfall: TemporalCrossSection +import Streamfall: StreamfallNetwork +import Streamfall.Analysis: TemporalCrossSection using Makie @@ -24,6 +26,15 @@ function Streamfall.Viz.quickplot(node::NetworkNode, climate::Climate) return f end +function Streamfall.Viz.quickplot(node::DamNode, climate::Climate) + date = timesteps(climate) + + @assert length(date) == length(node.level) || "Date length and result lengths do not match!" + + f, ax, sp = lines(date, node.level) + + return f +end function Streamfall.Viz.quickplot( obs::Vector, node::NetworkNode, climate::Climate; @@ -31,6 +42,12 @@ function Streamfall.Viz.quickplot( ) return quickplot(obs, node.outflow, climate; label, log, burn_in=burn_in, limit=limit, metric=metric) end +function Streamfall.Viz.quickplot( + obs::Vector, node::DamNode, climate::Climate; + label="Modeled", log=false, burn_in=1, limit=nothing, metric=Streamfall.mKGE +) + return quickplot(obs, node.level, climate; label, log, burn_in=burn_in, limit=limit, metric=metric) +end function Streamfall.Viz.quickplot( obs::DataFrame, node::NetworkNode, climate::Climate; @@ -38,6 +55,12 @@ function Streamfall.Viz.quickplot( ) return Streamfall.Viz.quickplot(obs[:, node.name], node.outflow, climate; label, log, burn_in, limit, metric) end +function Streamfall.Viz.quickplot( + obs::DataFrame, node::DamNode, climate::Climate; + label="", log=false, burn_in=1, limit=nothing, metric=Streamfall.mKGE +) + return Streamfall.Viz.quickplot(obs[:, node.name], node.level, climate; label, log, burn_in, limit, metric) +end function Streamfall.Viz.quickplot( obs::Vector, sim::Vector, climate::Climate; @@ -83,9 +106,6 @@ function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, xticklabels::Vector, return f end - -using LaTeXStrings - """ temporal_cross_section( dates, obs; @@ -377,4 +397,17 @@ function Streamfall.Viz.temporal_cross_section( return fig end +""" + save_figure(f::Figure, fn::String) + +Save a figure. +""" +function Streamfall.Viz.save_figure(f::Figure, fn::String) + save(fn, f, update=false) +end +function Streamfall.Viz.save_figure!(fn::String) + f = current_figure() + Streamfall.Viz.save_figure(f, fn) +end + end \ No newline at end of file diff --git a/src/Streamfall.jl b/src/Streamfall.jl index e471bd1d..d717269b 100644 --- a/src/Streamfall.jl +++ b/src/Streamfall.jl @@ -274,7 +274,7 @@ function run_node!( end include("Analysis/Analysis.jl") -include("plotting.jl") +# include("plotting.jl") include("viz/Viz.jl") include("viz/NetworkViz.jl") @@ -299,7 +299,10 @@ export best_candidate, best_fitness, best_params export extract_flow, extract_climate, align_time_frame # Plotting methods -export quickplot, plot_network, save_figure, temporal_cross_section +import .Viz: quickplot, temporal_cross_section, save_figure, save_figure! +import .NetworkViz: plot_network + +export plot_network, quickplot,temporal_cross_section, save_figure, save_figure! # Data interface (climate) export timesteps diff --git a/src/viz/NetworkViz.jl b/src/viz/NetworkViz.jl index e478fd5f..12fbe7c8 100644 --- a/src/viz/NetworkViz.jl +++ b/src/viz/NetworkViz.jl @@ -1,10 +1,7 @@ module NetworkViz -using GraphMakie -using Streamfall +function plot_network end -function plot end - -export plot +export plot_network end diff --git a/src/viz/Viz.jl b/src/viz/Viz.jl index 78175195..76aa0c7f 100644 --- a/src/viz/Viz.jl +++ b/src/viz/Viz.jl @@ -4,6 +4,9 @@ function quickplot end function temporal_cross_section end -export quickplot, temporal_cross_section +function save_figure end +function save_figure! end + +export quickplot, temporal_cross_section, save_figure, save_figure! end \ No newline at end of file From c564e16d17c2bc30ba5fff9ef6cb0174259c9670 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Sun, 5 Jan 2025 12:58:39 +1100 Subject: [PATCH 05/13] Code clean up --- src/Analysis/Analysis.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/Analysis/Analysis.jl b/src/Analysis/Analysis.jl index 1666dc14..ca6401aa 100644 --- a/src/Analysis/Analysis.jl +++ b/src/Analysis/Analysis.jl @@ -1,5 +1,7 @@ module Analysis + include("uncertainty.jl") + end @@ -12,19 +14,19 @@ Gini's Mean Difference. The absolute mean of all pairwise distances between elements in a given set. # References -1. La Haye, R., & Zizler, P. (2019). - The Gini mean difference and variance. - METRON, 77(1), 43-52. +1. La Haye, R., & Zizler, P. (2019). + The Gini mean difference and variance. + METRON, 77(1), 43-52. https://doi.org/10.1007/s40300-019-00149-2 -2. Yitzhaki, S. (2003). +2. Yitzhaki, S. (2003). Gini's Mean difference: A superior measure of variability for non-normal - distributions. + distributions. Metron - International Journal of Statistics, LXI(2), 285-316. https://ideas.repec.org/a/mtn/ancoec/030208.html 3. Kashif, M., Aslam, M., Al-Marshadi, A. H., & Jun, C.-H. (2016). - Capability Indices for Non-Normal Distribution Using Gini's Mean Difference as Measure of Variability. + Capability Indices for Non-Normal Distribution Using Gini's Mean Difference as Measure of Variability. IEEE Access, 4, 7322-7330. https://doi.org/10.1109/ACCESS.2016.2620241 """ From 5ab21aed422e350e76dfd6da2d4def0c83654c3b Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Sat, 25 Jan 2025 23:08:06 +1100 Subject: [PATCH 06/13] Remove Plots and GraphPlots dependencies --- Project.toml | 12 ++---------- src/Network.jl | 2 +- 2 files changed, 3 insertions(+), 11 deletions(-) diff --git a/Project.toml b/Project.toml index 4ce6a4bd..f44957c2 100644 --- a/Project.toml +++ b/Project.toml @@ -17,9 +17,6 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" -GR_jll = "d2c73de3-f751-5644-a686-071e5b155ba9" -GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" @@ -29,7 +26,6 @@ ModelParameters = "4744a3fa-6c31-4707-899e-a3298e4618ad" OnlineStats = "a15396b6-48d5-5d58-9928-6d29437db91e" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" @@ -43,12 +39,12 @@ YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [weakdeps] -Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" [extensions] -MakieExt = "Makie" GraphMakieExt = ["Makie", "GraphMakie"] +MakieExt = "Makie" [compat] BlackBoxOptim = "0.5, 0.6" @@ -62,9 +58,6 @@ DataFrames = "1.6" DataStructures = "0.18.10" Dates = "1.11.0" Distributions = "0.25.100" -GR = "0.72.9, 0.73" -GR_jll = "0.73" -GraphPlot = "0.5, 0.6" Graphs = "1.9" JSON = "0.20.1, 0.21" LaTeXStrings = "1.3" @@ -74,7 +67,6 @@ ModelParameters = "0.4.4" OnlineStats = "1.6" OrderedCollections = "1.6" Parameters = "0.12" -Plots = "1.39" PrettyTables = "2.3" Setfield = "1" Statistics = "1.10, 1.11.1" diff --git a/src/Network.jl b/src/Network.jl index 2d47f982..75b8480c 100644 --- a/src/Network.jl +++ b/src/Network.jl @@ -1,7 +1,7 @@ using OrderedCollections using Cairo, Compose -using Graphs, MetaGraphs, GraphPlot +using Graphs, MetaGraphs using ModelParameters import YAML: load_file, write_file From 50d490e85d4ddc0efea263af369e8014880e5ce5 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Sat, 22 Feb 2025 19:13:56 +1100 Subject: [PATCH 07/13] Formatting changes --- ext/MakieExt/MakieExt.jl | 14 ++++++------- src/plotting.jl | 44 +++++++++++++++++++++------------------- 2 files changed, 30 insertions(+), 28 deletions(-) diff --git a/ext/MakieExt/MakieExt.jl b/ext/MakieExt/MakieExt.jl index 652cc885..bad9bcb4 100644 --- a/ext/MakieExt/MakieExt.jl +++ b/ext/MakieExt/MakieExt.jl @@ -82,15 +82,15 @@ function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, xticklabels::Vector, # Add small constant in case of 0-flow obs = copy(obs) sim = copy(sim) - obs[obs .== 0.0] .+= 1e-4 - sim[sim .== 0.0] .+= 1e-4 + obs[obs.==0.0] .+= 1e-4 + sim[sim.==0.0] .+= 1e-4 end scale = log == false ? identity : log10 - f = Figure(size=(850,400)) - flow_ax = Axis(f[1,1]; xlabel="Date", ylabel="Streamflow", yscale=scale) - qq_ax = Axis(f[1,2]; xlabel="Observed", ylabel="Modeled", xscale=scale, yscale=scale) + f = Figure(size=(850, 400)) + flow_ax = Axis(f[1, 1]; xlabel="Date", ylabel="Streamflow", yscale=scale) + qq_ax = Axis(f[1, 2]; xlabel="Observed", ylabel="Modeled", xscale=scale, yscale=scale) label = "$(label) ($(metric_name): $(score))" lines!(flow_ax, xticklabels, obs, label="Observed") @@ -218,7 +218,7 @@ function Streamfall.Viz.temporal_cross_section( ax.xticks = (1:12:length(xlabels), xlabels[1:12:end]) end - ax.xticklabelrotation = π/4 + ax.xticklabelrotation = π / 4 # Apply log scale if specified if !isnothing(tmp) && (tmp in logscale) @@ -360,7 +360,7 @@ function Streamfall.Viz.temporal_cross_section( ax.xticks = (1:12:length(xlabels), xlabels[1:12:end]) end - ax.xticklabelrotation = π/4 + ax.xticklabelrotation = π / 4 # Apply log scale if specified if !isnothing(tmp) && (tmp in logscale) diff --git a/src/plotting.jl b/src/plotting.jl index a0c6fde6..29664af5 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -75,7 +75,7 @@ function quickplot(obs::Vector, sim::Vector, xticklabels::Vector, label="Modeled yaxis!(qqfig, :log10) end - combined = plot(fig, qqfig, size=(1000, 500), left_margin=10mm, bottom_margin=5mm, layout=(1,2)) + combined = plot(fig, qqfig, size=(1000, 500), left_margin=10mm, bottom_margin=5mm, layout=(1, 2)) return combined end @@ -96,9 +96,9 @@ Plot residual between two sequences. function plot_residuals(x::Array, y::Array; xlabel="", ylabel="", title="") # 1:1 Plot fig_1to1 = scatter(x, y, legend=false, - markerstrokewidth=0, markerstrokealpha=0, alpha=0.2) + markerstrokewidth=0, markerstrokealpha=0, alpha=0.2) plot!(x, y, color=:red, markersize=0.1, markerstrokewidth=0, - xlabel=xlabel, ylabel=ylabel, title=title) + xlabel=xlabel, ylabel=ylabel, title=title) return fig_1to1 end @@ -130,9 +130,9 @@ Filters out leap days. - `period::Function` : Method from `Dates` package to group (defaults to `monthday`) """ function temporal_cross_section(dates, obs; - title="", ylabel="ME", label=nothing, - period::Function=monthday, - kwargs...) # show_extremes::Bool=false, + title="", ylabel="ME", label=nothing, + period::Function=monthday, + kwargs...) # show_extremes::Bool=false, if isnothing(label) label = ylabel end @@ -190,20 +190,22 @@ function temporal_cross_section(dates, obs; fig = plot(xlabels, lower_95, fillrange=upper_95, color="lightblue", alpha=0.3, label="CI₉₅ μ: $(wr95_m_ind), σ: $(wr95_sd_ind)", linealpha=0) plot!(fig, xlabels, lower_75, fillrange=upper_75, color="lightblue", alpha=0.5, label="CI₇₅ μ: $(wr75_m_ind), σ: $(wr75_sd_ind)", linealpha=0) - plot!(fig, xlabels, x_section, - label="Mean of $(label) μ: $(m_ind), σ: $(sd_ind)", - color="black", - xlabel=nameof(period), - ylabel=ylabel, - legend=:bottomleft, - legendfont=Plots.font(10), - fg_legend=:transparent, - bg_legend=:transparent, - left_margin=5mm, - bottom_margin=5mm, - title=title, - yformatter=format_func; - kwargs...) + plot!( + fig, xlabels, x_section, + label="Mean of $(label) μ: $(m_ind), σ: $(sd_ind)", + color="black", + xlabel=nameof(period), + ylabel=ylabel, + legend=:bottomleft, + legendfont=Plots.font(10), + fg_legend=:transparent, + bg_legend=:transparent, + left_margin=5mm, + bottom_margin=5mm, + title=title, + yformatter=format_func; + kwargs... + ) # if show_extremes # scatter!(fig, xlabels, min_section, label="", alpha=0.5, color="lightblue", markerstrokewidth=0; kwargs...) @@ -346,7 +348,7 @@ function plot(node::DamNode, climate::Climate) levels = plot(timesteps(climate), node.level, xlabel="Date", ylabel="Dam Level", label=node.name) outflows = plot(timesteps(climate), node.outflow, xlabel="Date", ylabel="Discharge", label=node.name) - combined = plot(levels, outflows, size=(1000, 400), left_margin=10mm, bottom_margin=5mm, layout=(1,2)) + combined = plot(levels, outflows, size=(1000, 400), left_margin=10mm, bottom_margin=5mm, layout=(1, 2)) return combined end \ No newline at end of file From c37a964a7476eb64f62eaec618abba5836409e81 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Fri, 25 Apr 2025 17:22:24 +1000 Subject: [PATCH 08/13] Separate out all plotting methods and add support for Plots.jl --- Project.toml | 5 +- ext/GraphPlotExt/GraphPlotExt.jl | 34 +++++ ext/MakieExt/MakieExt.jl | 43 +++++- src/plotting.jl => ext/PlotsExt/PlotsExt.jl | 141 ++++++++++---------- src/Network.jl | 29 ---- src/Streamfall.jl | 3 +- src/viz/NetworkViz.jl | 3 +- src/viz/Viz.jl | 10 +- 8 files changed, 160 insertions(+), 108 deletions(-) create mode 100644 ext/GraphPlotExt/GraphPlotExt.jl rename src/plotting.jl => ext/PlotsExt/PlotsExt.jl (70%) diff --git a/Project.toml b/Project.toml index f44957c2..d04c4054 100644 --- a/Project.toml +++ b/Project.toml @@ -41,10 +41,14 @@ ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [weakdeps] GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231" [extensions] GraphMakieExt = ["Makie", "GraphMakie"] MakieExt = "Makie" +GraphPlotExt = ["Plots", "GraphPlot"] +PlotsExt = "Plots" [compat] BlackBoxOptim = "0.5, 0.6" @@ -72,7 +76,6 @@ Setfield = "1" Statistics = "1.10, 1.11.1" StatsBase = "0.34" StatsFuns = "1.3" -StatsPlots = "0.15" YAML = "0.4" ZipFile = "0.10" julia = "1" diff --git a/ext/GraphPlotExt/GraphPlotExt.jl b/ext/GraphPlotExt/GraphPlotExt.jl new file mode 100644 index 00000000..2ce4eda4 --- /dev/null +++ b/ext/GraphPlotExt/GraphPlotExt.jl @@ -0,0 +1,34 @@ +module GraphPlotExt + +using Plots, GraphPlot +using Streamfall +import Streamfall: vertices +import Streamfall: StreamfallNetwork + +""" + plot_network(sn::StreamfallNetwork) + +Simple plot of stream network. +""" +function Streamfall.NetworkViz.plot_network(sn::StreamfallNetwork; as_html::Bool=false) + node_labels = ["$(sn[i].name)\n" * string(nameof(typeof(sn[i]))) for i in vertices(sn.mg)] + + if as_html + plot_func = gplothtml + else + plot_func = gplot + end + + plot_func(sn.mg, nodelabel=node_labels) +end + +""" + save_figure(sn::StreamfallNetwork, fn::String) + +Save a figure of the network in SVG format. +""" +function Streamfall.NetworkViz.save_figure(sn::StreamfallNetwork, fn::String) + draw(SVG(fn, 16cm, 16cm), plot_network(sn)) +end + +end diff --git a/ext/MakieExt/MakieExt.jl b/ext/MakieExt/MakieExt.jl index bad9bcb4..e36bfd6f 100644 --- a/ext/MakieExt/MakieExt.jl +++ b/ext/MakieExt/MakieExt.jl @@ -97,7 +97,7 @@ function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, xticklabels::Vector, lines!(flow_ax, xticklabels, sim; label=label, alpha=0.5) qqplot!( - qq_ax, obs, sim; + qq_ax, obs[burn_in:end], sim[burn_in:end]; qqline=:identity, strokewidth=0.01, strokecolor=(:blue, 0.01), # legend=false, strokewidth=0.03, strokealpha=0.1, markercolor=(:blue, 0.02) ) @@ -106,6 +106,47 @@ function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, xticklabels::Vector, return f end +""" + plot_residuals(obs::AbstractVector, sim::AbstractVector; xlabel="", ylabel="", title="") + +Plot residual between two sequences. + +# Arguments +- x : x-axis data +- y : y-axis data +- xlabel : x-axis label +- ylabel : y-axis label +- title : title text +""" +function Streamfall.Viz.plot_residuals(x::AbstractVector, y::AbstractVector; xlabel="", ylabel="", title="") + # 1:1 Plot + f, ax, sp = scatter(x, y, strokewidth=0, alpha=0.2) + + # Calculate means + x_mean = mean(x) + y_mean = mean(y) + + # Calculate slope (β) + β = sum((x .- x_mean) .* (y .- y_mean)) / sum((x .- x_mean) .^ 2) + + # Calculate intercept (α) + α = y_mean - β * x_mean + + # Approximate x values for the line + x_line = range(minimum(x), maximum(x), length=length(x)) + + # Calculate corresponding y values + y_line = α .+ β .* x_line + + ax.xlabel = xlabel + ax.ylabel = ylabel + ax.title = title + + lines!(ax, y_line; color=(:black, 0.3)) + + return f +end + """ temporal_cross_section( dates, obs; diff --git a/src/plotting.jl b/ext/PlotsExt/PlotsExt.jl similarity index 70% rename from src/plotting.jl rename to ext/PlotsExt/PlotsExt.jl index 29664af5..5d1d8269 100644 --- a/src/plotting.jl +++ b/ext/PlotsExt/PlotsExt.jl @@ -1,40 +1,72 @@ +module PlotsExt + using Plots, StatsPlots using Plots.Measures -import Plots: plot, plot! + using DataFrames, Dates, Statistics, Distributions, LaTeXStrings import Bootstrap: bootstrap, BalancedSampling -import .Analysis: TemporalCrossSection +using Streamfall +import Streamfall: Analysis.TemporalCrossSection +function Streamfall.Viz.plot(node::NetworkNode, climate::Climate) + return Plots.plot( + timesteps(climate), + node.outflow, + label=node.name, + xlabel="Date", + ylabel="Outflows", + ) +end -function quickplot(node::NetworkNode) - fig = plot(node.outflow) +function Streamfall.Viz.plot!(node::NetworkNode, climate::Climate) + Plots.plot!( + timesteps(climate), + node.outflow, + label=node.name + ) +end +function Streamfall.Viz.plot!(f, node::NetworkNode, climate::Climate) + Plots.plot!(f, timesteps(climate), node.outflow, label=node.name) +end + +function Streamfall.Viz.plot(node::DamNode, climate::Climate) + levels = Plots.plot(timesteps(climate), node.level, xlabel="Date", ylabel="Dam Level", label=node.name) + outflows = Plots.plot(timesteps(climate), node.outflow, xlabel="Date", ylabel="Discharge", label=node.name) + + combined = Plots.plot(levels, outflows, size=(1000, 400), left_margin=10mm, bottom_margin=5mm, layout=(1, 2)) + + return combined +end + +function Streamfall.Viz.quickplot(node::NetworkNode) + fig = Plots.plot(node.outflow) return fig end -function quickplot(node::NetworkNode, climate::Climate) - date = timesteps(climate) +function Streamfall.Viz.quickplot(node::NetworkNode, climate::Climate) + all_dates = timesteps(climate) - @assert length(date) == length(node.outflow) || "Date length and result lengths do not match!" + @assert length(all_dates) == length(node.outflow) || "Date length and result lengths do not match!" - fig = plot(date, node.outflow) + fig = plot(all_dates, node.outflow) return fig end -function quickplot(obs, node::NetworkNode, climate::Climate, label="", log=false; burn_in=1, limit=nothing, metric=Streamfall.mKGE) - return quickplot(obs, node.outflow, climate, label, log; burn_in=burn_in, limit=limit, metric=metric) +function Streamfall.Viz.quickplot(obs, node::NetworkNode, climate::Climate, label="", log=false; burn_in=1, limit=nothing, metric=Streamfall.mKGE) + return Streamfall.Viz.quickplot(obs, node.outflow, climate, label, log; burn_in=burn_in, limit=limit, metric=metric) end -function quickplot(obs::DataFrame, sim::Vector, climate::Climate, label="", log=false; burn_in=1, limit=nothing, metric=Streamfall.mKGE) - return quickplot(Matrix(obs[:, Not("Date")])[:, 1], sim, climate, label, log; burn_in, limit, metric) +function Streamfall.Viz.quickplot(obs::DataFrame, sim::Vector, climate::Climate, label="", log=false; burn_in=1, limit=nothing, metric=Streamfall.mKGE) + return Streamfall.Viz.quickplot(Matrix(obs[:, Not("Date")])[:, 1], sim, climate, label, log; burn_in, limit, metric) end -function quickplot(obs::Vector, sim::Vector, climate::Climate, label="", log=false; burn_in=1, limit=nothing, metric=Streamfall.mKGE) +function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, climate::Climate, label="", log=false; burn_in=1, limit=nothing, metric=Streamfall.mKGE) date = timesteps(climate) last_e = !isnothing(limit) ? limit : lastindex(obs) show_range = burn_in:last_e return quickplot(obs[show_range], sim[show_range], date[show_range], label, log; metric=metric) end -function quickplot(obs::Vector, sim::Vector, xticklabels::Vector, label="Modeled", log=false; metric=Streamfall.mKGE) +function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, xticklabels::Vector, label="Modeled", log=false; metric=Streamfall.mKGE) @assert length(xticklabels) == length(obs) || "x-axis tick label length and observed lengths do not match!" @assert length(xticklabels) == length(sim) || "x-axis tick label length and simulated lengths do not match!" @@ -48,8 +80,8 @@ function quickplot(obs::Vector, sim::Vector, xticklabels::Vector, label="Modeled end label = "$(label) ($(metric_name): $(score))" - fig = plot( - xticklabels, obs, + fig = Plots.plot( + xticklabels, obs; label="Observed", legend=:best, ylabel="Streamflow", @@ -75,14 +107,14 @@ function quickplot(obs::Vector, sim::Vector, xticklabels::Vector, label="Modeled yaxis!(qqfig, :log10) end - combined = plot(fig, qqfig, size=(1000, 500), left_margin=10mm, bottom_margin=5mm, layout=(1, 2)) + combined = Plots.plot(fig, qqfig, size=(1000, 500), left_margin=10mm, bottom_margin=5mm, layout=(1, 2)) return combined end """ - plot_residuals(obs::Array, sim::Array; xlabel="", ylabel="", title="") + plot_residuals(obs::AbstractVector, sim::AbstractVector; xlabel="", ylabel="", title="") Plot residual between two sequences. @@ -93,27 +125,19 @@ Plot residual between two sequences. - ylabel : y-axis label - title : title text """ -function plot_residuals(x::Array, y::Array; xlabel="", ylabel="", title="") +function Streamfall.Viz.plot_residuals(x::AbstractVector, y::AbstractVector; xlabel="", ylabel="", title="") # 1:1 Plot - fig_1to1 = scatter(x, y, legend=false, - markerstrokewidth=0, markerstrokealpha=0, alpha=0.2) - plot!(x, y, color=:red, markersize=0.1, markerstrokewidth=0, - xlabel=xlabel, ylabel=ylabel, title=title) + fig_1to1 = scatter( + x, y; + legend=false, + markerstrokewidth=0, markerstrokealpha=0, alpha=0.2, smooth=true, + xlabel=xlabel, ylabel=ylabel, title=title + ) return fig_1to1 end -"""Symmetrical log values. - -https://kar.kent.ac.uk/32810/2/2012_Bi-symmetric-log-transformation_v5.pdf -https://discourse.julialang.org/t/symmetrical-log-plot/45709/3 -""" -function symlog(y) - return sign.(y) .* log10.(1.0 .+ abs.(y)) -end - - """ temporal_cross_section(dates, obs; ylabel=nothing, period::Function=monthday) @@ -129,7 +153,8 @@ Filters out leap days. - ylabel : Optional replacement ylabel. Uses name of `func` if not provided. - `period::Function` : Method from `Dates` package to group (defaults to `monthday`) """ -function temporal_cross_section(dates, obs; +function Streamfall.Viz.temporal_cross_section( + dates, obs; title="", ylabel="ME", label=nothing, period::Function=monthday, kwargs...) # show_extremes::Bool=false, @@ -188,10 +213,9 @@ function temporal_cross_section(dates, obs; wr95_m_ind = round(xsect_res.mean_95, digits=2) wr95_sd_ind = round(xsect_res.std_95, digits=2) - fig = plot(xlabels, lower_95, fillrange=upper_95, color="lightblue", alpha=0.3, label="CI₉₅ μ: $(wr95_m_ind), σ: $(wr95_sd_ind)", linealpha=0) - plot!(fig, xlabels, lower_75, fillrange=upper_75, color="lightblue", alpha=0.5, label="CI₇₅ μ: $(wr75_m_ind), σ: $(wr75_sd_ind)", linealpha=0) - plot!( - fig, xlabels, x_section, + fig = Plots.plot(xlabels, lower_95; fillrange=upper_95, color="lightblue", alpha=0.3, label="CI₉₅ μ: $(wr95_m_ind), σ: $(wr95_sd_ind)", linealpha=0) + Plots.plot!(fig, xlabels, lower_75; fillrange=upper_75, color="lightblue", alpha=0.5, label="CI₇₅ μ: $(wr75_m_ind), σ: $(wr75_sd_ind)", linealpha=0) + Plots.plot!(fig, xlabels, x_section; label="Mean of $(label) μ: $(m_ind), σ: $(sd_ind)", color="black", xlabel=nameof(period), @@ -203,9 +227,8 @@ function temporal_cross_section(dates, obs; left_margin=5mm, bottom_margin=5mm, title=title, - yformatter=format_func; - kwargs... - ) + yformatter=format_func, + kwargs...) # if show_extremes # scatter!(fig, xlabels, min_section, label="", alpha=0.5, color="lightblue", markerstrokewidth=0; kwargs...) @@ -237,7 +260,7 @@ Filters out leap days. - `label` : Optional legend label. Uses `ylabel` if not provided. - `period` : Method from `Dates` package to group (defaults to `month`) """ -function temporal_cross_section( +function Streamfall.Viz.temporal_cross_section( dates, obs, sim; title="", ylabel="Median Error", label=nothing, period::Function=monthday, kwargs... ) # show_extremes::Bool=false, @@ -296,9 +319,9 @@ function temporal_cross_section( wr95_m_ind = round(xsect_res.mean_95, digits=2) wr95_sd_ind = round(xsect_res.std_95, digits=2) - fig = plot(xlabels, lower_95, fillrange=upper_95, color="lightblue", alpha=0.3, label="CI₉₅ μ: $(wr95_m_ind), σ: $(wr95_sd_ind)", linealpha=0) - plot!(fig, xlabels, lower_75, fillrange=upper_75, color="lightblue", alpha=0.5, label="CI₇₅ μ: $(wr75_m_ind), σ: $(wr75_sd_ind)", linealpha=0) - plot!( + fig = Plots.plot(xlabels, lower_95, fillrange=upper_95, color="lightblue", alpha=0.3, label="CI₉₅ μ: $(wr95_m_ind), σ: $(wr95_sd_ind)", linealpha=0) + Plots.plot!(fig, xlabels, lower_75, fillrange=upper_75, color="lightblue", alpha=0.5, label="CI₇₅ μ: $(wr75_m_ind), σ: $(wr75_sd_ind)", linealpha=0) + Plots.plot!( fig, xlabels, x_section, label="$(label) μ: $(m_ind), σ: $(sd_ind)", color="black", @@ -323,32 +346,4 @@ function temporal_cross_section( return fig end -function plot(node::NetworkNode, climate::Climate) - return plot( - timesteps(climate), - node.outflow, - label=node.name, - xlabel="Date", - ylabel="Outflows", - ) -end - -function plot!(node::NetworkNode, climate::Climate) - plot!( - timesteps(climate), - node.outflow, - label=node.name - ) -end -function plot!(f, node::NetworkNode, climate::Climate) - plot!(f, timesteps(climate), node.outflow, label=node.name) end - -function plot(node::DamNode, climate::Climate) - levels = plot(timesteps(climate), node.level, xlabel="Date", ylabel="Dam Level", label=node.name) - outflows = plot(timesteps(climate), node.outflow, xlabel="Date", ylabel="Discharge", label=node.name) - - combined = plot(levels, outflows, size=(1000, 400), left_margin=10mm, bottom_margin=5mm, layout=(1, 2)) - - return combined -end \ No newline at end of file diff --git a/src/Network.jl b/src/Network.jl index 75b8480c..2965a5fd 100644 --- a/src/Network.jl +++ b/src/Network.jl @@ -334,35 +334,6 @@ function Base.show(io::IO, sn::StreamfallNetwork) end end - -""" - plot_network(sn::StreamfallNetwork) - -Simple plot of stream network. -""" -function plot_network(sn::StreamfallNetwork; as_html=false) - node_labels = ["$(sn[i].name)\n" * string(nameof(typeof(sn[i]))) for i in vertices(sn.mg)] - - if as_html - plot_func = gplothtml - else - plot_func = gplot - end - - plot_func(sn.mg, nodelabel=node_labels) -end - - -""" - save_figure(sn::StreamfallNetwork, fn::String) - -Save a figure of the network in SVG format. -""" -function save_figure(sn::StreamfallNetwork, fn::String) - draw(SVG(fn, 16cm, 16cm), plot_network(sn)) -end - - function Base.iterate(sn::StreamfallNetwork, state=1) if state > length(sn) return nothing diff --git a/src/Streamfall.jl b/src/Streamfall.jl index d717269b..9b46901c 100644 --- a/src/Streamfall.jl +++ b/src/Streamfall.jl @@ -274,7 +274,6 @@ function run_node!( end include("Analysis/Analysis.jl") -# include("plotting.jl") include("viz/Viz.jl") include("viz/NetworkViz.jl") @@ -302,7 +301,7 @@ export extract_flow, extract_climate, align_time_frame import .Viz: quickplot, temporal_cross_section, save_figure, save_figure! import .NetworkViz: plot_network -export plot_network, quickplot,temporal_cross_section, save_figure, save_figure! +export plot_network, quickplot, temporal_cross_section, save_figure, save_figure! # Data interface (climate) export timesteps diff --git a/src/viz/NetworkViz.jl b/src/viz/NetworkViz.jl index 12fbe7c8..3783678d 100644 --- a/src/viz/NetworkViz.jl +++ b/src/viz/NetworkViz.jl @@ -1,7 +1,8 @@ module NetworkViz function plot_network end +function save_figure end -export plot_network +export plot_network, save_figure end diff --git a/src/viz/Viz.jl b/src/viz/Viz.jl index 76aa0c7f..b9075715 100644 --- a/src/viz/Viz.jl +++ b/src/viz/Viz.jl @@ -1,12 +1,20 @@ module Viz + +function plot end +function plot! end function quickplot end function temporal_cross_section end +function plot_residuals end + function save_figure end function save_figure! end -export quickplot, temporal_cross_section, save_figure, save_figure! +export plot, plot!, quickplot +export temporal_cross_section, plot_residuals +export save_figure, save_figure! + end \ No newline at end of file From 1feaf9b80eb66d4284aafc71f7c6c17e4d7f84eb Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Fri, 25 Apr 2025 17:22:51 +1000 Subject: [PATCH 09/13] Committing deprecated file Marked as new file so committing so it is in history. --- src/viz/_plotting.jl | 354 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 354 insertions(+) create mode 100644 src/viz/_plotting.jl diff --git a/src/viz/_plotting.jl b/src/viz/_plotting.jl new file mode 100644 index 00000000..29664af5 --- /dev/null +++ b/src/viz/_plotting.jl @@ -0,0 +1,354 @@ +using Plots, StatsPlots +using Plots.Measures +import Plots: plot, plot! +using DataFrames, Dates, Statistics, Distributions, LaTeXStrings +import Bootstrap: bootstrap, BalancedSampling + +import .Analysis: TemporalCrossSection + + +function quickplot(node::NetworkNode) + fig = plot(node.outflow) + return fig +end + + +function quickplot(node::NetworkNode, climate::Climate) + date = timesteps(climate) + + @assert length(date) == length(node.outflow) || "Date length and result lengths do not match!" + + fig = plot(date, node.outflow) + + return fig +end +function quickplot(obs, node::NetworkNode, climate::Climate, label="", log=false; burn_in=1, limit=nothing, metric=Streamfall.mKGE) + return quickplot(obs, node.outflow, climate, label, log; burn_in=burn_in, limit=limit, metric=metric) +end +function quickplot(obs::DataFrame, sim::Vector, climate::Climate, label="", log=false; burn_in=1, limit=nothing, metric=Streamfall.mKGE) + return quickplot(Matrix(obs[:, Not("Date")])[:, 1], sim, climate, label, log; burn_in, limit, metric) +end +function quickplot(obs::Vector, sim::Vector, climate::Climate, label="", log=false; burn_in=1, limit=nothing, metric=Streamfall.mKGE) + date = timesteps(climate) + last_e = !isnothing(limit) ? limit : lastindex(obs) + show_range = burn_in:last_e + return quickplot(obs[show_range], sim[show_range], date[show_range], label, log; metric=metric) +end +function quickplot(obs::Vector, sim::Vector, xticklabels::Vector, label="Modeled", log=false; metric=Streamfall.mKGE) + @assert length(xticklabels) == length(obs) || "x-axis tick label length and observed lengths do not match!" + @assert length(xticklabels) == length(sim) || "x-axis tick label length and simulated lengths do not match!" + + score = round(metric(obs, sim), digits=4) + metric_name = String(Symbol(metric)) + + if log + # Add small constant in case of 0-flow + obs = obs .+ 1e-2 + sim = sim .+ 1e-2 + end + + label = "$(label) ($(metric_name): $(score))" + fig = plot( + xticklabels, obs, + label="Observed", + legend=:best, + ylabel="Streamflow", + xlabel="Date", + fg_legend=:transparent, + bg_legend=:transparent + ) + plot!(xticklabels, sim, label=label, alpha=0.5) + + if log + # modify yaxis + yaxis!(fig, :log10) + end + + qqfig = qqplot( + obs, sim, + legend=false, markerstrokewidth=0.03, markerstrokealpha=0.1, markeralpha=0.2, + xlabel="Observed", ylabel="Modeled" + ) + + if log + xaxis!(qqfig, :log10) + yaxis!(qqfig, :log10) + end + + combined = plot(fig, qqfig, size=(1000, 500), left_margin=10mm, bottom_margin=5mm, layout=(1, 2)) + + return combined +end + + +""" + plot_residuals(obs::Array, sim::Array; xlabel="", ylabel="", title="") + +Plot residual between two sequences. + +# Arguments +- x : x-axis data +- y : y-axis data +- xlabel : x-axis label +- ylabel : y-axis label +- title : title text +""" +function plot_residuals(x::Array, y::Array; xlabel="", ylabel="", title="") + # 1:1 Plot + fig_1to1 = scatter(x, y, legend=false, + markerstrokewidth=0, markerstrokealpha=0, alpha=0.2) + plot!(x, y, color=:red, markersize=0.1, markerstrokewidth=0, + xlabel=xlabel, ylabel=ylabel, title=title) + + return fig_1to1 +end + + +"""Symmetrical log values. + +https://kar.kent.ac.uk/32810/2/2012_Bi-symmetric-log-transformation_v5.pdf +https://discourse.julialang.org/t/symmetrical-log-plot/45709/3 +""" +function symlog(y) + return sign.(y) .* log10.(1.0 .+ abs.(y)) +end + + +""" + temporal_cross_section(dates, obs; ylabel=nothing, period::Function=monthday) + +Provides indication of temporal variation and uncertainty across time, grouped by `period`. + +Notes: +Assumes daily data. +Filters out leap days. + +# Arguments +- dates : Date of each observation +- obs : observed data +- ylabel : Optional replacement ylabel. Uses name of `func` if not provided. +- `period::Function` : Method from `Dates` package to group (defaults to `monthday`) +""" +function temporal_cross_section(dates, obs; + title="", ylabel="ME", label=nothing, + period::Function=monthday, + kwargs...) # show_extremes::Bool=false, + if isnothing(label) + label = ylabel + end + + arg_keys = keys(kwargs) + format_func = y -> y + logscale = [:log, :log10] + tmp = nothing + + xsect_res = TemporalCrossSection(dates, obs, period) + + if :yscale in arg_keys || :yaxis in arg_keys + tmp = (:yscale in arg_keys) ? kwargs[:yscale] : kwargs[:yaxis] + + if tmp in logscale + log_obs = symlog(copy(obs)) + + # Format function for y-axis tick labels (e.g., 10^x) + format_func = y -> (y != 0) ? L"%$(Int(round(sign(y)) * 10))^{%$(round(abs(y), digits=1))}" : L"0" + + log_xsect_res = TemporalCrossSection(dates, log_obs, period) + target = log_xsect_res.cross_section + else + target = xsect_res.cross_section + end + else + target = xsect_res.cross_section + end + + x_section = target.median + lower_75 = target.lower_75 + upper_75 = target.upper_75 + lower_95 = target.lower_95 + upper_95 = target.upper_95 + + sp = target.subperiod + xlabels = join.(sp, "-") + + if !isnothing(tmp) & (tmp in logscale) + # Remove keys so later plotting methods do not error + kwargs = Dict(kwargs) + delete!(kwargs, :yscale) + delete!(kwargs, :yaxis) + end + + # Display indicator values using original data instead of log-transformed data + med = xsect_res.cross_section.median + m_ind = round(mean(med), digits=2) + sd_ind = round(std(med), digits=2) + + wr75_m_ind = round(xsect_res.mean_75, digits=2) + wr75_sd_ind = round(xsect_res.std_75, digits=2) + wr95_m_ind = round(xsect_res.mean_95, digits=2) + wr95_sd_ind = round(xsect_res.std_95, digits=2) + + fig = plot(xlabels, lower_95, fillrange=upper_95, color="lightblue", alpha=0.3, label="CI₉₅ μ: $(wr95_m_ind), σ: $(wr95_sd_ind)", linealpha=0) + plot!(fig, xlabels, lower_75, fillrange=upper_75, color="lightblue", alpha=0.5, label="CI₇₅ μ: $(wr75_m_ind), σ: $(wr75_sd_ind)", linealpha=0) + plot!( + fig, xlabels, x_section, + label="Mean of $(label) μ: $(m_ind), σ: $(sd_ind)", + color="black", + xlabel=nameof(period), + ylabel=ylabel, + legend=:bottomleft, + legendfont=Plots.font(10), + fg_legend=:transparent, + bg_legend=:transparent, + left_margin=5mm, + bottom_margin=5mm, + title=title, + yformatter=format_func; + kwargs... + ) + + # if show_extremes + # scatter!(fig, xlabels, min_section, label="", alpha=0.5, color="lightblue", markerstrokewidth=0; kwargs...) + # scatter!(fig, xlabels, max_section, label="", alpha=0.5, color="lightblue", markerstrokewidth=0; kwargs...) + # end + + return fig +end + + +""" + temporal_cross_section( + dates, obs, sim; + title="", ylabel="Median Error", label=nothing, period::Function=monthday, kwargs... + ) + +Provides indication of temporal variation and uncertainty across time, grouped by `period`. + +Notes: +Assumes daily data. +Filters out leap days. + +# Arguments +- `dates` : Date of each observation +- `obs` : Observed data +- `sim` : Simulated data +- `title` : Optional plot title. Blank if not provided. +- `ylabel` : Optional replacement ylabel. Uses name of `period` if not provided. +- `label` : Optional legend label. Uses `ylabel` if not provided. +- `period` : Method from `Dates` package to group (defaults to `month`) +""" +function temporal_cross_section( + dates, obs, sim; + title="", ylabel="Median Error", label=nothing, period::Function=monthday, kwargs... +) # show_extremes::Bool=false, + if isnothing(label) + label = ylabel + end + + arg_keys = keys(kwargs) + format_func = y -> y + logscale = [:log, :log10] + tmp = nothing + + xsect_res = TemporalCrossSection(dates, obs, sim, period) + target = xsect_res.cross_section + + if :yscale in arg_keys || :yaxis in arg_keys + tmp = (:yscale in arg_keys) ? kwargs[:yscale] : kwargs[:yaxis] + + if tmp in logscale + log_obs = symlog(copy(obs)) + log_sim = symlog(copy(sim)) + + # Format function for y-axis tick labels (e.g., 10^x) + format_func = y -> (y != 0) ? L"%$(Int(round(sign(y)) * 10))^{%$(round(abs(y), digits=1))}" : L"0" + + log_xsect_res = TemporalCrossSection(dates, log_obs, log_sim, period) + target = log_xsect_res.cross_section + end + end + + x_section = target.median + lower_75 = target.lower_75 + upper_75 = target.upper_75 + lower_95 = target.lower_95 + upper_95 = target.upper_95 + + sp = target.subperiod + xlabels = join.(sp, "-") + + if !isnothing(tmp) & (tmp in logscale) + # Remove keys so later plotting methods do not error + kwargs = Dict(kwargs) + delete!(kwargs, :yscale) + delete!(kwargs, :yaxis) + end + + xsect_df = xsect_res.cross_section + + # Display indicator values using original data instead of log-transformed data + med = xsect_df.median + m_ind = round(mean(med), digits=2) + sd_ind = round(std(med), digits=2) + + wr75_m_ind = round(xsect_res.mean_75, digits=2) + wr75_sd_ind = round(xsect_res.std_75, digits=2) + wr95_m_ind = round(xsect_res.mean_95, digits=2) + wr95_sd_ind = round(xsect_res.std_95, digits=2) + + fig = plot(xlabels, lower_95, fillrange=upper_95, color="lightblue", alpha=0.3, label="CI₉₅ μ: $(wr95_m_ind), σ: $(wr95_sd_ind)", linealpha=0) + plot!(fig, xlabels, lower_75, fillrange=upper_75, color="lightblue", alpha=0.5, label="CI₇₅ μ: $(wr75_m_ind), σ: $(wr75_sd_ind)", linealpha=0) + plot!( + fig, xlabels, x_section, + label="$(label) μ: $(m_ind), σ: $(sd_ind)", + color="black", + xlabel=nameof(period), + ylabel=ylabel, + legend=:bottomleft, + legendfont=Plots.font(10), + fg_legend=:transparent, + bg_legend=:transparent, + left_margin=5mm, + bottom_margin=5mm, + title=title, + yformatter=format_func; + kwargs... + ) + + # if show_extremes + # scatter!(fig, xlabels, min_section, label="", alpha=0.5, color="lightblue", markerstrokewidth=0; kwargs...) + # scatter!(fig, xlabels, max_section, label="", alpha=0.5, color="lightblue", markerstrokewidth=0; kwargs...) + # end + + return fig +end + +function plot(node::NetworkNode, climate::Climate) + return plot( + timesteps(climate), + node.outflow, + label=node.name, + xlabel="Date", + ylabel="Outflows", + ) +end + +function plot!(node::NetworkNode, climate::Climate) + plot!( + timesteps(climate), + node.outflow, + label=node.name + ) +end +function plot!(f, node::NetworkNode, climate::Climate) + plot!(f, timesteps(climate), node.outflow, label=node.name) +end + +function plot(node::DamNode, climate::Climate) + levels = plot(timesteps(climate), node.level, xlabel="Date", ylabel="Dam Level", label=node.name) + outflows = plot(timesteps(climate), node.outflow, xlabel="Date", ylabel="Discharge", label=node.name) + + combined = plot(levels, outflows, size=(1000, 400), left_margin=10mm, bottom_margin=5mm, layout=(1, 2)) + + return combined +end \ No newline at end of file From a560076b5bc5e6f70a11371f3ee4a64583c36c37 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Fri, 25 Apr 2025 17:25:42 +1000 Subject: [PATCH 10/13] Explicit comment detailing import to activate extension --- docs/src/examples/calibration.md | 3 ++- docs/src/examples/network_loading.md | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/src/examples/calibration.md b/docs/src/examples/calibration.md index 11e77d74..9da7ec3b 100644 --- a/docs/src/examples/calibration.md +++ b/docs/src/examples/calibration.md @@ -16,7 +16,8 @@ using Statistics using CSV, DataFrames, YAML using Streamfall -using Plots +# Import visualization packages to compile extensions +using Plots, GraphPlot sn = load_network("Example Network", "../test/data/campaspe/campaspe_network.yml") diff --git a/docs/src/examples/network_loading.md b/docs/src/examples/network_loading.md index 7bbf396c..90930807 100644 --- a/docs/src/examples/network_loading.md +++ b/docs/src/examples/network_loading.md @@ -4,6 +4,7 @@ Loading a pre-defined network from a YAML file. ```julia using YAML +using Plots, GraphPlot using Streamfall From 4841f44d72159b45b9ce59cb296afa4e7ed1bf79 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Fri, 25 Apr 2025 23:09:35 +1000 Subject: [PATCH 11/13] Standardize interfaces --- Project.toml | 9 +++++---- README.md | 17 ++++++++++++----- ext/GraphPlotExt/GraphPlotExt.jl | 2 +- ext/MakieExt/MakieExt.jl | 2 +- ext/PlotsExt/PlotsExt.jl | 27 ++++++++++++++------------- src/Analysis/uncertainty.jl | 20 +++++++++++--------- src/Streamfall.jl | 2 -- src/viz/Viz.jl | 10 ++++++++++ src/viz/_plotting.jl | 14 +++++++------- 9 files changed, 61 insertions(+), 42 deletions(-) diff --git a/Project.toml b/Project.toml index d04c4054..1c39417f 100644 --- a/Project.toml +++ b/Project.toml @@ -32,7 +32,6 @@ Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" -StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" Tar = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" @@ -41,14 +40,14 @@ ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [weakdeps] GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231" [extensions] GraphMakieExt = ["Makie", "GraphMakie"] MakieExt = "Makie" -GraphPlotExt = ["Plots", "GraphPlot"] -PlotsExt = "Plots" +GraphPlotExt = ["StatsPlots", "GraphPlot"] +PlotsExt = "StatsPlots" [compat] BlackBoxOptim = "0.5, 0.6" @@ -79,3 +78,5 @@ StatsFuns = "1.3" YAML = "0.4" ZipFile = "0.10" julia = "1" +StatsPlots = "0.15" +GraphPlot = "0.6" diff --git a/README.md b/README.md index 54452715..3e9c4d4e 100644 --- a/README.md +++ b/README.md @@ -89,11 +89,17 @@ The examples below are run from the [examples](https://github.com/ConnectedSyste ```julia using Statistics using CSV, DataFrames, YAML -using Plots +using StatsPlots using Streamfall +example_data_dir = joinpath(dirname(dirname(pathof(Streamfall))), "test/data") + # Load data file which holds observed streamflow, precipitation and PET data -obs_data = CSV.read("../test/data/cotter/climate/CAMELS-AUS_410730.csv", DataFrame; comment="#") +obs_data = CSV.read( + joinpath(example_data_dir, "cotter/climate/CAMELS-AUS_410730.csv"), + DataFrame; + comment="#" +) # 18808×8 DataFrame # Row │ year month day Date 410730_P 410730_PET 410730_max_T 410730_Q # │ Int64 Int64 Int64 Date Float64 Float64 Float64 Float64 @@ -213,17 +219,18 @@ model. ```julia using CSV, DataFrames, YAML -using Plots +using StatsPlots using Streamfall # Load a network from a file, providing a name for the network and the file path. # Creates a graph representation of the stream with associated metadata. -sn = load_network("Example Network", "../test/data/campaspe/campaspe_network.yml") +example_data_dir = joinpath(dirname(dirname(pathof(Streamfall))), "test/data") +sn = load_network("Example Network", joinpath(example_data_dir, "campaspe/campaspe_network.yml")) # Load climate data, in this case from a CSV file with data for all nodes. # Indicate which columns are precipitation and evaporation data based on partial identifiers -climate = Climate("../test/data/campaspe/climate/climate.csv", "_rain", "_evap") +climate = Climate(joinpath(example_data_dir, "campaspe/climate/climate.csv"), "_rain", "_evap") # This runs an entire stream network @info "Running an example stream..." diff --git a/ext/GraphPlotExt/GraphPlotExt.jl b/ext/GraphPlotExt/GraphPlotExt.jl index 2ce4eda4..3240a990 100644 --- a/ext/GraphPlotExt/GraphPlotExt.jl +++ b/ext/GraphPlotExt/GraphPlotExt.jl @@ -1,6 +1,6 @@ module GraphPlotExt -using Plots, GraphPlot +using StatsPlots, GraphPlot using Streamfall import Streamfall: vertices import Streamfall: StreamfallNetwork diff --git a/ext/MakieExt/MakieExt.jl b/ext/MakieExt/MakieExt.jl index e36bfd6f..b315fe52 100644 --- a/ext/MakieExt/MakieExt.jl +++ b/ext/MakieExt/MakieExt.jl @@ -71,7 +71,7 @@ function Streamfall.Viz.quickplot( show_range = burn_in:last_e return Streamfall.Viz.quickplot(obs[show_range], sim[show_range], date[show_range], label, log; metric=metric) end -function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, xticklabels::Vector, label="Modeled", log=false; metric=Streamfall.mKGE) +function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, xticklabels::Vector; label="Modeled", log=false, metric=Streamfall.mKGE) @assert length(xticklabels) == length(obs) || "x-axis tick label length and observed lengths do not match!" @assert length(xticklabels) == length(sim) || "x-axis tick label length and simulated lengths do not match!" diff --git a/ext/PlotsExt/PlotsExt.jl b/ext/PlotsExt/PlotsExt.jl index 5d1d8269..aba90cab 100644 --- a/ext/PlotsExt/PlotsExt.jl +++ b/ext/PlotsExt/PlotsExt.jl @@ -1,13 +1,14 @@ module PlotsExt -using Plots, StatsPlots -using Plots.Measures +using StatsPlots +using StatsPlots.Plots.Measures using DataFrames, Dates, Statistics, Distributions, LaTeXStrings import Bootstrap: bootstrap, BalancedSampling using Streamfall import Streamfall: Analysis.TemporalCrossSection +import Streamfall.Viz: symlog function Streamfall.Viz.plot(node::NetworkNode, climate::Climate) return Plots.plot( @@ -50,23 +51,23 @@ function Streamfall.Viz.quickplot(node::NetworkNode, climate::Climate) @assert length(all_dates) == length(node.outflow) || "Date length and result lengths do not match!" - fig = plot(all_dates, node.outflow) + fig = Plots.plot(all_dates, node.outflow) return fig end -function Streamfall.Viz.quickplot(obs, node::NetworkNode, climate::Climate, label="", log=false; burn_in=1, limit=nothing, metric=Streamfall.mKGE) - return Streamfall.Viz.quickplot(obs, node.outflow, climate, label, log; burn_in=burn_in, limit=limit, metric=metric) +function Streamfall.Viz.quickplot(obs::Vector, node::NetworkNode, climate::Climate; label::String="", log::Bool=false, burn_in=1, limit=nothing, metric=Streamfall.mKGE) + return Streamfall.Viz.quickplot(obs, node.outflow, climate; label, log, burn_in, limit, metric) end -function Streamfall.Viz.quickplot(obs::DataFrame, sim::Vector, climate::Climate, label="", log=false; burn_in=1, limit=nothing, metric=Streamfall.mKGE) +function Streamfall.Viz.quickplot(obs::DataFrame, sim::Vector, climate::Climate; label::String="", log::Bool=false, burn_in=1, limit=nothing, metric=Streamfall.mKGE) return Streamfall.Viz.quickplot(Matrix(obs[:, Not("Date")])[:, 1], sim, climate, label, log; burn_in, limit, metric) end -function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, climate::Climate, label="", log=false; burn_in=1, limit=nothing, metric=Streamfall.mKGE) +function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, climate::Climate; label::String="", log::Bool=false, burn_in=1, limit=nothing, metric=Streamfall.mKGE) date = timesteps(climate) last_e = !isnothing(limit) ? limit : lastindex(obs) show_range = burn_in:last_e - return quickplot(obs[show_range], sim[show_range], date[show_range], label, log; metric=metric) + return quickplot(obs[show_range], sim[show_range], date[show_range]; label, log, metric) end -function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, xticklabels::Vector, label="Modeled", log=false; metric=Streamfall.mKGE) +function Streamfall.Viz.quickplot(obs::Vector, sim::Vector, xticklabels::Vector; label="Modeled", log=false, metric=Streamfall.mKGE) @assert length(xticklabels) == length(obs) || "x-axis tick label length and observed lengths do not match!" @assert length(xticklabels) == length(sim) || "x-axis tick label length and simulated lengths do not match!" @@ -167,7 +168,7 @@ function Streamfall.Viz.temporal_cross_section( logscale = [:log, :log10] tmp = nothing - xsect_res = TemporalCrossSection(dates, obs, period) + xsect_res = TemporalCrossSection(dates, obs; period) if :yscale in arg_keys || :yaxis in arg_keys tmp = (:yscale in arg_keys) ? kwargs[:yscale] : kwargs[:yaxis] @@ -178,7 +179,7 @@ function Streamfall.Viz.temporal_cross_section( # Format function for y-axis tick labels (e.g., 10^x) format_func = y -> (y != 0) ? L"%$(Int(round(sign(y)) * 10))^{%$(round(abs(y), digits=1))}" : L"0" - log_xsect_res = TemporalCrossSection(dates, log_obs, period) + log_xsect_res = TemporalCrossSection(dates, log_obs; period) target = log_xsect_res.cross_section else target = xsect_res.cross_section @@ -273,7 +274,7 @@ function Streamfall.Viz.temporal_cross_section( logscale = [:log, :log10] tmp = nothing - xsect_res = TemporalCrossSection(dates, obs, sim, period) + xsect_res = TemporalCrossSection(dates, obs, sim; period) target = xsect_res.cross_section if :yscale in arg_keys || :yaxis in arg_keys @@ -286,7 +287,7 @@ function Streamfall.Viz.temporal_cross_section( # Format function for y-axis tick labels (e.g., 10^x) format_func = y -> (y != 0) ? L"%$(Int(round(sign(y)) * 10))^{%$(round(abs(y), digits=1))}" : L"0" - log_xsect_res = TemporalCrossSection(dates, log_obs, log_sim, period) + log_xsect_res = TemporalCrossSection(dates, log_obs, log_sim; period) target = log_xsect_res.cross_section end end diff --git a/src/Analysis/uncertainty.jl b/src/Analysis/uncertainty.jl index da9e4bea..74035183 100644 --- a/src/Analysis/uncertainty.jl +++ b/src/Analysis/uncertainty.jl @@ -4,23 +4,23 @@ using OrderedCollections import ..ME, ..MAE -mutable struct TemporalCrossSection{A<:Dates.Date, B<:Real, D<:OrderedDict} +mutable struct TemporalCrossSection{A<:Dates.Date,B<:Real,D<:OrderedDict} dates::Array{A} ts::Array{B} subperiods::D - cross_section::Union{DataFrame, Nothing} + cross_section::Union{DataFrame,Nothing} end -function TemporalCrossSection(dates::Array, ts::Array, period::Function=monthday) +function TemporalCrossSection(dates::Array, ts::Array; period::Function=monthday) df = DataFrame(Date=dates, data=ts) sp = sort(unique(period.(dates))) # Remove leap days (when using monthday) - deleteat!(sp, findall(x -> x == (2,29), sp)) + deleteat!(sp, findall(x -> x == (2, 29), sp)) res = OrderedDict( - sp_i => df[period.(df.Date) .== [sp_i], :].data for sp_i in sp + sp_i => df[period.(df.Date).==[sp_i], :].data for sp_i in sp ) return TemporalCrossSection(dates, ts, res, nothing) @@ -30,10 +30,12 @@ end """ TemporalCrossSection constructor that handles two time series """ -function TemporalCrossSection(dates::Array, obs::Array, sim::Array, - period::Function=monthday) +function TemporalCrossSection( + dates::Array, obs::Array, sim::Array; + period::Function=monthday +) Y = ME.(obs, sim) - return TemporalCrossSection(dates, Y, period) + return TemporalCrossSection(dates, Y; period) end @@ -56,7 +58,7 @@ function cross_section(tr::TemporalCrossSection) end cols = [:abs_min, :lower_95, :lower_75, :median, - :upper_75, :upper_95, :abs_max, :mean, :std] + :upper_75, :upper_95, :abs_max, :mean, :std] x_section = DataFrame(boxframe, cols) x_section[:, :subperiod] = collect(keys(sp)) diff --git a/src/Streamfall.jl b/src/Streamfall.jl index 9b46901c..8b9adfdb 100644 --- a/src/Streamfall.jl +++ b/src/Streamfall.jl @@ -4,8 +4,6 @@ using Statistics using Graphs, MetaGraphs, Distributed, DataFrames -const MODPATH = @__DIR__ - include("Network.jl") include("Nodes/Node.jl") include("Climate.jl") diff --git a/src/viz/Viz.jl b/src/viz/Viz.jl index b9075715..3b017c59 100644 --- a/src/viz/Viz.jl +++ b/src/viz/Viz.jl @@ -12,6 +12,16 @@ function plot_residuals end function save_figure end function save_figure! end + +"""Symmetrical log values. + +https://kar.kent.ac.uk/32810/2/2012_Bi-symmetric-log-transformation_v5.pdf +https://discourse.julialang.org/t/symmetrical-log-plot/45709/3 +""" +function symlog(y) + return sign.(y) .* log10.(1.0 .+ abs.(y)) +end + export plot, plot!, quickplot export temporal_cross_section, plot_residuals export save_figure, save_figure! diff --git a/src/viz/_plotting.jl b/src/viz/_plotting.jl index 29664af5..5f82c9a7 100644 --- a/src/viz/_plotting.jl +++ b/src/viz/_plotting.jl @@ -1,6 +1,6 @@ -using Plots, StatsPlots -using Plots.Measures -import Plots: plot, plot! +using StatsPlots +using StatsPlots.Plots +import StatsPlots: plot, plot! using DataFrames, Dates, Statistics, Distributions, LaTeXStrings import Bootstrap: bootstrap, BalancedSampling @@ -142,7 +142,7 @@ function temporal_cross_section(dates, obs; logscale = [:log, :log10] tmp = nothing - xsect_res = TemporalCrossSection(dates, obs, period) + xsect_res = TemporalCrossSection(dates, obs; period) if :yscale in arg_keys || :yaxis in arg_keys tmp = (:yscale in arg_keys) ? kwargs[:yscale] : kwargs[:yaxis] @@ -153,7 +153,7 @@ function temporal_cross_section(dates, obs; # Format function for y-axis tick labels (e.g., 10^x) format_func = y -> (y != 0) ? L"%$(Int(round(sign(y)) * 10))^{%$(round(abs(y), digits=1))}" : L"0" - log_xsect_res = TemporalCrossSection(dates, log_obs, period) + log_xsect_res = TemporalCrossSection(dates, log_obs; period) target = log_xsect_res.cross_section else target = xsect_res.cross_section @@ -250,7 +250,7 @@ function temporal_cross_section( logscale = [:log, :log10] tmp = nothing - xsect_res = TemporalCrossSection(dates, obs, sim, period) + xsect_res = TemporalCrossSection(dates, obs, sim; period) target = xsect_res.cross_section if :yscale in arg_keys || :yaxis in arg_keys @@ -263,7 +263,7 @@ function temporal_cross_section( # Format function for y-axis tick labels (e.g., 10^x) format_func = y -> (y != 0) ? L"%$(Int(round(sign(y)) * 10))^{%$(round(abs(y), digits=1))}" : L"0" - log_xsect_res = TemporalCrossSection(dates, log_obs, log_sim, period) + log_xsect_res = TemporalCrossSection(dates, log_obs, log_sim; period) target = log_xsect_res.cross_section end end From 5915b5f5a9db9e3ac5f0423c3a10bb9d11ced27b Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Fri, 25 Apr 2025 23:16:29 +1000 Subject: [PATCH 12/13] Update tests --- src/Nodes/Ensembles/WeightedEnsembleNode.jl | 18 ++++---- test/test_calibration.jl | 6 +-- test/test_data_op.jl | 2 +- test/test_networks.jl | 47 ++++++++++++++++++--- 4 files changed, 56 insertions(+), 17 deletions(-) diff --git a/src/Nodes/Ensembles/WeightedEnsembleNode.jl b/src/Nodes/Ensembles/WeightedEnsembleNode.jl index eed17728..bf501351 100644 --- a/src/Nodes/Ensembles/WeightedEnsembleNode.jl +++ b/src/Nodes/Ensembles/WeightedEnsembleNode.jl @@ -1,10 +1,12 @@ -Base.@kwdef mutable struct WeightedEnsembleNode{N<:NetworkNode, P, A<:Real} <: EnsembleNode +import ..Analysis: TemporalCrossSection + +Base.@kwdef mutable struct WeightedEnsembleNode{N<:NetworkNode,P,A<:Real} <: EnsembleNode name::String area::A instances::Array{N} = NetworkNode[] weights::Array{P} = [ - Param(1.0; bounds=[0.0,1.0]) + Param(1.0; bounds=[0.0, 1.0]) ] # Default to normalized weighted sum @@ -25,8 +27,8 @@ function create_node( end function WeightedEnsembleNode(nodes::Vector{<:NetworkNode}; weights::Vector{Float64}, - bounds::Union{Nothing, Vector{Tuple}}=nothing, - comb_method::Union{Nothing, Function}=nothing)::WeightedEnsembleNode + bounds::Union{Nothing,Vector{Tuple}}=nothing, + comb_method::Union{Nothing,Function}=nothing)::WeightedEnsembleNode if isnothing(bounds) num_nodes = length(nodes) @assert(length(weights) == num_nodes, "Number of nodes do not match provided number of weights") @@ -39,7 +41,7 @@ function WeightedEnsembleNode(nodes::Vector{<:NetworkNode}; weights::Vector{Floa p_weights = [Param(w; bounds=b) for (w, b) in zip(weights, bounds)] n1 = nodes[1] - tmp = WeightedEnsembleNode{NetworkNode, Param, Float64}(; + tmp = WeightedEnsembleNode{NetworkNode,Param,Float64}(; name=n1.name, area=n1.area, instances=nodes, @@ -157,12 +159,12 @@ function calibrate!( ensemble::WeightedEnsembleNode, climate::Climate, calib_data::DataFrame, - metric::Union{AbstractDict{String, C}, C}; + metric::Union{AbstractDict{String,C},C}; kwargs... ) where {C<:Function} return invoke( calibrate!, - Tuple{NetworkNode, Climate, DataFrame, typeof(metric)}, + Tuple{NetworkNode,Climate,DataFrame,typeof(metric)}, ensemble, climate, calib_data, @@ -208,7 +210,7 @@ function apply_bias_correction( period=monthday ) where {T<:Real} dates = timesteps(climate) - tcs = TemporalCrossSection(dates, obs, ensemble.outflow) + tcs = TemporalCrossSection(dates, obs, ensemble.outflow; period) return ensemble.outflow .+ (-tcs.ts) end diff --git a/test/test_calibration.jl b/test/test_calibration.jl index 04198196..044aff40 100644 --- a/test/test_calibration.jl +++ b/test/test_calibration.jl @@ -3,7 +3,7 @@ using Dates, DataFrames, CSV, YAML using Streamfall -DATA_PATH = joinpath(@__DIR__, "data/hymod") +DATA_PATH = joinpath(dirname(dirname(pathof(Streamfall))), "test/data/hymod") @testset "Single node calibration" begin # Load and generate stream network @@ -13,8 +13,8 @@ DATA_PATH = joinpath(@__DIR__, "data/hymod") # Load climate data date_format = "YYYY-mm-dd" obs_data = CSV.File(joinpath(DATA_PATH, "leaf_river_data.csv"), - comment="#", - dateformat=date_format) |> DataFrame + comment="#", + dateformat=date_format) |> DataFrame # Column names must be the gauge name rename!(obs_data, ["leaf_river_outflow" => "leaf_river"]) diff --git a/test/test_data_op.jl b/test/test_data_op.jl index 9fd63222..5638dd3c 100644 --- a/test/test_data_op.jl +++ b/test/test_data_op.jl @@ -6,7 +6,7 @@ using Streamfall @testset "Data alignment" begin - here = @__DIR__ + here = joinpath(dirname(dirname(pathof(Streamfall))), "test") climate_data = joinpath(here, "data/campaspe/climate/climate_historic.csv") dam_data_loc = joinpath(here, "data/campaspe/gauges") diff --git a/test/test_networks.jl b/test/test_networks.jl index 5388b114..70f012a1 100644 --- a/test/test_networks.jl +++ b/test/test_networks.jl @@ -7,7 +7,7 @@ using DataFrames using Streamfall -DATA_PATH = joinpath(@__DIR__, "data/hymod") +DATA_PATH = joinpath(dirname(dirname(pathof(Streamfall))), "test/data/hymod") @testset "Network loading" begin @@ -71,7 +71,7 @@ end @testset "Reloading a network spec (HyMod)" begin - spec::Dict{String, Union{Dict{String, Any}, Any}} = Streamfall.extract_network_spec(sn) + spec::Dict{String,Union{Dict{String,Any},Any}} = Streamfall.extract_network_spec(sn) @test spec isa Dict @test haskey(spec, "leaf_river") || "Known node does not exist" @@ -97,7 +97,7 @@ end ) # Column name must match node name - rename!(obs_data, ["leaf_river_outflow"=>"leaf_river"]) + rename!(obs_data, ["leaf_river_outflow" => "leaf_river"]) climate_data = obs_data[:, ["Date", "leaf_river_P", "leaf_river_ET"]] climate = Climate(climate_data, "_P", "_ET") @@ -109,9 +109,46 @@ end @testset "Recursing IHACRESNode upstream" begin begin - include("../examples/run_nodes.jl") - # Ensure example does not error out + data_path = joinpath(dirname(dirname(pathof(Streamfall))), "test/data/campaspe/") + + # Load and generate stream network + sn = load_network("Example Network", joinpath(data_path, "campaspe_network.yml")) + + climate = Climate("../test/data/campaspe/climate/climate.csv", "_rain", "_evap") + + # Historic flows and dam level data + calib_data = CSV.read( + joinpath(data_path, "gauges", "outflow_and_level.csv"), + DataFrame; + comment="#" + ) + # Historic extractions from the dam + extraction_data = CSV.read( + joinpath(data_path, "gauges", "dam_extraction.csv"), + DataFrame; + comment="#" + ) + + @info "Running example stream..." + + dam_id, dam_node = sn["406000"] + Streamfall.run_node!(sn, dam_id, climate; extraction=extraction_data) + + # Extract data for comparison with 1-year burn-in period + dam_obs = calib_data[:, "406000"][366:end] + dam_sim = dam_node.level[366:end] + + nnse_score = Streamfall.NNSE(dam_obs, dam_sim) + nse_score = Streamfall.NSE(dam_obs, dam_sim) + rmse_score = Streamfall.RMSE(dam_obs, dam_sim) + + @info "Obj Func Scores:" rmse_score nnse_score nse_score + + nse = round(nse_score, digits=4) + rmse = round(rmse_score, digits=4) + + # Ensure example does not error out reset!(sn) run_basin!(sn, climate; extraction=extraction_data) @test Streamfall.RMSE(dam_obs, sn[dam_id].level[366:end]) < 2.5 From b93e291e179cc46fcb2be95318371d826e7fbc66 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Fri, 25 Apr 2025 23:21:44 +1000 Subject: [PATCH 13/13] Reorder includes --- src/Streamfall.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/Streamfall.jl b/src/Streamfall.jl index 8b9adfdb..501810e5 100644 --- a/src/Streamfall.jl +++ b/src/Streamfall.jl @@ -9,6 +9,8 @@ include("Nodes/Node.jl") include("Climate.jl") include("metrics.jl") include("calibration.jl") +include("Analysis/Analysis.jl") + include("Nodes/IHACRES/IHACRESNode.jl") # include("Nodes/IHACRES/IHACRESExpuhNode.jl") include("Nodes/GR4J/GR4JNode.jl") @@ -17,6 +19,9 @@ include("Nodes/SYMHYD/SYMHYDNode.jl") include("Nodes/DamNode.jl") include("Nodes/Ensembles/EnsembleNode.jl") +include("viz/Viz.jl") +include("viz/NetworkViz.jl") + function timestep_value(ts::Int, gauge_id::String, col_partial::String, dataset::DataFrame)::Float64 target_col = filter( @@ -271,10 +276,6 @@ function run_node!( ) end -include("Analysis/Analysis.jl") -include("viz/Viz.jl") -include("viz/NetworkViz.jl") - # Nodes export NetworkNode, GenericNode, GenericDirectNode export IHACRES, IHACRESNode, IHACRESBilinearNode, DamNode, Climate