Merge pull request #8 from CliMA/orad/river-model-infrastructure

WIP: Starting river-model infrastructure

.gitignore

@@ -22,3 +22,6 @@ docs/site/
 # committed for packages, but should be committed for applications that require a static
 # environment.
 Manifest.toml
+
+# Mini data set exception
+!mini_data/**/*.csv

Project.toml

@@ -3,5 +3,29 @@ uuid = "9363829b-25e6-4e43-a03e-985c67d112e3"
 authors = ["Mauricio Lima, Andrew Chiang, Katherine Deck, Ollie Dunbar"]
 version = "0.1.0"
 
+[deps]
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
 [compat]
-julia = "1.10"
+CSV = "0.10"
+DSP = "0.7"
+DataFrames = "1.7"
+Dates = "1"
+DocStringExtensions = "0.9"
+Glob = "1.3"
+JSON = "0.21"
+LinearAlgebra = "1"
+Random = "1"
+SpecialFunctions = "2.4"
+Statistics = "1"
+julia = "1.10"

examples/routing/compare_results.jl

@@ -0,0 +1,67 @@
+using CSV
+using DataFrames
+using Statistics
+
+# Define paths
+basin_path = "/groups/esm/achiang/ClimaRivers.jl/data/routing/simulations/simulations_lv05/simulation_gamma-IRF"
+channel_path = "/groups/esm/achiang/ClimaRivers.jl/data/routing/simulations/simulations_lv05/gamma_IRF"
+hillslope_path = "/groups/esm/achiang/ClimaRivers.jl/data/routing/simulations/simulations_lv05/gamma_IRF"
+output_path = "/groups/esm/achiang/ClimaRivers.jl/examples/routing/csv_files"  # Directory to save output CSVs
+
+# Function to check the streamflow relationship and save results
+function check_basin_streamflow(basin_id::String)
+    # Construct file paths
+    basin_file = joinpath(basin_path, "basin_$basin_id.csv")
+    channel_file = joinpath(channel_path, "channel_basin_$basin_id.csv")
+    hillslope_file = joinpath(hillslope_path, "hillslope_basin_$basin_id.csv")
+    output_file = joinpath(output_path, "comparison_basin_$basin_id.csv")
+
+    # Check if files exist
+    if !(isfile(basin_file) && isfile(channel_file) && isfile(hillslope_file))
+        println("One or more files for basin_id $basin_id are missing.")
+        return nothing
+    end
+
+    # Load data
+    basin_df = CSV.read(basin_file, DataFrame)
+    channel_df = CSV.read(channel_file, DataFrame)
+    hillslope_df = CSV.read(hillslope_file, DataFrame)
+
+    # Check if the date columns align
+    if !all(basin_df.date .== channel_df.date) ||
+       !all(basin_df.date .== hillslope_df.date)
+        println("Date mismatch for basin_id $basin_id.")
+        return nothing
+    end
+
+    # Compute the sum of hillslope and channel streamflows
+    summed_streamflow = hillslope_df.streamflow .+ channel_df.streamflow
+
+    # Compute the difference
+    diff = summed_streamflow .- basin_df.streamflow
+
+    # Create a new DataFrame with the results
+    result_df = DataFrame(
+        date = basin_df.date,
+        summed_streamflow = summed_streamflow,
+        basin_streamflow = basin_df.streamflow,
+        difference = diff,
+    )
+
+    # Save the DataFrame to a new CSV
+    CSV.write(output_file, result_df)
+    println("Results saved to $output_file")
+
+    # Return basic statistics
+    max_diff = maximum(abs.(diff))
+    mean_diff = mean(abs.(diff))
+    println("Basin ID: $basin_id")
+    println("Maximum absolute difference: $max_diff")
+    println("Mean absolute difference: $mean_diff")
+
+    return (max_diff, mean_diff)
+end
+
+# Example usage
+basin_id_to_check = "1051315930"  # Replace with the specific basin_id you want to check
+check_basin_streamflow(basin_id_to_check)

examples/routing/gamma_IRF.jl

@@ -0,0 +1,54 @@
+using ClimaRivers
+using JSON, Dates
+
+# build hillslope model
+hillslope = MizurouteHillslopeV1{Float64}()
+
+# build channel model
+channel = MizurouteChannelV1{Float64}()
+
+# build model
+river_model = HillslopeChannelRiverModel(hillslope, channel)
+
+# build environment
+data_file_path = joinpath(@__DIR__, "..", "..", "data", "routing")
+
+# build static environment
+@info "reading data files from $(data_file_path)"
+graph_dict =
+    JSON.parsefile(joinpath(data_file_path, "graphs", "graph_lv05.json"))
+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 =
+    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)
+
+env = Environment(static_env, dynamic_env)
+
+## evolutionary model, evolving a state over time
+model_types = ["instant"]
+model_type = model_types[1]
+
+start_date = Date("1996-01-01", "yyyy-mm-dd")
+end_date = Date("2014-12-31", "yyyy-mm-dd")
+dates = collect(start_date:Day(1):end_date)
+
+# 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 = RiverState(hillslope_data, channel_data)
+
+if model_type == "instant"
+    ## full-timeseries model, predicts all states at once
+    compute_streamflow!(river_state, river_model, env, start_date, end_date)
+end

examples/routing/mini_gamma_IRF.jl

@@ -0,0 +1,53 @@
+using ClimaRivers
+using JSON, Dates
+
+# build hillslope model
+hillslope = MizurouteHillslopeV1{Float64}()
+
+# build channel model
+channel = MizurouteChannelV1{Float64}()
+
+# build model
+river_model = HillslopeChannelRiverModel(hillslope, channel)
+
+# build environment
+data_file_path = joinpath(@__DIR__, "..", "..", "mini_data", "routing")
+
+# build static environment
+@info "reading data files from $(data_file_path)"
+graph_dict =
+    JSON.parsefile(joinpath(data_file_path, "graphs", "graph_lv05.json"))
+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 =
+    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)
+
+env = Environment(static_env, dynamic_env)
+
+## evolutionary model, evolving a state over time
+model_types = ["instant"]
+model_type = model_types[1]
+
+start_date = Date("1996-01-01", "yyyy-mm-dd")
+end_date = Date("2014-12-31", "yyyy-mm-dd")
+
+# 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)
+
+if model_type == "instant"
+    ## full-timeseries model, predicts all states at once
+    compute_streamflow!(river_state, river_model, env, start_date, end_date)
+end

mini_data/routing/attributes/attributes_lv05/attributes.csv

@@ -0,0 +1,6 @@
+HYBAS_ID,area,DIST_MAIN
+1051429650,16915.2,1593.8
+1051429580,21259.1,1594.0
+1051435100,658.7,1550.6
+1051460430,7361.1,1372.2
+1051435110,10896.7,1550.6

mini_data/routing/graphs/graph_lv05.json

@@ -0,0 +1 @@
+{"1051429650":[], "1051429580":[], "1051435110":[],"1051435100":[1051429650,1051429580], "1051460430":[1051435100,1051435110]}

mini_data/routing/routing_lvs/routing_lvs_lv05/all_basin_ids.txt

@@ -0,0 +1,5 @@
+1051429650
+1051429580
+1051435110
+1051435100
+1051460430