From 54f909ee80645f798689a08c41b5847fbda2477d Mon Sep 17 00:00:00 2001 From: achiang17 Date: Fri, 6 Dec 2024 12:48:31 -0800 Subject: [PATCH] move total streamflow script to riverModels.jl --- examples/routing/gamma_IRF.jl | 22 +++++++----- examples/routing/mini_gamma_IRF.jl | 10 +++--- src/ClimaRivers.jl | 1 + src/Environments.jl | 2 +- src/MizurouteV1.jl | 6 +--- src/RiverModels.jl | 57 +++++++++++++++++++++++++++--- src/Routing.jl | 49 ++++++++++++++----------- 7 files changed, 104 insertions(+), 43 deletions(-) diff --git a/examples/routing/gamma_IRF.jl b/examples/routing/gamma_IRF.jl index e9fcebf..c649584 100644 --- a/examples/routing/gamma_IRF.jl +++ b/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 diff --git a/examples/routing/mini_gamma_IRF.jl b/examples/routing/mini_gamma_IRF.jl index 2ced7c0..cff8b06 100644 --- a/examples/routing/mini_gamma_IRF.jl +++ b/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) diff --git a/src/ClimaRivers.jl b/src/ClimaRivers.jl index 40dfdf5..3a017c2 100644 --- a/src/ClimaRivers.jl +++ b/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 diff --git a/src/Environments.jl b/src/Environments.jl index 66debcf..4736f33 100644 --- a/src/Environments.jl +++ b/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" diff --git a/src/MizurouteV1.jl b/src/MizurouteV1.jl index 43f02f5..ea297dc 100644 --- a/src/MizurouteV1.jl +++ b/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 \ No newline at end of file diff --git a/src/RiverModels.jl b/src/RiverModels.jl index 69bdc6e..3dea4a7 100644 --- a/src/RiverModels.jl +++ b/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") \ No newline at end of file diff --git a/src/Routing.jl b/src/Routing.jl index 857b2b0..52b2091 100644 --- a/src/Routing.jl +++ b/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)