move total streamflow script to riverModels.jl

examples/routing/gamma_IRF.jl

@@ -10,21 +10,27 @@ channel = MizurouteChannelV1{Float64}()
 # build model
 river_model = HillslopeChannelRiverModel(hillslope, channel)
 
+# build environment
+data_file_path = joinpath(@__DIR__,"..","..","data","routing")
+
 # build static environment
-routing_lv_basins_dir = "/groups/esm/achiang/ClimaRivers.jl/data/routing/routing_lvs/routing_lvs_lv05"
-attributes_dir = "/groups/esm/achiang/ClimaRivers.jl/data/routing/attributes/attributes_lv05"
+@info "reading data files from $(data_file_path)"
 graph_dict = JSON.parsefile(
-    "/groups/esm/achiang/ClimaRivers.jl/data/routing/graphs/graph_lv05.json",
+    joinpath(data_file_path, "graphs","graph_lv05.json"),
 )
-static_env = StaticEnvironment(routing_lv_basins_dir, attributes_dir, graph_dict)
+basins_dir = joinpath(data_file_path, "routing_lvs","routing_lvs_lv05")
+attributes_dir = joinpath(data_file_path, "attributes","attributes_lv05")
+static_env = StaticEnvironment(basins_dir, attributes_dir, graph_dict)
 
 # build dynamic environment
-forcing_timeseries_dir = "/groups/esm/achiang/ClimaRivers.jl/data/routing/timeseries/timeseries_lv05"
-output_dir = "/groups/esm/achiang/ClimaRivers.jl/data/routing/simulations/simulations_lv05/gamma_IRF"
+forcing_timeseries_dir = joinpath(data_file_path, "timeseries","timeseries_lv05")
+output_dir = joinpath(data_file_path, "simulations","simulations_lv05","gamma_IRF")
+@info "creating output"
+if !isdir(output_dir)
+    mkpath(output_dir)
+end
 dynamic_env = DynamicEnvironment(forcing_timeseries_dir, output_dir)
 
-# build environment
-
 env = Environment(static_env, dynamic_env)
 
 ## evolutionary model, evolving a state over time

examples/routing/mini_gamma_IRF.jl

@@ -18,9 +18,9 @@ data_file_path = joinpath(@__DIR__,"..","..","mini_data","routing")
 graph_dict = JSON.parsefile(
     joinpath(data_file_path, "graphs","graph_lv05.json"),
 )
-routing_lv_basins_dir = joinpath(data_file_path, "routing_lvs","routing_lvs_lv05")
+basins_dir = joinpath(data_file_path, "routing_lvs","routing_lvs_lv05")
 attributes_dir = joinpath(data_file_path, "attributes","attributes_lv05")
-static_env = StaticEnvironment(routing_lv_basins_dir, attributes_dir, graph_dict)
+static_env = StaticEnvironment(basins_dir, attributes_dir, graph_dict)
 
 # build dynamic environment
 forcing_timeseries_dir = joinpath(data_file_path, "timeseries","timeseries_lv05")
@@ -40,9 +40,9 @@ model_type = model_types[1]
 start_date = Date("1996-01-01", "yyyy-mm-dd")
 end_date = Date("2014-12-31", "yyyy-mm-dd")
 
-# streamflow = zeros(dates,basins)
-hillslope_data = zeros(10, 10)  # Replace with actual data once implemented
-channel_data = zeros(10, 10)   # Replace with actual data once implemented
+# River state loaded into csv files currently, placehodler variable
+hillslope_data = zeros(10, 10)
+channel_data = zeros(10, 10)
 
 river_state = RiverState(hillslope_data, channel_data)
 

src/ClimaRivers.jl

@@ -9,4 +9,5 @@ include("Environments.jl")
 include("States.jl")
 include("RiverModels.jl")
 include("Routing.jl")
+include("MizurouteV1.jl")
 end # module ClimaRivers

src/Environments.jl

@@ -4,7 +4,7 @@ export StaticEnvironment, DynamicEnvironment, Environment
 # Static Data Objects
 struct StaticEnvironment
     "Directory containing list of all basins"
-    routing_lv_basins_dir::String
+    basins_dir::String
     "Directory containing attribute csv"
     attributes_dir::String
     "Basin mapping dictionary"

src/MizurouteV1.jl

@@ -1,6 +1,5 @@
 export MizurouteHillslopeV1, MizurouteChannelV1
 
-
 # Contains the Mizuroute model    
 struct MizurouteHillslopeV1{FT <: AbstractFloat} <: AbstractHillslopeModel
     "Shape factor, a (adjusted)"
@@ -35,7 +34,4 @@ function MizurouteChannelV1{FT}() where {FT <: AbstractFloat}
     diffusivity = FT(800 * 86400)
     t_max = FT(120.0)
     return MizurouteChannelV1(wave_velocity, diffusivity, t_max)
-end
-
-
-function update_state_from_channel!(river_state, hillslope_model, environment) end
+end

src/RiverModels.jl

@@ -1,5 +1,7 @@
 # Contains the abstract river model types, structs and methods
-export HillslopeChannelRiverModel
+export HillslopeChannelRiverModel, get_basin_list
+
+using CSV, DataFrames, Statistics
 
 #general model level
 abstract type AbstractRiverModel end
