diff --git a/docs/source/notebooks/N1_Demand_WaterRights_File_Modification.rst b/docs/source/notebooks/N1_Demand_WaterRights_File_Modification.rst index 6302a1d..43a7fee 100644 --- a/docs/source/notebooks/N1_Demand_WaterRights_File_Modification.rst +++ b/docs/source/notebooks/N1_Demand_WaterRights_File_Modification.rst @@ -58,10 +58,10 @@ conditions. import pickle from string import Template import subprocess - + import matplotlib.pyplot as plt import numpy as np - import pandas as pd + import pandas as pd import statemodify as stm .. container:: alert alert-block alert-info @@ -79,16 +79,16 @@ conditions. # statemod directory statemod_dir = "/usr/src/statemodify/statemod_gunnison_sjd" - + # root directory of statemod data for the target basin root_dir = os.path.join(statemod_dir, "src", "main", "fortran") - + # home directory of notebook instance home_dir = os.path.dirname(os.getcwd()) - + # path to the statemod executable statemod_exe = os.path.join(root_dir, "statemod-17.0.3-gfortran-lin-64bit-o3") - + # data directory and root name for the target basin data_dir = os.path.join( home_dir, @@ -97,25 +97,25 @@ conditions. "sj2015_StateMod_modified", "StateMod" ) - + # directory to the target basin input files with root name for the basin basin_path = os.path.join(data_dir, "sj2015B") - + # scenarios output directory scenarios_dir_ddm = os.path.join(data_dir, "scenarios_ddm") scenarios_dir_ddr = os.path.join(data_dir, "scenarios_ddr") - + # parquet files output directory parquet_dir_ddm = os.path.join(data_dir, "parquet_ddm") parquet_dir_ddr = os.path.join(data_dir, "parquet_ddr") - + # path to ddm and ddr template file ddm_template_file = os.path.join( home_dir, "data", "sj2015B_template_ddm.rsp" ) - + ddr_template_file = os.path.join( home_dir, "data", @@ -131,19 +131,19 @@ conditions. .. parsed-literal:: - Startup log file for messages to this point: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/sj2015B.rsp + Startup log file for messages to this point: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/sj2015B.rsp Closing startup log file: statem.log - Opening dataset log file: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/sj2015B.log + Opening dataset log file: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/sj2015B.log ________________________________________________________________________ - - StateMod - State of Colorado - Water Supply Planning Model - - Version: 17.0.3 + + StateMod + State of Colorado - Water Supply Planning Model + + Version: 17.0.3 Last revision date: 2021/09/12 - + ________________________________________________________________________ - + Subroutine Execut Subroutine Datinp @@ -151,7 +151,7 @@ conditions. ________________________________________________________________________ Execut; Successful Termination - Statem; See detailed messages in dataset log file: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/sj2015B.log + Statem; See detailed messages in dataset log file: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/sj2015B.log Stop 0 Once StateMod has run successfully, we can now extract user shortages @@ -167,16 +167,16 @@ and saved. .. code:: ipython3 - #Extract shortages using statemodify convert_xdd() function - - # create a directory to store the historical shortages + #Extract shortages using statemodify convert_xdd() function + + # create a directory to store the historical shortages output_dir = os.path.join(data_dir, "historic_shortages") - + # create a directory to store the new files in if it does not exist output_directory = os.path.join(data_dir, "historic_shortages") if not os.path.exists(output_directory): os.makedirs(output_directory) - + stm.xdd.convert_xdd( # path to a directory where output .parquet files should be written output_path=output_dir, @@ -212,11 +212,11 @@ and saved. .dataframe tbody tr th:only-of-type { vertical-align: middle; } - + .dataframe tbody tr th { vertical-align: top; } - + .dataframe thead th { text-align: right; } @@ -529,11 +529,11 @@ We can then take these shortages and plot them for our list of users. .. code:: ipython3 fig, ax = plt.subplots() - + for name, group in data.groupby('structure_id'): ax.scatter( group['year'], group['shortage_total'], label=name) - + plt.xlabel("Year") plt.ylabel("Shortage (AF)") plt.legend() @@ -585,29 +585,29 @@ the ``input_files`` directory. "ids": ["2900501", "2900519","2900555"], "bounds": [0.5, 1.5] } - + output_directory = output_dir = os.path.join(data_dir, "input_files") - + scenario = "1" - + # the number of samples you wish to generate n_samples = 2 - + # seed value for reproducibility if so desired seed_value = 1 - + # number of rows to skip in file after comment skip_rows = 1 - + # name of field to query query_field = "id" - + # number of jobs to launch in parallel; -1 is all but 1 processor used n_jobs = -1 - + # basin to process basin_name = "San_Juan" - + # generate a batch of files using generated LHS stm.modify_ddm( modify_dict=setup_dict, @@ -669,53 +669,53 @@ scenario will take approximately 4 minutes. # set realization and sample realization = 1 sample = np.arange(0, 2, 1) - + # read RSP template with open(ddm_template_file) as template_obj: - + # read in file template_rsp = Template(template_obj.read()) - + for i in sample: - + # create scenario name scenario = f"S{i}_{realization}" - + # dictionary holding search keys and replacement values to update the template file d = {"DDM": f"../../input_files/sj2015B_{scenario}.ddm"} - + # update the template new_rsp = template_rsp.safe_substitute(d) - + # construct simulated scenario directory simulated_scenario_dir = os.path.join(scenarios_dir_ddm, scenario) if not os.path.exists(simulated_scenario_dir): os.makedirs(simulated_scenario_dir) - + # target rsp file rsp_file = os.path.join(simulated_scenario_dir, f"sj2015B_{scenario}.rsp") - + # write updated rsp file with open(rsp_file, "w") as f1: f1.write(new_rsp) - + # construct simulated basin path simulated_basin_path = os.path.join(simulated_scenario_dir, f"sj2015B_{scenario}") - + # run StateMod print(f"Running: {scenario}") os.chdir(simulated_scenario_dir) - + subprocess.call([statemod_exe, simulated_basin_path, "-simulate"]) - - #Save output to parquet files + + #Save output to parquet files print('creating parquet for ' + scenario) - + output_directory = os.path.join(parquet_dir_ddm+"/scenario/"+ scenario) - + if not os.path.exists(output_directory): os.makedirs(output_directory) - + stm.xdd.convert_xdd( output_path=output_directory, allow_overwrite=False, @@ -730,19 +730,19 @@ scenario will take approximately 4 minutes. .. parsed-literal:: Running: S0_1 - Startup log file for messages to this point: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/scenarios_ddm/S0_1/sj2015B_S0_1.rsp + Startup log file for messages to this point: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/scenarios_ddm/S0_1/sj2015B_S0_1.rsp Closing startup log file: statem.log - Opening dataset log file: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/scenarios_ddm/S0_1/sj2015B_S0_1.log + Opening dataset log file: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/scenarios_ddm/S0_1/sj2015B_S0_1.log ________________________________________________________________________ - - StateMod - State of Colorado - Water Supply Planning Model - - Version: 17.0.3 + + StateMod + State of Colorado - Water Supply Planning Model + + Version: 17.0.3 Last revision date: 2021/09/12 - + ________________________________________________________________________ - + Subroutine Execut Subroutine Datinp @@ -750,7 +750,7 @@ scenario will take approximately 4 minutes. ________________________________________________________________________ Execut; Successful Termination - Statem; See detailed messages in dataset log file: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/scenarios_ddm/S0_1/sj2015B_S0_1.log + Statem; See detailed messages in dataset log file: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/scenarios_ddm/S0_1/sj2015B_S0_1.log Stop 0 creating parquet for S0_1 @@ -791,19 +791,19 @@ respect to the shortages received in the baseline case. baseline=pd.read_parquet(data_dir+'/historic_shortages/sj2015B.parquet',engine='pyarrow') SOW_1=pd.read_parquet(parquet_dir_ddm+'/scenario/S0_1/sj2015B_S0_1.parquet',engine='pyarrow') SOW_2=pd.read_parquet(parquet_dir_ddm+'/scenario/S1_1/sj2015B_S1_1.parquet',engine='pyarrow') - + # Subtract shortages with respect to the baseline subset_df=pd.concat([baseline['year'],baseline['shortage_total'],SOW_1['shortage_total'],SOW_2['shortage_total']],axis=1) subset_df = subset_df.set_axis(['Year', 'Baseline', 'SOW_1','SOW_2'], axis=1) subset_df['SOW_1_diff'] = subset_df['SOW_1']-subset_df['Baseline'] subset_df['SOW_2_diff'] = subset_df['SOW_2']-subset_df['Baseline'] - + # Plot shortages fig, ax = plt.subplots() - + ax.scatter(subset_df['Year'], subset_df['SOW_1_diff'],label='Decreased Demand') ax.scatter(subset_df['Year'], subset_df['SOW_2_diff'],label='Increased Demand') - + plt.xlabel("Year") plt.ylabel("Shortage (AF)") plt.title("Change in Shortages from the Baseline") @@ -846,36 +846,36 @@ rank to 1. setup_dict = { # ids can either be 'struct' or 'id' values "ids": ["2900501"], - + # turn id on or off completely or for a given period # if 0 = off, 1 = on, YYYY = on for years >= YYYY, -YYYY = off for years > YYYY; see file header "on_off": [1], - + # apply rank of administrative order where 0 is lowest (senior) and n is highest (junior); None is no change "admin": [1], } - + output_directory = os.path.join(data_dir, "input_files") scenario = "1" - + # the number of samples you wish to generate n_samples = 1 - + # seed value for reproducibility if so desired seed_value = 1 - + # number of rows to skip in file after comment skip_rows = 0 - + # name of field to query query_field = "struct" - + # number of jobs to launch in parallel; -1 is all but 1 processor used n_jobs = -1 - + # basin to process basin_name = "San_Juan" - + # generate a batch of files using generated LHS stm.modify_ddr( modify_dict=setup_dict, @@ -907,53 +907,53 @@ using the .\ ``ddr`` template file. # set realization and sample realization = 1 sample = np.arange(0, 1, 1) - + # read RSP template with open(ddr_template_file) as template_obj: - + # read in file template_rsp = Template(template_obj.read()) - + for i in sample: - + # create scenario name scenario = f"S{i}_{realization}" - + # dictionary holding search keys and replacement values to update the template file d = {"DDR": f"../../input_files/sj2015B_{scenario}.ddr"} - + # update the template new_rsp = template_rsp.safe_substitute(d) - + # construct simulated scenario directory simulated_scenario_dir = os.path.join(scenarios_dir_ddr, scenario) if not os.path.exists(simulated_scenario_dir): os.makedirs(simulated_scenario_dir) - + # target rsp file rsp_file = os.path.join(simulated_scenario_dir, f"sj2015B_{scenario}.rsp") - + # write updated rsp file with open(rsp_file, "w") as f1: f1.write(new_rsp) - + # construct simulated basin path simulated_basin_path = os.path.join(simulated_scenario_dir, f"sj2015B_{scenario}") - + # run StateMod print(f"Running: {scenario}") os.chdir(simulated_scenario_dir) - + subprocess.call([statemod_exe, simulated_basin_path, "-simulate"]) - - #Save output to parquet files + + #Save output to parquet files print('creating parquet for ' + scenario) - + output_directory = os.path.join(parquet_dir_ddr+"/scenario/"+ scenario) - + if not os.path.exists(output_directory): os.makedirs(output_directory) - + stm.xdd.convert_xdd( output_path=output_directory, allow_overwrite=False, @@ -968,27 +968,27 @@ using the .\ ``ddr`` template file. .. parsed-literal:: Running: S0_1 - Startup log file for messages to this point: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/scenarios_ddr/S0_1/sj2015B_S0_1.rsp + Startup log file for messages to this point: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/scenarios_ddr/S0_1/sj2015B_S0_1.rsp Closing startup log file: statem.log - Opening dataset log file: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/scenarios_ddr/S0_1/sj2015B_S0_1.log + Opening dataset log file: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/scenarios_ddr/S0_1/sj2015B_S0_1.log ________________________________________________________________________ - - StateMod - State of Colorado - Water Supply Planning Model - - Version: 17.0.3 + + StateMod + State of Colorado - Water Supply Planning Model + + Version: 17.0.3 Last revision date: 2021/09/12 - + ________________________________________________________________________ - + Subroutine Execut Subroutine Datinp - + ... ________________________________________________________________________ Execut; Successful Termination - Statem; See detailed messages in dataset log file: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/scenarios_ddr/S0_1/sj2015B_S0_1.log + Statem; See detailed messages in dataset log file: /home/jovyan/data/sj2015_StateMod_modified/sj2015_StateMod_modified/StateMod/scenarios_ddr/S0_1/sj2015B_S0_1.log Stop 0 creating parquet for S0_1 @@ -1001,17 +1001,17 @@ with respect to the baseline shortages. # Read in raw parquet files baseline=pd.read_parquet(data_dir+'/historic_shortages/sj2015B.parquet',engine='pyarrow') SOW_1=pd.read_parquet(parquet_dir_ddr+ '/scenario/S0_1/sj2015B_S0_1.parquet',engine='pyarrow') - + # Subtract shortages with respect to the baseline subset_df=pd.concat([baseline['year'],baseline['shortage_total'],SOW_1['shortage_total']],axis=1) subset_df = subset_df.set_axis(['Year', 'Baseline', 'SOW_1'], axis=1) subset_df['diff']=subset_df['SOW_1']-subset_df['Baseline'] - + # Plot shortages fig, ax = plt.subplots() - + ax.scatter(subset_df['Year'], subset_df['diff']) - + plt.xlabel("Year") plt.ylabel("Shortage (AF)") plt.title("Change in Shortages from the Baseline") @@ -1052,4 +1052,3 @@ reservoir evaporation modification fuction. :: 2. modify_ddr() - diff --git a/docs/source/notebooks/N2_Evaporation_File_Modification.rst b/docs/source/notebooks/N2_Evaporation_File_Modification.rst index fc96c12..7f347df 100644 --- a/docs/source/notebooks/N2_Evaporation_File_Modification.rst +++ b/docs/source/notebooks/N2_Evaporation_File_Modification.rst @@ -25,10 +25,10 @@ dataset for the Gunnison. import pickle from string import Template import subprocess - + import matplotlib.pyplot as plt import numpy as np - import pandas as pd + import pandas as pd import statemodify as stm .. container:: alert alert-block alert-info @@ -44,16 +44,16 @@ dataset for the Gunnison. # statemod directory statemod_dir = "/usr/src/statemodify/statemod_gunnison_sjd" - + # root directory of statemod data for the target basin root_dir = os.path.join(statemod_dir, "src", "main", "fortran") - + # home directory of notebook instance home_dir = os.path.dirname(os.getcwd()) - + # path to the statemod executable statemod_exe = os.path.join(root_dir, "statemod-17.0.3-gfortran-lin-64bit-o3") - + # data directory and root name for the target basin data_dir = os.path.join( home_dir, @@ -62,13 +62,13 @@ dataset for the Gunnison. "gm2015_StateMod_modified", "StateMod" ) - + # directory to the target basin input files with root name for the basin basin_path = os.path.join(data_dir, "gm2015B") - + # scenarios output directory scenarios_dir = os.path.join(data_dir, "scenarios") - + # path to eva template file eva_template_file = os.path.join( home_dir, @@ -85,27 +85,27 @@ dataset for the Gunnison. .. parsed-literal:: - Startup log file for messages to this point: /home/jovyan/data/gm2015_StateMod_modified/gm2015_StateMod_modified/StateMod/gm2015B.rsp + Startup log file for messages to this point: /home/jovyan/data/gm2015_StateMod_modified/gm2015_StateMod_modified/StateMod/gm2015B.rsp Closing startup log file: statem.log - Opening dataset log file: /home/jovyan/data/gm2015_StateMod_modified/gm2015_StateMod_modified/StateMod/gm2015B.log + Opening dataset log file: /home/jovyan/data/gm2015_StateMod_modified/gm2015_StateMod_modified/StateMod/gm2015B.log ________________________________________________________________________ - - StateMod - State of Colorado - Water Supply Planning Model - - Version: 17.0.3 + + StateMod + State of Colorado - Water Supply Planning Model + + Version: 17.0.3 Last revision date: 2021/09/12 - + ________________________________________________________________________ - + Subroutine Execut Subroutine Datinp - + ... ________________________________________________________________________ Execut; Successful Termination - Statem; See detailed messages in dataset log file: /home/jovyan/data/gm2015_StateMod_modified/gm2015_StateMod_modified/StateMod/gm2015B.log + Statem; See detailed messages in dataset log file: /home/jovyan/data/gm2015_StateMod_modified/gm2015_StateMod_modified/StateMod/gm2015B.log Stop 0 @@ -122,21 +122,21 @@ for supplying emergency water to Lake Powell. .. code:: ipython3 - # create a directory to store the historical reservoir levels at Blue Mesa + # create a directory to store the historical reservoir levels at Blue Mesa output_dir = os.path.join(data_dir, "historic_reservoir_levels") - + if not os.path.exists(output_dir): os.makedirs(output_dir) - + # path the the xre file xre_file = os.path.join(data_dir, "gm2015B.xre") - + # structure ID for reservoir of interest - structure_ID = '6203532' - + structure_ID = '6203532' + # name of the reservoir - structure_name = 'Blue_Mesa' - + structure_name = 'Blue_Mesa' + # extract the target info into a Pandas data frame df = stm.extract_xre_data(structure_name=structure_name, structure_id=structure_ID, @@ -153,16 +153,16 @@ reservoir storage. .. code:: ipython3 output_xre_file = os.path.join(output_dir, "Blue_Mesa_xre_data.csv") - + # read output data into a data frame df = pd.read_csv( - output_xre_file, + output_xre_file, usecols=['Year','Init. Storage'], - index_col=False) - + index_col=False) + # calculate the annual average df = df.groupby('Year').mean().reset_index() - + df @@ -176,11 +176,11 @@ reservoir storage. .dataframe tbody tr th:only-of-type { vertical-align: middle; } - + .dataframe tbody tr th { vertical-align: top; } - + .dataframe thead th { text-align: right; } @@ -264,9 +264,9 @@ the 1930s dustbowl and 1950s drought and the severe early 2002 drought). .. code:: ipython3 fig, ax = plt.subplots() - + plt.plot(df['Year'], df['Init. Storage']) - + plt.title("Blue Mesa Storage") plt.xlabel("Year") plt.ylabel("Reservoir Storage (AF)") @@ -304,38 +304,38 @@ create 2 alternative states of the world and store them in the .. code:: ipython3 - # a dictionary to describe what you want to modify and the bounds for the Latin hypercube sample. + # a dictionary to describe what you want to modify and the bounds for the Latin hypercube sample. setup_dict = { "ids": ['10011'], "bounds": [-0.5, 1.0] } - + # create a directory to store the new files in if it does not exist output_directory = os.path.join(data_dir, "input_files") if not os.path.exists(output_directory): os.makedirs(output_directory) - + # scenario name scenario = "1" - + # the number of samples you wish to generate n_samples = 2 - + # seed value for reproducibility if so desired seed_value = 1 - + # number of rows to skip in file after comment skip_rows = 1 - + # name of field to query query_field = "id" - + # number of jobs to launch in parallel; -1 is all but 1 processor used n_jobs = -1 - + # basin to process basin_name = "Gunnison" - + # generate a batch of files using generated LHS stm.modify_eva(modify_dict=setup_dict, query_field=query_field, @@ -364,10 +364,10 @@ termed SOW 1 and SOW 2 respectively. # path to the numpy file containing the samples eva_samples_file = os.path.join(output_directory, "eva_2-samples_scenario-1.npy") - - # load samples + + # load samples sample_array = np.load(eva_samples_file) - + sample_array @@ -395,43 +395,43 @@ scenarios and extract the reservoir levels for Blue Mesa. # set realization and sample realization = 1 sample = np.arange(0, 2, 1) - + # read RSP template with open(eva_template_file) as template_obj: - + # read in file template_rsp = Template(template_obj.read()) - + for i in sample: - + # create scenario name scenario = f"S{i}_{realization}" - + # dictionary holding search keys and replacement values to update the template file d = {"EVA": f"../../input_files/gm2015B_{scenario}.eva"} - + # update the template new_rsp = template_rsp.safe_substitute(d) - + # construct simulated scenario directory simulated_scenario_dir = os.path.join(scenarios_dir, scenario) if not os.path.exists(simulated_scenario_dir): os.makedirs(simulated_scenario_dir) - + # target rsp file rsp_file = os.path.join(simulated_scenario_dir, f"gm2015B_{scenario}.rsp") - + # write updated rsp file with open(rsp_file, "w") as f1: f1.write(new_rsp) - + # construct simulated basin path simulated_basin_path = os.path.join(simulated_scenario_dir, f"gm2015B_{scenario}") - + # run StateMod print(f"Running: {scenario}") os.chdir(simulated_scenario_dir) - + subprocess.call([statemod_exe, simulated_basin_path, "-simulate"]) @@ -439,27 +439,27 @@ scenarios and extract the reservoir levels for Blue Mesa. .. parsed-literal:: Running: S0_1 - Startup log file for messages to this point: /home/jovyan/data/gm2015_StateMod_modified/gm2015_StateMod_modified/StateMod/scenarios/S0_1/gm2015B_S0_1.rsp + Startup log file for messages to this point: /home/jovyan/data/gm2015_StateMod_modified/gm2015_StateMod_modified/StateMod/scenarios/S0_1/gm2015B_S0_1.rsp Closing startup log file: statem.log - Opening dataset log file: /home/jovyan/data/gm2015_StateMod_modified/gm2015_StateMod_modified/StateMod/scenarios/S0_1/gm2015B_S0_1.log + Opening dataset log file: /home/jovyan/data/gm2015_StateMod_modified/gm2015_StateMod_modified/StateMod/scenarios/S0_1/gm2015B_S0_1.log ________________________________________________________________________ - - StateMod - State of Colorado - Water Supply Planning Model - - Version: 17.0.3 + + StateMod + State of Colorado - Water Supply Planning Model + + Version: 17.0.3 Last revision date: 2021/09/12 - + ________________________________________________________________________ - + Subroutine Execut Subroutine Datinp - + ... ________________________________________________________________________ Execut; Successful Termination - Statem; See detailed messages in dataset log file: /home/jovyan/data/gm2015_StateMod_modified/gm2015_StateMod_modified/StateMod/scenarios/S1_1/gm2015B_S1_1.log + Statem; See detailed messages in dataset log file: /home/jovyan/data/gm2015_StateMod_modified/gm2015_StateMod_modified/StateMod/scenarios/S1_1/gm2015B_S1_1.log Stop 0 @@ -473,16 +473,16 @@ at Blue Mesa in the two new SOWs. # SOW 1 output_dir= os.path.join(scenarios_dir, "S0_1") - + # path the the xre file - xre_file = os.path.join(output_dir, "gm2015B_S0_1.xre") - + xre_file = os.path.join(output_dir, "gm2015B_S0_1.xre") + # structure ID for reservoir of interest - structure_ID = '6203532' - + structure_ID = '6203532' + # name of the reservoir - structure_name = 'Blue_Mesa' - + structure_name = 'Blue_Mesa' + # extract the target info into a Pandas data frame df = stm.extract_xre_data(structure_name=structure_name, structure_id=structure_ID, @@ -492,14 +492,14 @@ at Blue Mesa in the two new SOWs. write_csv=True, write_parquet=None ) - + # SOW 2 output_dir= os.path.join(scenarios_dir, "S1_1") - + # path the the xre file - xre_file = os.path.join(output_dir, "gm2015B_S1_1.xre") - - + xre_file = os.path.join(output_dir, "gm2015B_S1_1.xre") + + # extract the target info into a Pandas data frame df = stm.extract_xre_data(structure_name=structure_name, structure_id=structure_ID, @@ -519,32 +519,32 @@ world and alternative states of the world. # historic reservoir directory historic_res_dir = os.path.join(data_dir, "historic_reservoir_levels") blue_mesa_file = os.path.join(historic_res_dir, "Blue_Mesa_xre_data.csv") - - # Import baseline dataframe - baseline = pd.read_csv(blue_mesa_file, index_col=False, usecols=['Year','Init. Storage']) + + # Import baseline dataframe + baseline = pd.read_csv(blue_mesa_file, index_col=False, usecols=['Year','Init. Storage']) baseline = baseline.groupby('Year').mean().reset_index() - + # Import SOW1 s0_1_file = os.path.join(scenarios_dir, "S0_1", "Blue_Mesa_xre_data.csv") - SOW1 = pd.read_csv(s0_1_file, index_col=False, usecols=['Year','Init. Storage']) + SOW1 = pd.read_csv(s0_1_file, index_col=False, usecols=['Year','Init. Storage']) SOW1 = SOW1.groupby('Year').mean().reset_index() - + # Import SOW2 s1_1_file = os.path.join(scenarios_dir, "S1_1", "Blue_Mesa_xre_data.csv") - SOW2 = pd.read_csv(s1_1_file, index_col=False, usecols=['Year','Init. Storage']) + SOW2 = pd.read_csv(s1_1_file, index_col=False, usecols=['Year','Init. Storage']) SOW2 = SOW2.groupby('Year').mean().reset_index() - - # Plot reservoir levels + + # Plot reservoir levels fig, ax = plt.subplots() - + plt.plot(baseline['Year'], baseline['Init. Storage'],label='Baseline') plt.plot(SOW1['Year'], SOW1['Init. Storage'],label='Reduced Evaporation') plt.plot(SOW2['Year'], SOW2['Init. Storage'],label='Increased Evaporation') - + plt.title("Blue Mesa Storage") plt.xlabel("Year") plt.ylabel("Reservoir Storage (AF)") - + plt.legend() @@ -589,4 +589,3 @@ e2020EF001503. :: 1. modify_eva() - diff --git a/docs/source/notebooks/N3_Reservoir_File_Modification.rst b/docs/source/notebooks/N3_Reservoir_File_Modification.rst index 9495edd..f61f26b 100644 --- a/docs/source/notebooks/N3_Reservoir_File_Modification.rst +++ b/docs/source/notebooks/N3_Reservoir_File_Modification.rst @@ -21,10 +21,10 @@ Step 1: Run a Historical Simulation in StateMod for the Uppper Colorado Subbasin import pickle from string import Template import subprocess - + import matplotlib.pyplot as plt import numpy as np - import pandas as pd + import pandas as pd import statemodify as stm .. container:: alert alert-block alert-info @@ -42,16 +42,16 @@ StateMod in a baseline simulation. # statemod directory statemod_dir = "/usr/src/statemodify/statemod_upper_co" - + # root directory of statemod data for the target basin root_dir = os.path.join(statemod_dir, "src", "main", "fortran") - + # home directory of notebook instance home_dir = os.path.dirname(os.getcwd()) - + # path to the statemod executable statemod_exe = os.path.join(root_dir, "statemod") - + # data directory and root name for the target basin data_dir = os.path.join( home_dir, @@ -59,17 +59,17 @@ StateMod in a baseline simulation. "cm2015_StateMod", "StateMod" ) - + # directory to the target basin input files with root name for the basin basin_path = os.path.join(data_dir, "cm2015B") - + # scenarios output directory scenarios_dir_res = os.path.join(data_dir, "scenarios_res") - + # parquet files output directory parquet_dir_res = os.path.join(data_dir, "parquet_res") - - + + # path to res template file res_template_file = os.path.join( home_dir, @@ -86,7 +86,7 @@ StateMod in a baseline simulation. .. code:: ipython3 - # Change directories first + # Change directories first os.chdir(data_dir) #This is needed specific to the Upper Colorado model as the path name is too long for the model to accept subprocess.call([statemod_exe, "cm2015B", "-simulate"]) @@ -94,28 +94,28 @@ StateMod in a baseline simulation. .. parsed-literal:: - Parse; Command line argument: - cm2015B -simulate + Parse; Command line argument: + cm2015B -simulate ________________________________________________________________________ - - StateMod - State of Colorado - Water Supply Planning Model - + + StateMod + State of Colorado - Water Supply Planning Model + Version: 15.00.01 Last revision date: 2015/10/28 - + ________________________________________________________________________ - - Opening log file cm2015B.log - + + Opening log file cm2015B.log + Subroutine Execut Subroutine Datinp - + ... ________________________________________________________________________ Execut; Successful Termination - Statem; See detailed messages in file: cm2015B.log + Statem; See detailed messages in file: cm2015B.log Stop 0 @@ -127,16 +127,16 @@ median water rights (47483.00000) and a smaller decree of 2.90 cfs. .. code:: ipython3 - #Extract shortages using statemodify convert_xdd() function - - # create a directory to store the historic shortages + #Extract shortages using statemodify convert_xdd() function + + # create a directory to store the historic shortages output_dir = os.path.join(data_dir, "historic_shortages") - + # create a directory to store the new files in if it does not exist output_directory = os.path.join(data_dir, "historic_shortages") if not os.path.exists(output_directory): os.makedirs(output_directory) - + stm.xdd.convert_xdd( # path to a directory where output .parquet files should be written output_path=output_dir, @@ -150,7 +150,7 @@ median water rights (47483.00000) and a smaller decree of 2.90 cfs. parallel_jobs=2, # convert to natural data types preserve_string_dtype=False - + ) @@ -164,13 +164,13 @@ Next we plot the shortages for Breckenridge. .. code:: ipython3 data=pd.read_parquet(output_dir +'/cm2015B.parquet',engine='pyarrow') - + fig, ax = plt.subplots() - + for name, group in data.groupby('structure_id'): ax.scatter( group['year'], group['shortage_total'], label=name) - + plt.xlabel("Year") plt.ylabel("Shortage (AF)") plt.title("Baseline Shortages for Breckenridge") @@ -220,16 +220,16 @@ storage at all reservoirs in the basin. scenario = "1" # basin name to process basin_name = "Upper_Colorado" - + # seed value for reproducibility if so desired seed_value = 1 - + # number of jobs to launch in parallel; -1 is all but 1 processor used n_jobs = 2 - + # number of samples to generate n_samples = 1 - + stm.modify_res(output_dir=output_directory, scenario=scenario, basin_name=basin_name, @@ -275,53 +275,53 @@ scenarios and extract the shortages for Breckenridge. # set realization and sample realization = 1 sample = np.arange(0, 2, 1) - + # read RSP template with open(res_template_file) as template_obj: - + # read in file template_rsp = Template(template_obj.read()) - + for i in sample: - + # create scenario name scenario = f"S{i}_{realization}" - + # dictionary holding search keys and replacement values to update the template file d = {"RES": f"../../input_files/cm2015B_{scenario}.res"} - + # update the template new_rsp = template_rsp.safe_substitute(d) - + # construct simulated scenario directory simulated_scenario_dir = os.path.join(scenarios_dir_res, scenario) if not os.path.exists(simulated_scenario_dir): os.makedirs(simulated_scenario_dir) - + # target rsp file rsp_file = os.path.join(simulated_scenario_dir, f"cm2015B_{scenario}.rsp") - + # write updated rsp file with open(rsp_file, "w") as f1: f1.write(new_rsp) - + # construct simulated basin path simulated_basin_path = f"cm2015B_{scenario}" - + # run StateMod print(f"Running: {scenario}") os.chdir(simulated_scenario_dir) - + subprocess.call([statemod_exe, simulated_basin_path, "-simulate"]) - - #Save output to parquet files + + #Save output to parquet files print('creating parquet for ' + scenario) - + output_directory = os.path.join(parquet_dir_res+"/scenario/"+ scenario) - + if not os.path.exists(output_directory): os.makedirs(output_directory) - + stm.xdd.convert_xdd( output_path=output_directory, allow_overwrite=True, @@ -335,28 +335,28 @@ scenarios and extract the shortages for Breckenridge. .. parsed-literal:: Running: S0_1 - Parse; Command line argument: - cm2015B_S0_1 -simulate + Parse; Command line argument: + cm2015B_S0_1 -simulate ________________________________________________________________________ - - StateMod - State of Colorado - Water Supply Planning Model - + + StateMod + State of Colorado - Water Supply Planning Model + Version: 15.00.01 Last revision date: 2015/10/28 - + ________________________________________________________________________ - - Opening log file cm2015B_S0_1.log - + + Opening log file cm2015B_S0_1.log + Subroutine Execut Subroutine Datinp - + ... ________________________________________________________________________ Execut; Successful Termination - Statem; See detailed messages in file: cm2015B_S0_1.log + Statem; See detailed messages in file: cm2015B_S0_1.log Stop 0 creating parquet for S0_1 @@ -369,32 +369,32 @@ scenarios and extract the shortages for Breckenridge. .. parsed-literal:: Running: S1_1 - Parse; Command line argument: - cm2015B_S1_1 -simulate + Parse; Command line argument: + cm2015B_S1_1 -simulate ________________________________________________________________________ - - StateMod - State of Colorado - Water Supply Planning Model - + + StateMod + State of Colorado - Water Supply Planning Model + Version: 15.00.01 Last revision date: 2015/10/28 - + ________________________________________________________________________ - - Opening log file cm2015B_S1_1.log - + + Opening log file cm2015B_S1_1.log + Subroutine Execut Subroutine Datinp - + ________________________________________________________________________ - Datinp; Control File (*.ctl) - + Datinp; Control File (*.ctl) + ________________________________________________________________________ Datinp; River Network File (*.rin) - + ________________________________________________________________________ Datinp; Reservoir Station File (*.res) - + ________________________________________________________________________ GetFile; Stopped in GetFile, see the log file (*.log) Stop 1 @@ -429,7 +429,7 @@ and alternative states of the world and plot the resulting shortages. 283922 0. 283923 220. 283924 201. - ... + ... 285280 0. 285281 0. 285282 0. @@ -443,17 +443,17 @@ and alternative states of the world and plot the resulting shortages. baseline=pd.read_parquet(data_dir+'/'+'historic_shortages/cm2015B.parquet',engine='pyarrow') SOW_1=pd.read_parquet(parquet_dir_res+'/scenario/S0_1/cm2015B_S0_1.parquet',engine='pyarrow') - + #Subtract shortages with respect to the baseline subset_df=pd.concat([baseline['year'],baseline['shortage_total'],SOW_1['shortage_total']],axis=1) subset_df = subset_df.set_axis(['Year', 'Baseline', 'SOW_1'], axis=1) subset_df['SOW_1_diff']=subset_df['SOW_1']-subset_df['Baseline'] - + #Plot shortages fig, ax = plt.subplots() - + ax.scatter(subset_df['Year'], subset_df['SOW_1_diff']) - + plt.xlabel("Year") plt.ylabel("Shortage (AF)") plt.title("Change in Breckenridge Shortages from the Baseline") @@ -516,4 +516,3 @@ e2020EF001503. :: 1. modify_res() - diff --git a/docs/source/notebooks/N4_Streamflow_File_Modification.rst b/docs/source/notebooks/N4_Streamflow_File_Modification.rst index bbee9ba..3f6ce84 100644 --- a/docs/source/notebooks/N4_Streamflow_File_Modification.rst +++ b/docs/source/notebooks/N4_Streamflow_File_Modification.rst @@ -57,10 +57,10 @@ Step 1: Fit Multi-Site HMM import pickle from string import Template import subprocess - + import matplotlib.pyplot as plt import numpy as np - import pandas as pd + import pandas as pd import statemodify as stm First we define directories and associated paths. @@ -69,16 +69,16 @@ First we define directories and associated paths. # statemod directory statemod_dir = "/usr/src/statemodify/statemod_upper_co" - + # root directory of statemod data for the target basin root_dir = os.path.join(statemod_dir, "src", "main", "fortran") - + # home directory of notebook instance home_dir = os.path.dirname(os.getcwd()) - + # path to the statemod executable statemod_exe = os.path.join(root_dir, "statemod") - + # data directory and root name for the target basin data_dir = os.path.join( home_dir, @@ -86,14 +86,14 @@ First we define directories and associated paths. "cm2015_StateMod", "StateMod" ) - + # directory to the target basin input files with root name for the basin basin_path = os.path.join(data_dir, "cm2015B") - + # scenarios output directory scenarios_dir = os.path.join(data_dir, "scenarios") - - # path to iwr/xbm file + + # path to iwr/xbm file xbm_iwr_template_file = os.path.join( home_dir, "data", @@ -108,21 +108,21 @@ store the HMM parameters. .. code:: ipython3 #Make directory to store HMM parameters - + output_dir = os.path.join(data_dir, "HMM_parameters") - + if not os.path.exists(output_dir): os.makedirs(output_dir) - + n_basins = 5 - + # choice to save parameters to NumPy array files save_parameters = True - + fit_array_dict = stm.hmm_multisite_fit(n_basins=n_basins, save_parameters=save_parameters, output_directory=output_dir) - + # unpack output dictionary unconditional_dry = fit_array_dict["unconditional_dry"] unconditional_wet = fit_array_dict["unconditional_wet"] @@ -145,12 +145,12 @@ in the ``HMM_Runs`` folder. .. code:: ipython3 #Create a folder to store the runs - + output_dir = os.path.join(data_dir, "HMM_Runs") - + if not os.path.exists(output_dir): os.makedirs(output_dir) - + # using the outputs of the fit function above; this function write output sample files to the output directory stm.hmm_multisite_sample(logAnnualQ_h, transition_matrix, @@ -209,30 +209,30 @@ function. .. code:: ipython3 - #Make directory to store input files - + #Make directory to store input files + output_dir = os.path.join(data_dir, "input_files") - + if not os.path.exists(output_dir): os.makedirs(output_dir) - - + + flow_realizations_directory = os.path.join(data_dir, "HMM_Runs") - + scenario = "1" - + # basin name to process basin_name = "Upper_Colorado" - + # seed value for reproducibility if so desired seed_value = 123 - + # number of jobs to launch in parallel; -1 is all but 1 processor used n_jobs = 2 - + # number of samples to generate (how many new xbm and iwr files); produces an IWR multiplier n_samples = 1 - + # generate a batch of files using generated LHS stm.modify_xbm_iwr(output_dir=output_dir, flow_realizations_directory=flow_realizations_directory, @@ -242,7 +242,7 @@ function. n_jobs=n_jobs, n_samples=n_samples, save_sample=True, - randomly_select_flow_sample=True) + randomly_select_flow_sample=True) Step 4: Read in the New Input Files and Run StateMod : Streamflow Example ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -268,43 +268,43 @@ levels. # set realization and sample realization = 1 sample = np.arange(0, 1, 1) - + # read RSP template with open(xbm_iwr_template_file) as template_obj: - + # read in file template_rsp = Template(template_obj.read()) - + for i in sample: - + # create scenario name scenario = f"S{i}_{realization}" - + # dictionary holding search keys and replacement values to update the template file d = {"XBM": f"../../input_files/cm2015B_{scenario}.xbm","IWR": f"../../input_files/cm2015B_{scenario}.iwr"} - + # update the template new_rsp = template_rsp.safe_substitute(d) - + # construct simulated scenario directory simulated_scenario_dir = os.path.join(scenarios_dir, scenario) if not os.path.exists(simulated_scenario_dir): os.makedirs(simulated_scenario_dir) - + # target rsp file rsp_file = os.path.join(simulated_scenario_dir, f"cm2015B_{scenario}.rsp") - + # write updated rsp file with open(rsp_file, "w") as f1: f1.write(new_rsp) - + # construct simulated basin path simulated_basin_path = f"cm2015B_{scenario}" - + # run StateMod print(f"Running: {scenario}") os.chdir(simulated_scenario_dir) - + subprocess.call([statemod_exe, simulated_basin_path, "-simulate"]) @@ -312,28 +312,28 @@ levels. .. parsed-literal:: Running: S0_1 - Parse; Command line argument: - cm2015B_S0_1 -simulate + Parse; Command line argument: + cm2015B_S0_1 -simulate ________________________________________________________________________ - - StateMod - State of Colorado - Water Supply Planning Model - + + StateMod + State of Colorado - Water Supply Planning Model + Version: 15.00.01 Last revision date: 2015/10/28 - + ________________________________________________________________________ - - Opening log file cm2015B_S0_1.log - + + Opening log file cm2015B_S0_1.log + Subroutine Execut Subroutine Datinp - + ... ________________________________________________________________________ Execut; Successful Termination - Statem; See detailed messages in file: cm2015B_S0_1.log + Statem; See detailed messages in file: cm2015B_S0_1.log Stop 0 @@ -351,10 +351,10 @@ the historical 105-year period. # Example with Granby Lake zip_file_path = os.path.join(home_dir, 'data', 'Granby_Dataset.zip') final_directory = os.path.join(home_dir, 'data/') - - !unzip $zip_file_path -d $final_directory + + !unzip $zip_file_path -d $final_directory granby_hmm, granby_hist, granby_hist_mean, granby_hist_1p = stm.read_xre(os.path.join(home_dir,"data/Upper_Colorado/"), 'Granby') - + # Plot quantiles stm.plot_res_quantiles(granby_hmm, granby_hist_mean, 'Lake Granby') @@ -415,13 +415,13 @@ assessments. :: - 2. Fitting Hidden Markov Models: Sample Scripts + 2. Fitting Hidden Markov Models: Sample Scripts .. container:: :: - 3. A Hidden-Markov Modeling Approach to Creating Synthetic Streamflow Scenarios Tutorial + 3. A Hidden-Markov Modeling Approach to Creating Synthetic Streamflow Scenarios Tutorial Notebook Specific References ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -447,4 +447,3 @@ e2020EF001503. :: 1. modify_xbm_iwr() - diff --git a/docs/source/notebooks/N5_Batch_Modification.rst b/docs/source/notebooks/N5_Batch_Modification.rst index 7389d78..46465e1 100644 --- a/docs/source/notebooks/N5_Batch_Modification.rst +++ b/docs/source/notebooks/N5_Batch_Modification.rst @@ -21,10 +21,10 @@ using the ``modify_batch()`` function. import pickle from string import Template import subprocess - + import matplotlib.pyplot as plt import numpy as np - import pandas as pd + import pandas as pd import statemodify as stm .. container:: alert alert-block alert-info @@ -39,16 +39,16 @@ using the ``modify_batch()`` function. # statemod directory statemod_dir = "/usr/src/statemodify/statemod_upper_co" - + # root directory of statemod data for the target basin root_dir = os.path.join(statemod_dir, "src", "main", "fortran") - + # home directory of notebook instance home_dir = os.path.dirname(os.getcwd()) - + # path to the statemod executable statemod_exe = os.path.join(root_dir, "statemod") - + # data directory and root name for the target basin data_dir = os.path.join( home_dir, @@ -56,17 +56,17 @@ using the ``modify_batch()`` function. "cm2015_StateMod", "StateMod" ) - + # directory to the target basin input files with root name for the basin basin_path = os.path.join(data_dir, "cm2015B") - + # scenarios output directory scenarios_dir = os.path.join(data_dir, "scenarios") - + #parquet output directory parquet_dir=os.path.join(data_dir, "parquet") - - + + # path to template file multi_template_file = os.path.join( home_dir, @@ -99,13 +99,13 @@ additives are applied to the target IDs listed. .. code:: ipython3 import statemodify as stm - + # variables that apply to multiple functions output_dir = os.path.join(data_dir, "input_files") basin_name = "Upper_Colorado" scenario = "1" seed_value = 77 - + # problem dictionary problem_dict = { "n_samples": 1, @@ -144,7 +144,7 @@ additives are applied to the target IDs listed. "ids": ["3600507", "3600603"] } } - + # run in batch fn_parameter_dict = stm.modify_batch(problem_dict=problem_dict) @@ -178,53 +178,53 @@ a specific user (ID: 3601008). # set realization and sample realization = 1 sample = np.arange(0, 1, 1) - + # read RSP template with open(multi_template_file) as template_obj: - + # read in file template_rsp = Template(template_obj.read()) - + for i in sample: - + # create scenario name scenario = f"S{i}_{realization}" - + # dictionary holding search keys and replacement values to update the template file d = {"EVA": f"../../input_files/cm2015B_{scenario}.eva","DDM": f"../../input_files/cm2015B_{scenario}.ddm","DDR": f"../../input_files/cm2015B_{scenario}.ddr"} - + # update the template new_rsp = template_rsp.safe_substitute(d) - + # construct simulated scenario directory simulated_scenario_dir = os.path.join(scenarios_dir, scenario) if not os.path.exists(simulated_scenario_dir): os.makedirs(simulated_scenario_dir) - + # target rsp file rsp_file = os.path.join(simulated_scenario_dir, f"cm2015B_{scenario}.rsp") - + # write updated rsp file with open(rsp_file, "w") as f1: f1.write(new_rsp) - + # construct simulated basin path simulated_basin_path = f"cm2015B_{scenario}" - + # run StateMod print(f"Running: {scenario}") os.chdir(simulated_scenario_dir) - + subprocess.call([statemod_exe, simulated_basin_path, "-simulate"]) - - #Save output to parquet files + + #Save output to parquet files print('creating parquet for ' + scenario) - + output_directory = os.path.join(parquet_dir+ "/scenario/"+ scenario) - + if not os.path.exists(output_directory): os.makedirs(output_directory) - + stm.xdd.convert_xdd(output_path=output_directory,allow_overwrite=False,xdd_files=scenarios_dir + "/"+ scenario + "/cm2015B_"+scenario+".xdd",id_subset=['3601008'],parallel_jobs=2) @@ -232,23 +232,23 @@ a specific user (ID: 3601008). .. parsed-literal:: Running: S0_1 - Parse; Command line argument: - cm2015B_S0_1 -simulate + Parse; Command line argument: + cm2015B_S0_1 -simulate ________________________________________________________________________ - - StateMod - State of Colorado - Water Supply Planning Model - + + StateMod + State of Colorado - Water Supply Planning Model + Version: 15.00.01 Last revision date: 2015/10/28 - + ________________________________________________________________________ - - Opening log file cm2015B_S0_1.log - + + Opening log file cm2015B_S0_1.log + Subroutine Execut Subroutine Datinp - + ... @@ -270,4 +270,3 @@ the prior notebooks. :: 1. modify_batch() - diff --git a/notebooks/N1_Demand_WaterRights_File_Modification.ipynb b/notebooks/N1_Demand_WaterRights_File_Modification.ipynb index ebed839..3cc40e0 100644 --- a/notebooks/N1_Demand_WaterRights_File_Modification.ipynb +++ b/notebooks/N1_Demand_WaterRights_File_Modification.ipynb @@ -77,12 +77,13 @@ "import logging\n", "import os\n", "import pickle\n", - "from string import Template\n", "import subprocess\n", + "from string import Template\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", - "import pandas as pd \n", + "import pandas as pd\n", + "\n", "import statemodify as stm" ] }, @@ -122,11 +123,7 @@ "\n", "# data directory and root name for the target basin\n", "data_dir = os.path.join(\n", - " home_dir,\n", - " \"data\",\n", - " \"sj2015_StateMod_modified\",\n", - " \"sj2015_StateMod_modified\",\n", - " \"StateMod\"\n", + " home_dir, \"data\", \"sj2015_StateMod_modified\", \"sj2015_StateMod_modified\", \"StateMod\"\n", ")\n", "\n", "# directory to the target basin input files with root name for the basin\n", @@ -141,17 +138,9 @@ "parquet_dir_ddr = os.path.join(data_dir, \"parquet_ddr\")\n", "\n", "# path to ddm and ddr template file\n", - "ddm_template_file = os.path.join(\n", - " home_dir,\n", - " \"data\",\n", - " \"sj2015B_template_ddm.rsp\"\n", - ")\n", + "ddm_template_file = os.path.join(home_dir, \"data\", \"sj2015B_template_ddm.rsp\")\n", "\n", - "ddr_template_file = os.path.join(\n", - " home_dir,\n", - " \"data\",\n", - " \"sj2015B_template_ddr.rsp\"\n", - ")" + "ddr_template_file = os.path.join(home_dir, \"data\", \"sj2015B_template_ddr.rsp\")" ] }, { @@ -1596,7 +1585,7 @@ ], "source": [ "# run statemod\n", - "subprocess.call([statemod_exe, basin_path, \"-simulate\"])\n" + "subprocess.call([statemod_exe, basin_path, \"-simulate\"])" ] }, { @@ -1624,9 +1613,9 @@ } ], "source": [ - "#Extract shortages using statemodify convert_xdd() function \n", + "# Extract shortages using statemodify convert_xdd() function\n", "\n", - "# create a directory to store the historical shortages \n", + "# create a directory to store the historical shortages\n", "output_dir = os.path.join(data_dir, \"historic_shortages\")\n", "\n", "# create a directory to store the new files in if it does not exist\n", @@ -1642,12 +1631,12 @@ " # path, glob, or a list of paths/globs to the .xdd files you want to convert\n", " xdd_files=os.path.join(data_dir, \"*.xdd\"),\n", " # if the output .parquet files should only contain a subset of structure ids, list them here; None for all\n", - " id_subset=['2900501','2900519','2900555'],\n", + " id_subset=[\"2900501\", \"2900519\", \"2900555\"],\n", " # how many .xdd files to convert in parallel; optimally you will want 2-4 CPUs per parallel process\n", " parallel_jobs=4,\n", " # convert to natural data types\n", - " preserve_string_dtype=False\n", - ")\n" + " preserve_string_dtype=False,\n", + ")" ] }, { @@ -2086,7 +2075,7 @@ } ], "source": [ - "data=pd.read_parquet(os.path.join(output_dir,'sj2015B.parquet'),engine='pyarrow')\n", + "data = pd.read_parquet(os.path.join(output_dir, \"sj2015B.parquet\"), engine=\"pyarrow\")\n", "data" ] }, @@ -2140,9 +2129,8 @@ "source": [ "fig, ax = plt.subplots()\n", "\n", - "for name, group in data.groupby('structure_id'):\n", - " ax.scatter(\n", - " group['year'], group['shortage_total'], label=name)\n", + "for name, group in data.groupby(\"structure_id\"):\n", + " ax.scatter(group[\"year\"], group[\"shortage_total\"], label=name)\n", "\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"Shortage (AF)\")\n", @@ -2187,10 +2175,7 @@ "outputs": [], "source": [ "# a dictionary to describe what users you want to modify and the bounds for the LHS\n", - "setup_dict = {\n", - " \"ids\": [\"2900501\", \"2900519\",\"2900555\"],\n", - " \"bounds\": [0.5, 1.5]\n", - "}\n", + "setup_dict = {\"ids\": [\"2900501\", \"2900519\", \"2900555\"], \"bounds\": [0.5, 1.5]}\n", "\n", "output_directory = output_dir = os.path.join(data_dir, \"input_files\")\n", "\n", @@ -2231,8 +2216,8 @@ " data_specification_file=None,\n", " min_bound_value=-0.5,\n", " max_bound_value=1.5,\n", - " save_sample=True\n", - ")\n" + " save_sample=True,\n", + ")" ] }, { @@ -2264,7 +2249,7 @@ } ], "source": [ - "sample_array = np.load(output_directory+'/ddm_2-samples_scenario-1.npy')\n", + "sample_array = np.load(output_directory + \"/ddm_2-samples_scenario-1.npy\")\n", "sample_array" ] }, @@ -3740,58 +3725,63 @@ "\n", "# read RSP template\n", "with open(ddm_template_file) as template_obj:\n", - " \n", " # read in file\n", " template_rsp = Template(template_obj.read())\n", "\n", " for i in sample:\n", - " \n", " # create scenario name\n", " scenario = f\"S{i}_{realization}\"\n", - " \n", + "\n", " # dictionary holding search keys and replacement values to update the template file\n", " d = {\"DDM\": f\"../../input_files/sj2015B_{scenario}.ddm\"}\n", - " \n", + "\n", " # update the template\n", " new_rsp = template_rsp.safe_substitute(d)\n", - " \n", + "\n", " # construct simulated scenario directory\n", " simulated_scenario_dir = os.path.join(scenarios_dir_ddm, scenario)\n", " if not os.path.exists(simulated_scenario_dir):\n", " os.makedirs(simulated_scenario_dir)\n", - " \n", + "\n", " # target rsp file\n", " rsp_file = os.path.join(simulated_scenario_dir, f\"sj2015B_{scenario}.rsp\")\n", - " \n", + "\n", " # write updated rsp file\n", " with open(rsp_file, \"w\") as f1:\n", " f1.write(new_rsp)\n", - " \n", + "\n", " # construct simulated basin path\n", - " simulated_basin_path = os.path.join(simulated_scenario_dir, f\"sj2015B_{scenario}\")\n", + " simulated_basin_path = os.path.join(\n", + " simulated_scenario_dir, f\"sj2015B_{scenario}\"\n", + " )\n", "\n", " # run StateMod\n", " print(f\"Running: {scenario}\")\n", " os.chdir(simulated_scenario_dir)\n", "\n", " subprocess.call([statemod_exe, simulated_basin_path, \"-simulate\"])\n", - " \n", - " #Save output to parquet files \n", - " print('creating parquet for ' + scenario)\n", - " \n", - " output_directory = os.path.join(parquet_dir_ddm+\"/scenario/\"+ scenario)\n", - " \n", + "\n", + " # Save output to parquet files\n", + " print(\"creating parquet for \" + scenario)\n", + "\n", + " output_directory = os.path.join(parquet_dir_ddm + \"/scenario/\" + scenario)\n", + "\n", " if not os.path.exists(output_directory):\n", " os.makedirs(output_directory)\n", - " \n", + "\n", " stm.xdd.convert_xdd(\n", " output_path=output_directory,\n", " allow_overwrite=False,\n", - " xdd_files=scenarios_dir_ddm + \"/\"+ scenario + \"/sj2015B_\"+scenario+\".xdd\",\n", - " id_subset=['2900501','2900519','2900555'],\n", + " xdd_files=scenarios_dir_ddm\n", + " + \"/\"\n", + " + scenario\n", + " + \"/sj2015B_\"\n", + " + scenario\n", + " + \".xdd\",\n", + " id_subset=[\"2900501\", \"2900519\", \"2900555\"],\n", " parallel_jobs=4,\n", - " preserve_string_dtype=False\n", - " )\n" + " preserve_string_dtype=False,\n", + " )" ] }, { @@ -3836,7 +3826,7 @@ ], "source": [ "output_directory = os.path.join(data_dir, \"input_files\")\n", - "sample_array = np.load(output_directory+'/ddm_2-samples_scenario-1.npy')\n", + "sample_array = np.load(output_directory + \"/ddm_2-samples_scenario-1.npy\")\n", "sample_array" ] }, @@ -3879,21 +3869,35 @@ ], "source": [ "# Read in raw parquet files\n", - "baseline=pd.read_parquet(data_dir+'/historic_shortages/sj2015B.parquet',engine='pyarrow')\n", - "SOW_1=pd.read_parquet(parquet_dir_ddm+'/scenario/S0_1/sj2015B_S0_1.parquet',engine='pyarrow')\n", - "SOW_2=pd.read_parquet(parquet_dir_ddm+'/scenario/S1_1/sj2015B_S1_1.parquet',engine='pyarrow')\n", + "baseline = pd.read_parquet(\n", + " data_dir + \"/historic_shortages/sj2015B.parquet\", engine=\"pyarrow\"\n", + ")\n", + "SOW_1 = pd.read_parquet(\n", + " parquet_dir_ddm + \"/scenario/S0_1/sj2015B_S0_1.parquet\", engine=\"pyarrow\"\n", + ")\n", + "SOW_2 = pd.read_parquet(\n", + " parquet_dir_ddm + \"/scenario/S1_1/sj2015B_S1_1.parquet\", engine=\"pyarrow\"\n", + ")\n", "\n", "# Subtract shortages with respect to the baseline\n", - "subset_df=pd.concat([baseline['year'],baseline['shortage_total'],SOW_1['shortage_total'],SOW_2['shortage_total']],axis=1)\n", - "subset_df = subset_df.set_axis(['Year', 'Baseline', 'SOW_1','SOW_2'], axis=1)\n", - "subset_df['SOW_1_diff'] = subset_df['SOW_1']-subset_df['Baseline']\n", - "subset_df['SOW_2_diff'] = subset_df['SOW_2']-subset_df['Baseline']\n", + "subset_df = pd.concat(\n", + " [\n", + " baseline[\"year\"],\n", + " baseline[\"shortage_total\"],\n", + " SOW_1[\"shortage_total\"],\n", + " SOW_2[\"shortage_total\"],\n", + " ],\n", + " axis=1,\n", + ")\n", + "subset_df = subset_df.set_axis([\"Year\", \"Baseline\", \"SOW_1\", \"SOW_2\"], axis=1)\n", + "subset_df[\"SOW_1_diff\"] = subset_df[\"SOW_1\"] - subset_df[\"Baseline\"]\n", + "subset_df[\"SOW_2_diff\"] = subset_df[\"SOW_2\"] - subset_df[\"Baseline\"]\n", "\n", "# Plot shortages\n", "fig, ax = plt.subplots()\n", "\n", - "ax.scatter(subset_df['Year'], subset_df['SOW_1_diff'],label='Decreased Demand')\n", - "ax.scatter(subset_df['Year'], subset_df['SOW_2_diff'],label='Increased Demand')\n", + "ax.scatter(subset_df[\"Year\"], subset_df[\"SOW_1_diff\"], label=\"Decreased Demand\")\n", + "ax.scatter(subset_df[\"Year\"], subset_df[\"SOW_2_diff\"], label=\"Increased Demand\")\n", "\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"Shortage (AF)\")\n", @@ -3938,11 +3942,9 @@ "setup_dict = {\n", " # ids can either be 'struct' or 'id' values\n", " \"ids\": [\"2900501\"],\n", - "\n", " # turn id on or off completely or for a given period\n", " # if 0 = off, 1 = on, YYYY = on for years >= YYYY, -YYYY = off for years > YYYY; see file header\n", " \"on_off\": [1],\n", - "\n", " # apply rank of administrative order where 0 is lowest (senior) and n is highest (junior); None is no change\n", " \"admin\": [1],\n", "}\n", @@ -3985,8 +3987,8 @@ " data_specification_file=None,\n", " min_bound_value=-0.5,\n", " max_bound_value=1.5,\n", - " save_sample=True\n", - ")\n" + " save_sample=True,\n", + ")" ] }, { @@ -5449,58 +5451,63 @@ "\n", "# read RSP template\n", "with open(ddr_template_file) as template_obj:\n", - " \n", " # read in file\n", " template_rsp = Template(template_obj.read())\n", "\n", " for i in sample:\n", - " \n", " # create scenario name\n", " scenario = f\"S{i}_{realization}\"\n", - " \n", + "\n", " # dictionary holding search keys and replacement values to update the template file\n", " d = {\"DDR\": f\"../../input_files/sj2015B_{scenario}.ddr\"}\n", - " \n", + "\n", " # update the template\n", " new_rsp = template_rsp.safe_substitute(d)\n", - " \n", + "\n", " # construct simulated scenario directory\n", " simulated_scenario_dir = os.path.join(scenarios_dir_ddr, scenario)\n", " if not os.path.exists(simulated_scenario_dir):\n", " os.makedirs(simulated_scenario_dir)\n", - " \n", + "\n", " # target rsp file\n", " rsp_file = os.path.join(simulated_scenario_dir, f\"sj2015B_{scenario}.rsp\")\n", - " \n", + "\n", " # write updated rsp file\n", " with open(rsp_file, \"w\") as f1:\n", " f1.write(new_rsp)\n", - " \n", + "\n", " # construct simulated basin path\n", - " simulated_basin_path = os.path.join(simulated_scenario_dir, f\"sj2015B_{scenario}\")\n", + " simulated_basin_path = os.path.join(\n", + " simulated_scenario_dir, f\"sj2015B_{scenario}\"\n", + " )\n", "\n", " # run StateMod\n", " print(f\"Running: {scenario}\")\n", " os.chdir(simulated_scenario_dir)\n", "\n", " subprocess.call([statemod_exe, simulated_basin_path, \"-simulate\"])\n", - " \n", - " #Save output to parquet files \n", - " print('creating parquet for ' + scenario)\n", - " \n", - " output_directory = os.path.join(parquet_dir_ddr+\"/scenario/\"+ scenario)\n", - " \n", + "\n", + " # Save output to parquet files\n", + " print(\"creating parquet for \" + scenario)\n", + "\n", + " output_directory = os.path.join(parquet_dir_ddr + \"/scenario/\" + scenario)\n", + "\n", " if not os.path.exists(output_directory):\n", " os.makedirs(output_directory)\n", - " \n", + "\n", " stm.xdd.convert_xdd(\n", " output_path=output_directory,\n", " allow_overwrite=False,\n", - " xdd_files=scenarios_dir_ddr + \"/\"+ scenario + \"/sj2015B_\"+scenario+\".xdd\",\n", - " id_subset=['2900501'],\n", + " xdd_files=scenarios_dir_ddr\n", + " + \"/\"\n", + " + scenario\n", + " + \"/sj2015B_\"\n", + " + scenario\n", + " + \".xdd\",\n", + " id_subset=[\"2900501\"],\n", " parallel_jobs=2,\n", - " preserve_string_dtype=False\n", - " )\n" + " preserve_string_dtype=False,\n", + " )" ] }, { @@ -5542,18 +5549,24 @@ ], "source": [ "# Read in raw parquet files\n", - "baseline=pd.read_parquet(data_dir+'/historic_shortages/sj2015B.parquet',engine='pyarrow')\n", - "SOW_1=pd.read_parquet(parquet_dir_ddr+ '/scenario/S0_1/sj2015B_S0_1.parquet',engine='pyarrow')\n", + "baseline = pd.read_parquet(\n", + " data_dir + \"/historic_shortages/sj2015B.parquet\", engine=\"pyarrow\"\n", + ")\n", + "SOW_1 = pd.read_parquet(\n", + " parquet_dir_ddr + \"/scenario/S0_1/sj2015B_S0_1.parquet\", engine=\"pyarrow\"\n", + ")\n", "\n", "# Subtract shortages with respect to the baseline\n", - "subset_df=pd.concat([baseline['year'],baseline['shortage_total'],SOW_1['shortage_total']],axis=1)\n", - "subset_df = subset_df.set_axis(['Year', 'Baseline', 'SOW_1'], axis=1)\n", - "subset_df['diff']=subset_df['SOW_1']-subset_df['Baseline']\n", + "subset_df = pd.concat(\n", + " [baseline[\"year\"], baseline[\"shortage_total\"], SOW_1[\"shortage_total\"]], axis=1\n", + ")\n", + "subset_df = subset_df.set_axis([\"Year\", \"Baseline\", \"SOW_1\"], axis=1)\n", + "subset_df[\"diff\"] = subset_df[\"SOW_1\"] - subset_df[\"Baseline\"]\n", "\n", "# Plot shortages\n", "fig, ax = plt.subplots()\n", "\n", - "ax.scatter(subset_df['Year'], subset_df['diff'])\n", + "ax.scatter(subset_df[\"Year\"], subset_df[\"diff\"])\n", "\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"Shortage (AF)\")\n", diff --git a/notebooks/N2_Evaporation_File_Modification.ipynb b/notebooks/N2_Evaporation_File_Modification.ipynb index 52d2e91..a951ea3 100644 --- a/notebooks/N2_Evaporation_File_Modification.ipynb +++ b/notebooks/N2_Evaporation_File_Modification.ipynb @@ -47,12 +47,13 @@ "import logging\n", "import os\n", "import pickle\n", - "from string import Template\n", "import subprocess\n", + "from string import Template\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", - "import pandas as pd \n", + "import pandas as pd\n", + "\n", "import statemodify as stm" ] }, @@ -96,11 +97,7 @@ "\n", "# data directory and root name for the target basin\n", "data_dir = os.path.join(\n", - " home_dir,\n", - " \"data\",\n", - " \"gm2015_StateMod_modified\",\n", - " \"gm2015_StateMod_modified\",\n", - " \"StateMod\"\n", + " home_dir, \"data\", \"gm2015_StateMod_modified\", \"gm2015_StateMod_modified\", \"StateMod\"\n", ")\n", "\n", "# directory to the target basin input files with root name for the basin\n", @@ -110,11 +107,7 @@ "scenarios_dir = os.path.join(data_dir, \"scenarios\")\n", "\n", "# path to eva template file\n", - "eva_template_file = os.path.join(\n", - " home_dir,\n", - " \"data\",\n", - " \"gm2015B_template_eva.rsp\"\n", - ")" + "eva_template_file = os.path.join(home_dir, \"data\", \"gm2015B_template_eva.rsp\")" ] }, { @@ -1570,7 +1563,7 @@ ], "source": [ "# run statemod\n", - "subprocess.call([statemod_exe, basin_path, \"-simulate\"])\n" + "subprocess.call([statemod_exe, basin_path, \"-simulate\"])" ] }, { @@ -1590,7 +1583,7 @@ }, "outputs": [], "source": [ - "# create a directory to store the historical reservoir levels at Blue Mesa \n", + "# create a directory to store the historical reservoir levels at Blue Mesa\n", "output_dir = os.path.join(data_dir, \"historic_reservoir_levels\")\n", "\n", "if not os.path.exists(output_dir):\n", @@ -1600,19 +1593,20 @@ "xre_file = os.path.join(data_dir, \"gm2015B.xre\")\n", "\n", "# structure ID for reservoir of interest\n", - "structure_ID = '6203532' \n", + "structure_ID = \"6203532\"\n", "\n", "# name of the reservoir\n", - "structure_name = 'Blue_Mesa' \n", + "structure_name = \"Blue_Mesa\"\n", "\n", "# extract the target info into a Pandas data frame\n", - "df = stm.extract_xre_data(structure_name=structure_name,\n", - " structure_id=structure_ID,\n", - " input_file=xre_file,\n", - " basin_name=None,\n", - " output_directory=output_dir,\n", - " write_csv=True,\n", - " write_parquet=None\n", + "df = stm.extract_xre_data(\n", + " structure_name=structure_name,\n", + " structure_id=structure_ID,\n", + " input_file=xre_file,\n", + " basin_name=None,\n", + " output_directory=output_dir,\n", + " write_csv=True,\n", + " write_parquet=None,\n", ")" ] }, @@ -1744,15 +1738,12 @@ "output_xre_file = os.path.join(output_dir, \"Blue_Mesa_xre_data.csv\")\n", "\n", "# read output data into a data frame\n", - "df = pd.read_csv(\n", - " output_xre_file, \n", - " usecols=['Year','Init. Storage'],\n", - " index_col=False) \n", + "df = pd.read_csv(output_xre_file, usecols=[\"Year\", \"Init. Storage\"], index_col=False)\n", "\n", "# calculate the annual average\n", - "df = df.groupby('Year').mean().reset_index()\n", + "df = df.groupby(\"Year\").mean().reset_index()\n", "\n", - "df\n" + "df" ] }, { @@ -1795,11 +1786,11 @@ "source": [ "fig, ax = plt.subplots()\n", "\n", - "plt.plot(df['Year'], df['Init. Storage'])\n", + "plt.plot(df[\"Year\"], df[\"Init. Storage\"])\n", "\n", "plt.title(\"Blue Mesa Storage\")\n", "plt.xlabel(\"Year\")\n", - "plt.ylabel(\"Reservoir Storage (AF)\")\n" + "plt.ylabel(\"Reservoir Storage (AF)\")" ] }, { @@ -1827,11 +1818,8 @@ }, "outputs": [], "source": [ - "# a dictionary to describe what you want to modify and the bounds for the Latin hypercube sample. \n", - "setup_dict = {\n", - " \"ids\": ['10011'],\n", - " \"bounds\": [-0.5, 1.0]\n", - "}\n", + "# a dictionary to describe what you want to modify and the bounds for the Latin hypercube sample.\n", + "setup_dict = {\"ids\": [\"10011\"], \"bounds\": [-0.5, 1.0]}\n", "\n", "# create a directory to store the new files in if it does not exist\n", "output_directory = os.path.join(data_dir, \"input_files\")\n", @@ -1860,22 +1848,24 @@ "basin_name = \"Gunnison\"\n", "\n", "# generate a batch of files using generated LHS\n", - "stm.modify_eva(modify_dict=setup_dict,\n", - " query_field=query_field,\n", - " output_dir=output_directory,\n", - " scenario=scenario,\n", - " basin_name=basin_name,\n", - " sampling_method=\"LHS\",\n", - " n_samples=n_samples,\n", - " skip_rows=skip_rows,\n", - " n_jobs=n_jobs,\n", - " seed_value=seed_value,\n", - " template_file=None,\n", - " factor_method=\"add\",\n", - " data_specification_file=None,\n", - " min_bound_value=-0.5,\n", - " max_bound_value=1.0,\n", - " save_sample=True)\n" + "stm.modify_eva(\n", + " modify_dict=setup_dict,\n", + " query_field=query_field,\n", + " output_dir=output_directory,\n", + " scenario=scenario,\n", + " basin_name=basin_name,\n", + " sampling_method=\"LHS\",\n", + " n_samples=n_samples,\n", + " skip_rows=skip_rows,\n", + " n_jobs=n_jobs,\n", + " seed_value=seed_value,\n", + " template_file=None,\n", + " factor_method=\"add\",\n", + " data_specification_file=None,\n", + " min_bound_value=-0.5,\n", + " max_bound_value=1.0,\n", + " save_sample=True,\n", + ")" ] }, { @@ -1910,10 +1900,10 @@ "# path to the numpy file containing the samples\n", "eva_samples_file = os.path.join(output_directory, \"eva_2-samples_scenario-1.npy\")\n", "\n", - "# load samples \n", + "# load samples\n", "sample_array = np.load(eva_samples_file)\n", "\n", - "sample_array\n" + "sample_array" ] }, { @@ -4807,41 +4797,41 @@ "\n", "# read RSP template\n", "with open(eva_template_file) as template_obj:\n", - " \n", " # read in file\n", " template_rsp = Template(template_obj.read())\n", "\n", " for i in sample:\n", - " \n", " # create scenario name\n", " scenario = f\"S{i}_{realization}\"\n", - " \n", + "\n", " # dictionary holding search keys and replacement values to update the template file\n", " d = {\"EVA\": f\"../../input_files/gm2015B_{scenario}.eva\"}\n", - " \n", + "\n", " # update the template\n", " new_rsp = template_rsp.safe_substitute(d)\n", - " \n", + "\n", " # construct simulated scenario directory\n", " simulated_scenario_dir = os.path.join(scenarios_dir, scenario)\n", " if not os.path.exists(simulated_scenario_dir):\n", " os.makedirs(simulated_scenario_dir)\n", - " \n", + "\n", " # target rsp file\n", " rsp_file = os.path.join(simulated_scenario_dir, f\"gm2015B_{scenario}.rsp\")\n", - " \n", + "\n", " # write updated rsp file\n", " with open(rsp_file, \"w\") as f1:\n", " f1.write(new_rsp)\n", - " \n", + "\n", " # construct simulated basin path\n", - " simulated_basin_path = os.path.join(simulated_scenario_dir, f\"gm2015B_{scenario}\")\n", + " simulated_basin_path = os.path.join(\n", + " simulated_scenario_dir, f\"gm2015B_{scenario}\"\n", + " )\n", "\n", " # run StateMod\n", " print(f\"Running: {scenario}\")\n", " os.chdir(simulated_scenario_dir)\n", "\n", - " subprocess.call([statemod_exe, simulated_basin_path, \"-simulate\"])\n" + " subprocess.call([statemod_exe, simulated_basin_path, \"-simulate\"])" ] }, { @@ -4870,43 +4860,45 @@ "outputs": [], "source": [ "# SOW 1\n", - "output_dir= os.path.join(scenarios_dir, \"S0_1\")\n", + "output_dir = os.path.join(scenarios_dir, \"S0_1\")\n", "\n", "# path the the xre file\n", - "xre_file = os.path.join(output_dir, \"gm2015B_S0_1.xre\") \n", + "xre_file = os.path.join(output_dir, \"gm2015B_S0_1.xre\")\n", "\n", "# structure ID for reservoir of interest\n", - "structure_ID = '6203532' \n", + "structure_ID = \"6203532\"\n", "\n", "# name of the reservoir\n", - "structure_name = 'Blue_Mesa' \n", + "structure_name = \"Blue_Mesa\"\n", "\n", "# extract the target info into a Pandas data frame\n", - "df = stm.extract_xre_data(structure_name=structure_name,\n", - " structure_id=structure_ID,\n", - " input_file=xre_file,\n", - " basin_name=None,\n", - " output_directory=output_dir,\n", - " write_csv=True,\n", - " write_parquet=None\n", + "df = stm.extract_xre_data(\n", + " structure_name=structure_name,\n", + " structure_id=structure_ID,\n", + " input_file=xre_file,\n", + " basin_name=None,\n", + " output_directory=output_dir,\n", + " write_csv=True,\n", + " write_parquet=None,\n", ")\n", "\n", "# SOW 2\n", - "output_dir= os.path.join(scenarios_dir, \"S1_1\")\n", + "output_dir = os.path.join(scenarios_dir, \"S1_1\")\n", "\n", "# path the the xre file\n", - "xre_file = os.path.join(output_dir, \"gm2015B_S1_1.xre\") \n", + "xre_file = os.path.join(output_dir, \"gm2015B_S1_1.xre\")\n", "\n", "\n", "# extract the target info into a Pandas data frame\n", - "df = stm.extract_xre_data(structure_name=structure_name,\n", - " structure_id=structure_ID,\n", - " input_file=xre_file,\n", - " basin_name=None,\n", - " output_directory=output_dir,\n", - " write_csv=True,\n", - " write_parquet=None\n", - ")\n" + "df = stm.extract_xre_data(\n", + " structure_name=structure_name,\n", + " structure_id=structure_ID,\n", + " input_file=xre_file,\n", + " basin_name=None,\n", + " output_directory=output_dir,\n", + " write_csv=True,\n", + " write_parquet=None,\n", + ")" ] }, { @@ -4951,32 +4943,34 @@ "historic_res_dir = os.path.join(data_dir, \"historic_reservoir_levels\")\n", "blue_mesa_file = os.path.join(historic_res_dir, \"Blue_Mesa_xre_data.csv\")\n", "\n", - "# Import baseline dataframe \n", - "baseline = pd.read_csv(blue_mesa_file, index_col=False, usecols=['Year','Init. Storage']) \n", - "baseline = baseline.groupby('Year').mean().reset_index()\n", + "# Import baseline dataframe\n", + "baseline = pd.read_csv(\n", + " blue_mesa_file, index_col=False, usecols=[\"Year\", \"Init. Storage\"]\n", + ")\n", + "baseline = baseline.groupby(\"Year\").mean().reset_index()\n", "\n", "# Import SOW1\n", "s0_1_file = os.path.join(scenarios_dir, \"S0_1\", \"Blue_Mesa_xre_data.csv\")\n", - "SOW1 = pd.read_csv(s0_1_file, index_col=False, usecols=['Year','Init. Storage']) \n", - "SOW1 = SOW1.groupby('Year').mean().reset_index()\n", - " \n", + "SOW1 = pd.read_csv(s0_1_file, index_col=False, usecols=[\"Year\", \"Init. Storage\"])\n", + "SOW1 = SOW1.groupby(\"Year\").mean().reset_index()\n", + "\n", "# Import SOW2\n", "s1_1_file = os.path.join(scenarios_dir, \"S1_1\", \"Blue_Mesa_xre_data.csv\")\n", - "SOW2 = pd.read_csv(s1_1_file, index_col=False, usecols=['Year','Init. Storage']) \n", - "SOW2 = SOW2.groupby('Year').mean().reset_index()\n", - " \n", - "# Plot reservoir levels \n", + "SOW2 = pd.read_csv(s1_1_file, index_col=False, usecols=[\"Year\", \"Init. Storage\"])\n", + "SOW2 = SOW2.groupby(\"Year\").mean().reset_index()\n", + "\n", + "# Plot reservoir levels\n", "fig, ax = plt.subplots()\n", "\n", - "plt.plot(baseline['Year'], baseline['Init. Storage'],label='Baseline')\n", - "plt.plot(SOW1['Year'], SOW1['Init. Storage'],label='Reduced Evaporation')\n", - "plt.plot(SOW2['Year'], SOW2['Init. Storage'],label='Increased Evaporation')\n", + "plt.plot(baseline[\"Year\"], baseline[\"Init. Storage\"], label=\"Baseline\")\n", + "plt.plot(SOW1[\"Year\"], SOW1[\"Init. Storage\"], label=\"Reduced Evaporation\")\n", + "plt.plot(SOW2[\"Year\"], SOW2[\"Init. Storage\"], label=\"Increased Evaporation\")\n", "\n", "plt.title(\"Blue Mesa Storage\")\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"Reservoir Storage (AF)\")\n", "\n", - "plt.legend()\n" + "plt.legend()" ] }, { diff --git a/notebooks/N3_Reservoir_File_Modification.ipynb b/notebooks/N3_Reservoir_File_Modification.ipynb index f440f62..1ea925a 100644 --- a/notebooks/N3_Reservoir_File_Modification.ipynb +++ b/notebooks/N3_Reservoir_File_Modification.ipynb @@ -37,12 +37,13 @@ "import logging\n", "import os\n", "import pickle\n", - "from string import Template\n", "import subprocess\n", + "from string import Template\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", - "import pandas as pd \n", + "import pandas as pd\n", + "\n", "import statemodify as stm" ] }, @@ -88,12 +89,7 @@ "statemod_exe = os.path.join(root_dir, \"statemod\")\n", "\n", "# data directory and root name for the target basin\n", - "data_dir = os.path.join(\n", - " home_dir,\n", - " \"data\",\n", - " \"cm2015_StateMod\",\n", - " \"StateMod\"\n", - ")\n", + "data_dir = os.path.join(home_dir, \"data\", \"cm2015_StateMod\", \"StateMod\")\n", "\n", "# directory to the target basin input files with root name for the basin\n", "basin_path = os.path.join(data_dir, \"cm2015B\")\n", @@ -106,11 +102,7 @@ "\n", "\n", "# path to res template file\n", - "res_template_file = os.path.join(\n", - " home_dir,\n", - " \"data\",\n", - " \"cm2015B_template_res.rsp\"\n", - ")" + "res_template_file = os.path.join(home_dir, \"data\", \"cm2015B_template_res.rsp\")" ] }, { @@ -2911,9 +2903,11 @@ } ], "source": [ - "# Change directories first \n", - "os.chdir(data_dir) #This is needed specific to the Upper Colorado model as the path name is too long for the model to accept\n", - "subprocess.call([statemod_exe, \"cm2015B\", \"-simulate\"])\n" + "# Change directories first\n", + "os.chdir(\n", + " data_dir\n", + ") # This is needed specific to the Upper Colorado model as the path name is too long for the model to accept\n", + "subprocess.call([statemod_exe, \"cm2015B\", \"-simulate\"])" ] }, { @@ -2943,9 +2937,9 @@ } ], "source": [ - "#Extract shortages using statemodify convert_xdd() function \n", + "# Extract shortages using statemodify convert_xdd() function\n", "\n", - "# create a directory to store the historic shortages \n", + "# create a directory to store the historic shortages\n", "output_dir = os.path.join(data_dir, \"historic_shortages\")\n", "\n", "# create a directory to store the new files in if it does not exist\n", @@ -2961,12 +2955,11 @@ " # path, glob, or a list of paths/globs to the .xdd files you want to convert\n", " xdd_files=os.path.join(data_dir, \"*.xdd\"),\n", " # if the output .parquet files should only contain a subset of structure ids, list them here; None for all\n", - " id_subset=['3601008'],\n", + " id_subset=[\"3601008\"],\n", " # how many .xdd files to convert in parallel; optimally you will want 2-4 CPUs per parallel process\n", " parallel_jobs=2,\n", " # convert to natural data types\n", - " preserve_string_dtype=False\n", - " \n", + " preserve_string_dtype=False,\n", ")" ] }, @@ -3008,13 +3001,12 @@ } ], "source": [ - "data=pd.read_parquet(output_dir +'/cm2015B.parquet',engine='pyarrow')\n", + "data = pd.read_parquet(output_dir + \"/cm2015B.parquet\", engine=\"pyarrow\")\n", "\n", "fig, ax = plt.subplots()\n", "\n", - "for name, group in data.groupby('structure_id'):\n", - " ax.scatter(\n", - " group['year'], group['shortage_total'], label=name)\n", + "for name, group in data.groupby(\"structure_id\"):\n", + " ax.scatter(group[\"year\"], group[\"shortage_total\"], label=name)\n", "\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"Shortage (AF)\")\n", @@ -3069,14 +3061,16 @@ "# number of samples to generate\n", "n_samples = 1\n", "\n", - "stm.modify_res(output_dir=output_directory,\n", - " scenario=scenario,\n", - " basin_name=basin_name,\n", - " target_structure_id_list=['3603575','3604512'],\n", - " seed_value=seed_value,\n", - " n_jobs=n_jobs,\n", - " n_samples=n_samples,\n", - " save_sample=True)\n" + "stm.modify_res(\n", + " output_dir=output_directory,\n", + " scenario=scenario,\n", + " basin_name=basin_name,\n", + " target_structure_id_list=[\"3603575\", \"3604512\"],\n", + " seed_value=seed_value,\n", + " n_jobs=n_jobs,\n", + " n_samples=n_samples,\n", + " save_sample=True,\n", + ")" ] }, { @@ -3108,7 +3102,8 @@ ], "source": [ "import numpy as np\n", - "sample_array = np.load(output_directory+'/res_1-samples_scenario-1.npy')\n", + "\n", + "sample_array = np.load(output_directory + \"/res_1-samples_scenario-1.npy\")\n", "sample_array" ] }, @@ -5963,33 +5958,31 @@ "\n", "# read RSP template\n", "with open(res_template_file) as template_obj:\n", - " \n", " # read in file\n", " template_rsp = Template(template_obj.read())\n", "\n", " for i in sample:\n", - " \n", " # create scenario name\n", " scenario = f\"S{i}_{realization}\"\n", - " \n", + "\n", " # dictionary holding search keys and replacement values to update the template file\n", " d = {\"RES\": f\"../../input_files/cm2015B_{scenario}.res\"}\n", - " \n", + "\n", " # update the template\n", " new_rsp = template_rsp.safe_substitute(d)\n", - " \n", + "\n", " # construct simulated scenario directory\n", " simulated_scenario_dir = os.path.join(scenarios_dir_res, scenario)\n", " if not os.path.exists(simulated_scenario_dir):\n", " os.makedirs(simulated_scenario_dir)\n", - " \n", + "\n", " # target rsp file\n", " rsp_file = os.path.join(simulated_scenario_dir, f\"cm2015B_{scenario}.rsp\")\n", - " \n", + "\n", " # write updated rsp file\n", " with open(rsp_file, \"w\") as f1:\n", " f1.write(new_rsp)\n", - " \n", + "\n", " # construct simulated basin path\n", " simulated_basin_path = f\"cm2015B_{scenario}\"\n", "\n", @@ -5998,22 +5991,27 @@ " os.chdir(simulated_scenario_dir)\n", "\n", " subprocess.call([statemod_exe, simulated_basin_path, \"-simulate\"])\n", - " \n", - " #Save output to parquet files \n", - " print('creating parquet for ' + scenario)\n", - " \n", - " output_directory = os.path.join(parquet_dir_res+\"/scenario/\"+ scenario)\n", - " \n", + "\n", + " # Save output to parquet files\n", + " print(\"creating parquet for \" + scenario)\n", + "\n", + " output_directory = os.path.join(parquet_dir_res + \"/scenario/\" + scenario)\n", + "\n", " if not os.path.exists(output_directory):\n", " os.makedirs(output_directory)\n", - " \n", + "\n", " stm.xdd.convert_xdd(\n", " output_path=output_directory,\n", " allow_overwrite=True,\n", - " xdd_files=scenarios_dir_res + \"/\"+ scenario + \"/cm2015B_\"+scenario+\".xdd\",\n", - " id_subset=['3601008'],\n", + " xdd_files=scenarios_dir_res\n", + " + \"/\"\n", + " + scenario\n", + " + \"/cm2015B_\"\n", + " + scenario\n", + " + \".xdd\",\n", + " id_subset=[\"3601008\"],\n", " parallel_jobs=2,\n", - " preserve_string_dtype=False\n", + " preserve_string_dtype=False,\n", " )" ] }, @@ -6034,8 +6032,12 @@ }, "outputs": [], "source": [ - "baseline=pd.read_parquet(data_dir+'/'+'historic_shortages/cm2015B.parquet',engine='pyarrow')\n", - "SOW_1=pd.read_parquet(parquet_dir_res+'/scenario/S0_1/cm2015B_S0_1.parquet',engine='pyarrow')\n" + "baseline = pd.read_parquet(\n", + " data_dir + \"/\" + \"historic_shortages/cm2015B.parquet\", engine=\"pyarrow\"\n", + ")\n", + "SOW_1 = pd.read_parquet(\n", + " parquet_dir_res + \"/scenario/S0_1/cm2015B_S0_1.parquet\", engine=\"pyarrow\"\n", + ")" ] }, { @@ -6109,18 +6111,24 @@ } ], "source": [ - "baseline=pd.read_parquet(data_dir+'/'+'historic_shortages/cm2015B.parquet',engine='pyarrow')\n", - "SOW_1=pd.read_parquet(parquet_dir_res+'/scenario/S0_1/cm2015B_S0_1.parquet',engine='pyarrow')\n", + "baseline = pd.read_parquet(\n", + " data_dir + \"/\" + \"historic_shortages/cm2015B.parquet\", engine=\"pyarrow\"\n", + ")\n", + "SOW_1 = pd.read_parquet(\n", + " parquet_dir_res + \"/scenario/S0_1/cm2015B_S0_1.parquet\", engine=\"pyarrow\"\n", + ")\n", "\n", - "#Subtract shortages with respect to the baseline\n", - "subset_df=pd.concat([baseline['year'],baseline['shortage_total'],SOW_1['shortage_total']],axis=1)\n", - "subset_df = subset_df.set_axis(['Year', 'Baseline', 'SOW_1'], axis=1)\n", - "subset_df['SOW_1_diff']=subset_df['SOW_1']-subset_df['Baseline']\n", + "# Subtract shortages with respect to the baseline\n", + "subset_df = pd.concat(\n", + " [baseline[\"year\"], baseline[\"shortage_total\"], SOW_1[\"shortage_total\"]], axis=1\n", + ")\n", + "subset_df = subset_df.set_axis([\"Year\", \"Baseline\", \"SOW_1\"], axis=1)\n", + "subset_df[\"SOW_1_diff\"] = subset_df[\"SOW_1\"] - subset_df[\"Baseline\"]\n", "\n", - "#Plot shortages\n", + "# Plot shortages\n", "fig, ax = plt.subplots()\n", "\n", - "ax.scatter(subset_df['Year'], subset_df['SOW_1_diff'])\n", + "ax.scatter(subset_df[\"Year\"], subset_df[\"SOW_1_diff\"])\n", "\n", "plt.xlabel(\"Year\")\n", "plt.ylabel(\"Shortage (AF)\")\n", diff --git a/notebooks/N4_Streamflow_File_Modification.ipynb b/notebooks/N4_Streamflow_File_Modification.ipynb index 0d34272..4aaaf4e 100644 --- a/notebooks/N4_Streamflow_File_Modification.ipynb +++ b/notebooks/N4_Streamflow_File_Modification.ipynb @@ -69,12 +69,13 @@ "import logging\n", "import os\n", "import pickle\n", - "from string import Template\n", "import subprocess\n", + "from string import Template\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", - "import pandas as pd \n", + "import pandas as pd\n", + "\n", "import statemodify as stm" ] }, @@ -108,12 +109,7 @@ "statemod_exe = os.path.join(root_dir, \"statemod\")\n", "\n", "# data directory and root name for the target basin\n", - "data_dir = os.path.join(\n", - " home_dir,\n", - " \"data\",\n", - " \"cm2015_StateMod\",\n", - " \"StateMod\"\n", - ")\n", + "data_dir = os.path.join(home_dir, \"data\", \"cm2015_StateMod\", \"StateMod\")\n", "\n", "# directory to the target basin input files with root name for the basin\n", "basin_path = os.path.join(data_dir, \"cm2015B\")\n", @@ -121,12 +117,8 @@ "# scenarios output directory\n", "scenarios_dir = os.path.join(data_dir, \"scenarios\")\n", "\n", - "# path to iwr/xbm file \n", - "xbm_iwr_template_file = os.path.join(\n", - " home_dir,\n", - " \"data\",\n", - " \"cm2015B_template_xbm_iwr.rsp\"\n", - ")" + "# path to iwr/xbm file\n", + "xbm_iwr_template_file = os.path.join(home_dir, \"data\", \"cm2015B_template_xbm_iwr.rsp\")" ] }, { @@ -146,7 +138,7 @@ }, "outputs": [], "source": [ - "#Make directory to store HMM parameters\n", + "# Make directory to store HMM parameters\n", "\n", "output_dir = os.path.join(data_dir, \"HMM_parameters\")\n", "\n", @@ -158,9 +150,9 @@ "# choice to save parameters to NumPy array files\n", "save_parameters = True\n", "\n", - "fit_array_dict = stm.hmm_multisite_fit(n_basins=n_basins,\n", - " save_parameters=save_parameters,\n", - " output_directory=output_dir)\n", + "fit_array_dict = stm.hmm_multisite_fit(\n", + " n_basins=n_basins, save_parameters=save_parameters, output_directory=output_dir\n", + ")\n", "\n", "# unpack output dictionary\n", "unconditional_dry = fit_array_dict[\"unconditional_dry\"]\n", @@ -170,7 +162,7 @@ "covariance_matrix_wet = fit_array_dict[\"covariance_matrix_wet\"]\n", "covariance_matrix_dry = fit_array_dict[\"covariance_matrix_dry\"]\n", "wet_state_means = fit_array_dict[\"wet_state_means\"]\n", - "dry_state_means = fit_array_dict[\"dry_state_means\"]\n" + "dry_state_means = fit_array_dict[\"dry_state_means\"]" ] }, { @@ -198,7 +190,7 @@ }, "outputs": [], "source": [ - "#Create a folder to store the runs\n", + "# Create a folder to store the runs\n", "\n", "output_dir = os.path.join(data_dir, \"HMM_Runs\")\n", "\n", @@ -206,15 +198,17 @@ " os.makedirs(output_dir)\n", "\n", "# using the outputs of the fit function above; this function write output sample files to the output directory\n", - "stm.hmm_multisite_sample(logAnnualQ_h,\n", - " transition_matrix,\n", - " unconditional_dry,\n", - " dry_state_means,\n", - " wet_state_means,\n", - " covariance_matrix_dry,\n", - " covariance_matrix_wet,\n", - " n_basins=n_basins,\n", - " output_directory=output_dir)" + "stm.hmm_multisite_sample(\n", + " logAnnualQ_h,\n", + " transition_matrix,\n", + " unconditional_dry,\n", + " dry_state_means,\n", + " wet_state_means,\n", + " covariance_matrix_dry,\n", + " covariance_matrix_wet,\n", + " n_basins=n_basins,\n", + " output_directory=output_dir,\n", + ")" ] }, { @@ -245,10 +239,13 @@ } ], "source": [ - "stm.plot_flow_duration_curves(flow_realizations_directory=output_dir,\n", - " save_figure=True,output_directory=output_dir,\n", - " figure_name= 'FDC',\n", - " dpi= 300)" + "stm.plot_flow_duration_curves(\n", + " flow_realizations_directory=output_dir,\n", + " save_figure=True,\n", + " output_directory=output_dir,\n", + " figure_name=\"FDC\",\n", + " dpi=300,\n", + ")" ] }, { @@ -278,16 +275,16 @@ }, "outputs": [], "source": [ - "#Make directory to store input files \n", + "# Make directory to store input files\n", "\n", "output_dir = os.path.join(data_dir, \"input_files\")\n", "\n", "if not os.path.exists(output_dir):\n", " os.makedirs(output_dir)\n", - " \n", + "\n", "\n", "flow_realizations_directory = os.path.join(data_dir, \"HMM_Runs\")\n", - " \n", + "\n", "scenario = \"1\"\n", "\n", "# basin name to process\n", @@ -303,15 +300,17 @@ "n_samples = 1\n", "\n", "# generate a batch of files using generated LHS\n", - "stm.modify_xbm_iwr(output_dir=output_dir,\n", - " flow_realizations_directory=flow_realizations_directory,\n", - " scenario=scenario,\n", - " basin_name=basin_name,\n", - " seed_value=seed_value,\n", - " n_jobs=n_jobs,\n", - " n_samples=n_samples,\n", - " save_sample=True,\n", - " randomly_select_flow_sample=True) " + "stm.modify_xbm_iwr(\n", + " output_dir=output_dir,\n", + " flow_realizations_directory=flow_realizations_directory,\n", + " scenario=scenario,\n", + " basin_name=basin_name,\n", + " seed_value=seed_value,\n", + " n_jobs=n_jobs,\n", + " n_samples=n_samples,\n", + " save_sample=True,\n", + " randomly_select_flow_sample=True,\n", + ")" ] }, { @@ -3128,33 +3127,34 @@ "\n", "# read RSP template\n", "with open(xbm_iwr_template_file) as template_obj:\n", - " \n", " # read in file\n", " template_rsp = Template(template_obj.read())\n", "\n", " for i in sample:\n", - " \n", " # create scenario name\n", " scenario = f\"S{i}_{realization}\"\n", - " \n", + "\n", " # dictionary holding search keys and replacement values to update the template file\n", - " d = {\"XBM\": f\"../../input_files/cm2015B_{scenario}.xbm\",\"IWR\": f\"../../input_files/cm2015B_{scenario}.iwr\"}\n", - " \n", + " d = {\n", + " \"XBM\": f\"../../input_files/cm2015B_{scenario}.xbm\",\n", + " \"IWR\": f\"../../input_files/cm2015B_{scenario}.iwr\",\n", + " }\n", + "\n", " # update the template\n", " new_rsp = template_rsp.safe_substitute(d)\n", - " \n", + "\n", " # construct simulated scenario directory\n", " simulated_scenario_dir = os.path.join(scenarios_dir, scenario)\n", " if not os.path.exists(simulated_scenario_dir):\n", " os.makedirs(simulated_scenario_dir)\n", - " \n", + "\n", " # target rsp file\n", " rsp_file = os.path.join(simulated_scenario_dir, f\"cm2015B_{scenario}.rsp\")\n", - " \n", + "\n", " # write updated rsp file\n", " with open(rsp_file, \"w\") as f1:\n", " f1.write(new_rsp)\n", - " \n", + "\n", " # construct simulated basin path\n", " simulated_basin_path = f\"cm2015B_{scenario}\"\n", "\n", @@ -3162,7 +3162,7 @@ " print(f\"Running: {scenario}\")\n", " os.chdir(simulated_scenario_dir)\n", "\n", - " subprocess.call([statemod_exe, simulated_basin_path, \"-simulate\"])\n" + " subprocess.call([statemod_exe, simulated_basin_path, \"-simulate\"])" ] }, { @@ -3194,14 +3194,16 @@ ], "source": [ "# Example with Granby Lake\n", - "zip_file_path = os.path.join(home_dir, 'data', 'Granby_Dataset.zip')\n", - "final_directory = os.path.join(home_dir, 'data/')\n", + "zip_file_path = os.path.join(home_dir, \"data\", \"Granby_Dataset.zip\")\n", + "final_directory = os.path.join(home_dir, \"data/\")\n", "\n", - "!unzip $zip_file_path -d $final_directory \n", - "granby_hmm, granby_hist, granby_hist_mean, granby_hist_1p = stm.read_xre(os.path.join(home_dir,\"data/Upper_Colorado/\"), 'Granby')\n", + "!unzip $zip_file_path -d $final_directory\n", + "granby_hmm, granby_hist, granby_hist_mean, granby_hist_1p = stm.read_xre(\n", + " os.path.join(home_dir, \"data/Upper_Colorado/\"), \"Granby\"\n", + ")\n", "\n", "# Plot quantiles\n", - "stm.plot_res_quantiles(granby_hmm, granby_hist_mean, 'Lake Granby')" + "stm.plot_res_quantiles(granby_hmm, granby_hist_mean, \"Lake Granby\")" ] }, { @@ -3241,7 +3243,7 @@ ], "source": [ "# using the output of the above `read_xre` function as inputs\n", - "stm.plot_reservoir_boxes(granby_hmm, granby_hist, 'Lake Granby')" + "stm.plot_reservoir_boxes(granby_hmm, granby_hist, \"Lake Granby\")" ] }, { diff --git a/notebooks/N5_Batch_Modification.ipynb b/notebooks/N5_Batch_Modification.ipynb index 6ef477a..5ed2e07 100644 --- a/notebooks/N5_Batch_Modification.ipynb +++ b/notebooks/N5_Batch_Modification.ipynb @@ -33,12 +33,13 @@ "import logging\n", "import os\n", "import pickle\n", - "from string import Template\n", "import subprocess\n", + "from string import Template\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", - "import pandas as pd \n", + "import pandas as pd\n", + "\n", "import statemodify as stm" ] }, @@ -80,12 +81,7 @@ "statemod_exe = os.path.join(root_dir, \"statemod\")\n", "\n", "# data directory and root name for the target basin\n", - "data_dir = os.path.join(\n", - " home_dir,\n", - " \"data\",\n", - " \"cm2015_StateMod\",\n", - " \"StateMod\"\n", - ")\n", + "data_dir = os.path.join(home_dir, \"data\", \"cm2015_StateMod\", \"StateMod\")\n", "\n", "# directory to the target basin input files with root name for the basin\n", "basin_path = os.path.join(data_dir, \"cm2015B\")\n", @@ -93,16 +89,12 @@ "# scenarios output directory\n", "scenarios_dir = os.path.join(data_dir, \"scenarios\")\n", "\n", - "#parquet output directory\n", - "parquet_dir=os.path.join(data_dir, \"parquet\")\n", + "# parquet output directory\n", + "parquet_dir = os.path.join(data_dir, \"parquet\")\n", "\n", "\n", "# path to template file\n", - "multi_template_file = os.path.join(\n", - " home_dir,\n", - " \"data\",\n", - " \"cm2015B_template_multi.rsp\"\n", - ")" + "multi_template_file = os.path.join(home_dir, \"data\", \"cm2015B_template_multi.rsp\")" ] }, { @@ -167,13 +159,9 @@ "# problem dictionary\n", "problem_dict = {\n", " \"n_samples\": 1,\n", - " 'num_vars': 3,\n", - " 'names': ['modify_eva', 'modify_ddr', 'modify_ddm'],\n", - " 'bounds': [\n", - " [-0.5, 1.0],\n", - " None,\n", - " [0.5, 1.0]\n", - " ],\n", + " \"num_vars\": 3,\n", + " \"names\": [\"modify_eva\", \"modify_ddr\", \"modify_ddm\"],\n", + " \"bounds\": [[-0.5, 1.0], None, [0.5, 1.0]],\n", " # additional settings for each function\n", " \"modify_eva\": {\n", " \"seed_value\": seed_value,\n", @@ -181,7 +169,7 @@ " \"scenario\": scenario,\n", " \"basin_name\": basin_name,\n", " \"query_field\": \"id\",\n", - " \"ids\": [\"10008\", \"10009\"]\n", + " \"ids\": [\"10008\", \"10009\"],\n", " },\n", " \"modify_ddr\": {\n", " \"seed_value\": seed_value,\n", @@ -191,7 +179,7 @@ " \"query_field\": \"id\",\n", " \"ids\": [\"3600507.01\", \"3600507.02\"],\n", " \"admin\": [1, None],\n", - " \"on_off\": [1, 1]\n", + " \"on_off\": [1, 1],\n", " },\n", " \"modify_ddm\": {\n", " \"seed_value\": seed_value,\n", @@ -199,8 +187,8 @@ " \"scenario\": scenario,\n", " \"basin_name\": basin_name,\n", " \"query_field\": \"id\",\n", - " \"ids\": [\"3600507\", \"3600603\"]\n", - " }\n", + " \"ids\": [\"3600507\", \"3600603\"],\n", + " },\n", "}\n", "\n", "# run in batch\n", @@ -511,33 +499,35 @@ "\n", "# read RSP template\n", "with open(multi_template_file) as template_obj:\n", - " \n", " # read in file\n", " template_rsp = Template(template_obj.read())\n", "\n", " for i in sample:\n", - " \n", " # create scenario name\n", " scenario = f\"S{i}_{realization}\"\n", - " \n", + "\n", " # dictionary holding search keys and replacement values to update the template file\n", - " d = {\"EVA\": f\"../../input_files/cm2015B_{scenario}.eva\",\"DDM\": f\"../../input_files/cm2015B_{scenario}.ddm\",\"DDR\": f\"../../input_files/cm2015B_{scenario}.ddr\"}\n", - " \n", + " d = {\n", + " \"EVA\": f\"../../input_files/cm2015B_{scenario}.eva\",\n", + " \"DDM\": f\"../../input_files/cm2015B_{scenario}.ddm\",\n", + " \"DDR\": f\"../../input_files/cm2015B_{scenario}.ddr\",\n", + " }\n", + "\n", " # update the template\n", " new_rsp = template_rsp.safe_substitute(d)\n", - " \n", + "\n", " # construct simulated scenario directory\n", " simulated_scenario_dir = os.path.join(scenarios_dir, scenario)\n", " if not os.path.exists(simulated_scenario_dir):\n", " os.makedirs(simulated_scenario_dir)\n", - " \n", + "\n", " # target rsp file\n", " rsp_file = os.path.join(simulated_scenario_dir, f\"cm2015B_{scenario}.rsp\")\n", - " \n", + "\n", " # write updated rsp file\n", " with open(rsp_file, \"w\") as f1:\n", " f1.write(new_rsp)\n", - " \n", + "\n", " # construct simulated basin path\n", " simulated_basin_path = f\"cm2015B_{scenario}\"\n", "\n", @@ -546,16 +536,22 @@ " os.chdir(simulated_scenario_dir)\n", "\n", " subprocess.call([statemod_exe, simulated_basin_path, \"-simulate\"])\n", - " \n", - " #Save output to parquet files \n", - " print('creating parquet for ' + scenario)\n", - " \n", - " output_directory = os.path.join(parquet_dir+ \"/scenario/\"+ scenario)\n", - " \n", + "\n", + " # Save output to parquet files\n", + " print(\"creating parquet for \" + scenario)\n", + "\n", + " output_directory = os.path.join(parquet_dir + \"/scenario/\" + scenario)\n", + "\n", " if not os.path.exists(output_directory):\n", " os.makedirs(output_directory)\n", - " \n", - " stm.xdd.convert_xdd(output_path=output_directory,allow_overwrite=False,xdd_files=scenarios_dir + \"/\"+ scenario + \"/cm2015B_\"+scenario+\".xdd\",id_subset=['3601008'],parallel_jobs=2)\n" + "\n", + " stm.xdd.convert_xdd(\n", + " output_path=output_directory,\n", + " allow_overwrite=False,\n", + " xdd_files=scenarios_dir + \"/\" + scenario + \"/cm2015B_\" + scenario + \".xdd\",\n", + " id_subset=[\"3601008\"],\n", + " parallel_jobs=2,\n", + " )" ] }, { diff --git a/statemodify/xdd.py b/statemodify/xdd.py index 4355581..e59a3cb 100644 --- a/statemodify/xdd.py +++ b/statemodify/xdd.py @@ -21,7 +21,7 @@ def __init__( xdd_files: Union[str, Path, list[Union[str, Path]]] = "**/*.xdd", id_subset: Union[None, list[str]] = None, parallel_jobs: int = 4, - preserve_string_dtype: bool = True + preserve_string_dtype: bool = True, ): """Convert object for transforming StateMod output .xdd files into compressed, columnar .parquet files. @@ -272,7 +272,7 @@ def _parse_file( output_path: str, id_subset: Union[None, list[str]], preserve_string_dtype: bool, - field_dtypes: dict + field_dtypes: dict, ): data = [] with open(file) as f: @@ -297,7 +297,7 @@ def _parse_file( if preserve_string_dtype is False: df = df.astype(field_dtypes) - df.to_parquet(f"{output_path}/{Path(file).stem}.parquet") + df.to_parquet(f"{output_path}/{Path(file).stem}.parquet") def convert_xdd( @@ -307,7 +307,7 @@ def convert_xdd( xdd_files: Union[str, Path, list[Union[str, Path]]] = "**/*.xdd", id_subset: Union[None, list[str]] = None, parallel_jobs: int = 4, - preserve_string_dtype: bool = True + preserve_string_dtype: bool = True, ): """Convert StateMod output .xdd files to compressed, columnar .parquet files.