@@ -16,11 +18,58 @@ struct HillslopeChannelRiverModel{
     channel_model::CH
 end
 
+# Auxiliary function
+## function for reading basins from txt file into Vector{Int}
+function get_basin_list(basins_file::String)
+    basins_list = Int64[]
+    file = open(joinpath(basins_file))
+    for line in eachline(file)
+        push!(basins_list, parse(Int64, line))
+    end
+    close(file)
+
+    return basins_list
+end
+
+function calculate_streamflow(
+    river_state::RS,
+    env::E,
+) where {RS <: RiverState, E <: Environment}
 
-function calculate_streamflow(river_state::RS) where {RS <: RiverState}
+    output_dir = env.dynamic_env.output_dir
 
-    return RS.hillslope + RS.channel
+    all_basin_ids = get_basin_list(
+        joinpath(env.static_env.basins_dir, "all_basin_ids.txt"),
+    )
 
+    for basin_id in all_basin_ids
+        channel_file = joinpath(output_dir, "channel_basin_$basin_id.csv")
+        hillslope_file = joinpath(output_dir, "hillslope_basin_$basin_id.csv")
+        total_file = joinpath(output_dir, "basin_$basin_id.csv")
+
+        # Check if files exist
+        if !(isfile(channel_file) && isfile(hillslope_file))
+            println("One or more files for basin_id $basin_id are missing.")
+            return nothing
+        end
+
+        channel_df = CSV.read(channel_file, DataFrame)
+        hillslope_df = CSV.read(hillslope_file, DataFrame)
+
+        # Compute the sum of hillslope and channel streamflows
+        total_streamflow = hillslope_df.streamflow .+ channel_df.streamflow
+
+        # Create a new DataFrame with the results
+        result_df = DataFrame(
+            date = channel_df.date,
+            total_streamflow = total_streamflow,
+        )
+
+        # Save the DataFrame to a new CSV
+        CSV.write(total_file, result_df)
+        println("Results saved to $total_file")
+    end
 end
+
 # Specific hillslope-channel models loaded here:
-include("MizurouteV1.jl")
+# include("MizurouteV1.jl")

src/Routing.jl

@@ -12,7 +12,17 @@ function compute_streamflow!(
     start_date::Date,
     end_date::Date,
 ) where {RS <: RiverState, HCM <: HillslopeChannelRiverModel, E <: Environment}
+    update_state!(river_state, river_model, environment, start_date, end_date)
+    return calculate_streamflow(river_state, environment)
+end
 
+function update_state!(
+    river_state::RS,
+    river_model::HCM,
+    environment::E,
+    start_date::Date,
+    end_date::Date,
+) where {RS <: RiverState, HCM <: HillslopeChannelRiverModel, E <: Environment}
     update_state_from_hillslope!(
         river_state,
         river_model.hillslope_model,
@@ -27,22 +37,20 @@ function compute_streamflow!(
         start_date,
         end_date,
     )
-    # return calculate_streamflow(river_state)
 end
 
-
 # Auxiliary function
 ## function for reading basins from txt file into Vector{Int}
-function read_routing_lv_basins(routing_lv_file::String)
-    routing_lv_basins = Int64[]
-    file = open(joinpath(routing_lv_file))
-    for line in eachline(file)
-        push!(routing_lv_basins, parse(Int64, line))
-    end
-    close(file)
+# function get_basin_list(basins_file::String)
+#     basins_file = Int64[]
+#     file = open(joinpath(basins_file))
+#     for line in eachline(file)
+#         push!(basins_file, parse(Int64, line))
+#     end
+#     close(file)
 
-    return routing_lv_basins
-end
+#     return basins_file
+# end
 
 # only for one routing lv
 function update_state_from_hillslope!(
@@ -62,7 +70,7 @@ function update_state_from_hillslope!(
     output_dir = env.dynamic_env.output_dir
     dates = collect(start_date:Day(1):end_date)
 
-    routing_lv_basins_dir = env.static_env.routing_lv_basins_dir
+    basins_dir = env.static_env.basins_dir
     attributes_df = CSV.read(
         joinpath(env.static_env.attributes_dir, "attributes.csv"),
         DataFrame,
@@ -76,11 +84,11 @@ function update_state_from_hillslope!(
         t in 0:(t_max - 1)
     ]
 
-    all_basin_ids = read_routing_lv_basins(
-        joinpath(routing_lv_basins_dir, "all_basin_ids.txt"),
+    all_basin_ids = get_basin_list(
+        joinpath(basins_dir, "all_basin_ids.txt"),
     )
-    # make new txt file 
-    for basin_id in all_basin_ids # routing_lv_basins
+
+    for basin_id in all_basin_ids
         timeseries_df =
             CSV.read(joinpath(forcing_timeseries_dir, "basin_$basin_id.csv"), DataFrame)
         filtered_df =
@@ -143,7 +151,7 @@ function update_state_from_channel!(
     forcing_timeseries_dir = env.dynamic_env.forcing_timeseries_dir
 
     graph_dict = env.static_env.graph_dict
-    routing_lv_basins_dir = env.static_env.routing_lv_basins_dir
+    basins_dir = env.static_env.basins_dir
     attributes_df = CSV.read(
         joinpath(env.static_env.attributes_dir, "attributes.csv"),
         DataFrame,
@@ -152,8 +160,8 @@ function update_state_from_channel!(
     C, D = channel_model.wave_velocity, channel_model.diffusivity
     t_max = channel_model.t_max
 
-    all_basin_ids = read_routing_lv_basins(
-        joinpath(routing_lv_basins_dir, "all_basin_ids.txt"),
+    all_basin_ids = get_basin_list(
+        joinpath(basins_dir, "all_basin_ids.txt"),
     )
 
     # Iterate over basins in the given routing level
@@ -185,7 +193,8 @@ function update_state_from_channel!(
                 ][1] -
                 attributes_df[attributes_df.HYBAS_ID .== basin_id, :DIST_MAIN][1]
             dist *= km_to_m
-            # Sum the runoff and sub-surface runoff columns
+
+            # Get upstream streamflow
             up_q = up_timeseries_df[:, :streamflow]
             up_q = reshape(up_q, :, 1)