Content:

  • Installation

  • Preparing model input parameters and dataset

  • Preparing DRYP simulations

  • Post processing model outputs

Preparing DRYP model simulations

  1. Preparing model parameter and setting files

  • Running DRYP

  • Preparing multiple simulation

  • Runing pipeline DRYP

[1]:
print("Local_directory\n|__training\n      |__basin\n      |   |__datasets\n      |   |    |__csv\n      |   |    |__shp\n      |   |__model\n      |   |    |__inputs\n      |   |__outputs\n      |__regional")
Local_directory
|__training
      |__basin
      |   |__datasets
      |   |    |__csv
      |   |    |__shp
      |   |__model
      |   |    |__inputs
      |   |__outputs
      |__regional

Specifiy the location of training files, once the training files are downloaded change the paths below to access them through this traning material.

[24]:

training_general_path = "D:/HAD/training/historical/" regional_path = "D:/HAD/training/historical/regional/" basin_path = "D:/HAD/training/historical/basin/"

3.1. Preparing model parameter and setting files

DRYP requires at leas two files for running, this is the input-paramter-file and the setting-paramter file. Additional files may be required depending on the model settings and dataset properties.

The folowing are the list of DRYP model files:

  • input parameter file (essential)

  • setting-paramter file (essential)

  • GW parameter file (optional)

  • riparian input paramter file (optional)

  • data projection setting file (optional)

Model input and settins files are simple plain text format files, so they can be easily updated using any text editor.

For conveniency we use here the library pandas to edit these files.

[3]:
import pandas as pd
  • Input parameter file

The input-parameter-file is used for specifying path names of all names model files, including input parameters, forcings, and setting files, as well as model outputs. You can also specified the model name as well as any additional model file.

[ ]:
default_input = {

   "model_name": "DRYP_model_sim",
   "path_input": "",

   "TERRAIN": {
      "path_dem": "required", # Digital elevation model
      "path_Qo": None, # initial channel storage
      "path_fdl": None, # flow direction in landlab format
      "path_riv_decay": None, # river flow velocity
      "path_mask": None, # basin mask. model active domain
      "path_riv_len": None, # river lenght
      "path_riv_width": None, # river width
      "path_riv_elev": None, # tiver stream bottom elevation
      "path_of_bc_flux": None, # flux boundary condition overland flow

   },

   "VEGETATION": {
      "path_veg_kc": None,
      "path_veg_lulc": None,
      #"path_veg_nn": None
   },

   "UNSATURATED": {
      "path_uz_theta_sat": None,
      "path_uz_theta_res": None,
      "path_uz_theta_awc": None,
      "path_uz_theta_wp": None,
      "path_uz_rootdepth": None,
      "path_uz_lambda": None,
      "path_uz_psi": None,
      "path_uz_ksat": None,
      "path_uz_sigmaksat": None,
      "path_uz_theta": None, # initial conditions
      "path_riv_ksat": None,
      "path_uz_bottomksat": None,
      "path_uz_bc_flux": None,
   },

   "SATURATED": {
      "path_sz_mask": None,
      "path_sz_ksat": None,
      "path_sz_sy": None,
      "path_sz_wte": None,
      "path_sz_bc_flux": None,
      "path_sz_bc_head": None,
      "path_sz_bottom": None,
      "path_sz_depth": None,
      "path_sz_bdd": None,
      "path_sz_type": None, # aquifer type
   },

   "METEO": {
      "path_pre": None,
      "path_pet": None,
      "path_aof": None,
            "path_lai": None,
            "path_savi": None,
            "path_kc": None,
    "path_TSOF": None,
    "path_TSUZ": None,
    "path_TSSZ": None,
    "path_TSav": None, # vegetation cover fraction
      "path_TSWB": None, # Variable flux for water bodies and dams
   },

   "OUTPUT": {
      "path_out_sz": None,
      "path_out_uz": None,
      "path_out_oz": None,
      "path_output": None,
      "path_setting": None,
      "path_store_settings": None
   },

   "RIPARIAN": {
      "path_rp_theta_sat": None,
      "path_rp_theta_res": None,
      "path_rp_theta_awc": None,
      "path_rp_theta_wp": None,
      "path_rp_rootdepth": None,
      "path_rp_lambda": None,
      "path_rp_psi": None,
      "path_rp_ksat": None,
      "path_rp_sigmaksat": None,
      "path_rp_theta": None,
      "path_rp_width": None
   },

   "INTERCEPTION": {
      "path_veg_lulc_frac": None,
      "path_veg_hs_savi": None,
      "path_veg_hs_par_a": None,
      "path_veg_hs_par_b": None,
      "path_veg_hs_savi_min": None,
      "path_veg_hs_savi_max": None,
      "path_veg_hs_lai": None,
      "path_veg_hs_tap_depth": None,
      "path_veg_hs_extinction_depth": None,
      "path_veg_hs_fcw": None,
      "path_veg_hs_sca": None,
      "path_veg_rp_par_a": None,
      "path_veg_rp_par_b": None,
      "path_veg_rp_savi_min": None,
      "path_veg_rp_savi_max": None,
      "path_veg_rp_lai": None,
      "path_veg_rp_tap_depth": None,
      "path_veg_rp_extinction_depth": None,
      "path_veg_rp_fcw": None,
      "path_veg_rp_sca": None
   },

   "GROUNDWATER": {
   #   "path_gw_depth": None,
   #   "path_gw_bdd": None,
      "path_gw_2l_bottom": None,
      "path_gw_2l_ksat": None,
      "path_gw_2l_sy": None,
      "path_gw_2l_ss": None,
      "path_gw_2l_wte": None,
   #   "path_gw_type": None, # aquifer type
   #   "path_gw_lake_elev": None, # lakes bathymetry
   #   "path_pnds_hmax": None, # ponds max depth
   #   "path_pnds_Amax": None, # ponds maximum extend
   },

   "WATER_BODIES": {
      "path_lake_depth": None, # lakes bathymetry
      "path_pnd_hmax": None, # ponds max depth
      "path_pnd_Amax": None, # ponds maximum extend
      "path_pnd_Vo": None, # ponds volume of water
      "path_wb_bc_flux": None, # water bodies boundary conditions
   }

}

  • Model setting file: Reference and projection system

[ ]:

default_settings = { "SIMULATION_PERIOD": { "start_date": "2000 1 1", "end_date": "2002 1 3", }, "PROJECTION" : None, "TIMESTEP_SETTINGS": { "dt_of": 60, "dt_uz": 60, "dt_gw": 60, }, "READING": { "data_reading": { "pre": 0, "pet": 0, "abs": 0, "kc": 0, "fluxOF": 0, "fluxUZ": 0, "fluxSZ": 0, "fluxWB": 0, "savi": 0, "savi_min": 0, "savi_max": 0, "lai": 0, "av": 0, }, "data_step": { "pre": 60, "pet": 60, "abs": 60, "kc": 60, "fluxOF": 60, "fluxUZ": 60, "fluxSZ": 60, "fluxWB": 60, "savi": 60, "savi_min": 60, "savi_max": 60, "lai": 60, "av": 60, }, "data_reproject": { "pre": True, "pet": True, "abs": True, "kc": True, "fluxOF": False, "fluxUZ": False, "fluxSZ": False, "fluxWB": False, "savi": True, "savi_min": True, "savi_max": True, "lai": True, "av": False, }, "data_interp": { "pre": True, "pet": True, "abs": True, "kc": True, "fluxOF": False, "fluxUZ": False, "fluxSZ": False, "fluxWB": False, "savi": True, "savi_min": True, "savi_max": True, "lai": True, "av": True, }, "data_projection" : { "pre": "EPSG:4326", "pet": "EPSG:4326", "abs": "EPSG:4326", "kc": "EPSG:4326", "fluxOF": "EPSG:4326", "fluxUZ": "EPSG:4326", "fluxSZ": "EPSG:4326", "fluxWB": "EPSG:4326", "savi":"EPSG:4326", "savi_min": "EPSG:4326", "savi_max": "EPSG:4326", "lai": "EPSG:4326", "av": "EPSG:4326", }, }, "COMPONENTS": { "method_inf": 1, # choose infiltration method "run_gw" : True, # activate groundwater component "method_gw": 0, # choose groundwater transimissivity approach }, "OUTPUT": { "output_csv": True, # activate save model outputs (only csv files) "output_grid": False, # activate save model outputs (grid files) "output_dt_csv": "1M", "output_dt": "1M", }, "GLOBAL_FACTORS": { "uz_kdt": 1.0, "uz_kdroot": 1.0, "uz_kawc": 1.0, "uz_kkast": 1.0, "uz_ksigma": 1.0, "riv_kksat": 1.0, "riv_kdecay": 1.0, "riv_kwidth": 1.0, "sz_kksat": 1.0, "sz_ksy": 1.0, "of_kflow": 1.0 } }
  • Changing model input paths

Creating a new model requires to specify the new paths of the model parameter files, it is a task that takes time and can be prone to errors. To avoid any mistake when changing paths, we can use the replace text option availble in python.

[8]:
from shutil import copyfile
import fileinput
import sys
[9]:
def replaceAll(file, searchExp, replaceExp):
    for line in fileinput.input(file, inplace=1):
        if searchExp in line:
            line = line.replace(searchExp,replaceExp)
        sys.stdout.write(line)
[10]:
fname_file = regional_path + "/model/HAD_IMERGcv_input_sim85.json"
fname_file_update = regional_path + "/model/HAD_IMERGcv_input_sim85_cuwalid.json"

Create a copy avoid the lost of information

[12]:
copyfile(fname_file, fname_file_update)
[12]:
'D:/HAD/training/historical/regional//model/HAD_IMERGcv_input_sim85_cuwalid.json'
[13]:
#pd.read_csv(fname_file_update)

Specify the which text you want to change and the new text

[14]:
searchExpresion = "/home/c1755103/HAD/input/"
replaceExpresion = "/home/cuwalid/HAD/input/"
#replaceAll(fname_file_update, searchExpresion, replaceExpresion)
[15]:
searchExpresion = "/home/c1755103/HAD/HAD_"
replaceExpresion = "/home/cuwalid/HAD/HAD_"
#replaceAll(fname_file_update, searchExpresion, replaceExpresion)

3.2. Runing DRYP

WARNING: If you are using Windows you have to use the docker for running DRYP

There are two main ways of runing DRYP, you can even add more option that are more convenient for yor analysis:

  1. Running DRYP from the command line

[ ]:
# python run_model_input.py <finput>
  1. Running dryp within the python environment

[ ]:
import sys
sys.path.append('C:/Users/Edisson/Documents/GitHub/DRYPv2.0.1')

from cuwalid.dryp.main_DRYP import run_DRYP
#run_DRYP(inputfile)

3.3. Preparing and running multiple DRYP simulations

This is useful for performing a Monte Carlo analysis or stochastic forecasting.

We can create a script to modify the model files that changes specific fiels of each of the all model files, below ther is an example of the script to genereate multiple simulation files for a calibration pourposes.

[16]:
import os
import json
import numpy as np
import calendar
import pandas as pd
[20]:
def write_JSON_dryp_files(json_template, destination, model_name=None, path_pre=None, path_pet=None,
                                               start_date=None, end_date=None, new_setting_file=None,
                                               path_Qo=None, path_uz_theta=None, path_sz_wte=None, path_rp_theta=None,
                                               path_pnd_Vo=None, path_outputs=None, parameter_factors=None):
    """ This function create the simulation and setting file for running DRYP. New
    files are created  based on files provided as original files, this function
    changes the model name, precipitation and potential evapotranspiration
    paths.
    WARNING: if no new filename is provided it will replace the original file

    Parameters:
    -----------
    json_template : string
                    model input, as a dictionary
    model_name : string
                    model name for the new file
    path_pre : string
                    precipitation dataset name, including path
    path_pet : string
                    potential evapotranspiration dataset name, including path
    start_date : string
                    date in the following format "YYYY-MM-DD" (e.g. 2002-01-01)
    end_date : string
                    date in the following format "YYYY-MM-DD" (e.g. 2002-03-01)
    new_setting_file : bool
                    If True it create the setting dryp file
    """

    # Change necesarry variables in template
    if model_name is not None:
            json_template["model_name"] = model_name
    if path_pre is not None:
            json_template["METEO"]["path_pre"] = path_pre
    if path_pet is not None:
            json_template["METEO"]["path_pet"] = path_pet

    if path_Qo is not None:
            json_template["TERRAIN"]["path_Qo"] = path_Qo
    if path_uz_theta is not None:
            json_template["UNSATURATED"]["path_uz_theta"] = path_uz_theta
    if path_sz_wte is not None:
            json_template["SATURATED"]["path_sz_wte"] = path_sz_wte
    if path_rp_theta is not None:
            json_template["RIPARIAN"]["path_rp_theta"] = path_rp_theta
    if path_pnd_Vo is not None:
            json_template["WATER_BODIES"]["path_pnd_Vo"] = path_pnd_Vo
    if path_outputs is not None:
            json_template["OUTPUT"]["path_output"] = path_outputs

    # create new settings file only if new setting file does not exist
    if new_setting_file is not None:

            # Get input file as dictionary
            settings_file_path = json_template["OUTPUT"]["path_setting"]
            with open(settings_file_path, 'r') as file:
                    settings_file_template = json.load(file)
            # Change settings file location
            json_template["OUTPUT"]["path_setting"] = new_setting_file

            ## Change variables in settings file
            if start_date is not None:
                    settings_file_template["SIMULATION_PERIOD"]["start_date"] = start_date
            if end_date is not None:
                    settings_file_template["SIMULATION_PERIOD"]["end_date"] = end_date

            # update parameter factors only if provided
            if parameter_factors is not None:
                    for ikey in parameter_factors.keys():
                            settings_file_template["GLOBAL_FACTORS"][ikey] = parameter_factors[ikey]

    # Save the `dryp` part to the destination file
    with open(destination, "w") as dest_file:
            json.dump(json_template, dest_file, indent=4)

    # Save the `dryp_settings` part to the new settings file
    if new_setting_file is not None:
            with open(new_setting_file, "w") as settings_file:
                    json.dump(settings_file_template, settings_file, indent=4)
[41]:
def gen_array_input_files(fname_input, fname_parameter_sets, model_name, nmax=100):
    # open input file
    with open(fname_input, 'r') as file:
            input_template = json.load(file)
    # read parameter sets
    parameter = pd.read_csv(fname_parameter_sets)

    #Create a copy of inputfile
    fname_root = fname_input.split('.')[0]
    fname_ext = fname_input.split('.')[1]
    for npar in range(0, nmax):
            # new input file name
            newfname_input = fname_root+'_'+str(npar)+'.json'
            imodel_name = model_name + "_" + str(npar)
            # replace all new values in dataset
            print(imodel_name)
            write_JSON_dryp_files(input_template, newfname_input, model_name=imodel_name,
                                            parameter_factors=dict(parameter.loc[npar]))

Write multiple input model files

[43]:
# ========================================================
# LOOP FOR CREATING MULTIPLE IMPORT FILES FOR RUNNING IN HPC
fname = [
#basin_path +'model/HAD_IMERGcv_input_sim85.json',
regional_path +'model/HAD_IMERGcv_input_sim85.json',
]
fname_parameter_sets = training_general_path + "/basin/dataset/csv/test_parameter_set.csv"
for ifname in fname:
    gen_array_input_files(ifname, fname_parameter_sets, "HAD_IMERGcv_input_sim85", nmax=2)
HAD_IMERGcv_input_sim85_0
HAD_IMERGcv_input_sim85_1

For running the script above you need to create a parameter set file, an example is shown below:

[ ]:
fname_parameter_sets = training_general_path + "/regional/datasets/csv/test_parameter_set.csv"
[ ]:
pd.read_csv(fname_parameter_sets).head(5)

TASK: modify the script to create model files for running multiples simulation with different forcing datasets (e.g. stochastic forecasting).

HINT: lines 66 and 68 needs to be modified in the input model parameters file, add lines in script to modify this lines

3.4. Runing pipeline DRYP simulations

For large simulations output files can become very large which may result in consuming cosiderable amount of memory that can even stop the model. This problem is very challenging to address in python, requiring significan changes in the code or even considering parallelisation. Here, we simply split a continuos simulation (in time, not in space) to store model output files in specified simulation periods (e.g. years).

For this, we can create a series of model files that read and save inital conditions for the previous and subsequente simulations. An example of a python script to generate pipeline simulaiton is presented below:

[45]:
def gen_inital_end_simulation_dates(year, month, day, dtyear=1, dtmonth=0, dtday=0):
    """This function gets the initial and end date of a specified period
    """
    date_ini = str(year+dtyear) + ' ' + str(month+dtmonth) + ' ' + str(day+dtday)
    date_end = str(year+dtyear+1) + ' ' + str(month+dtmonth) + ' ' + str(day+dtday)

    name = str(year)
    name_end = str(year+dtyear)

    return date_ini, date_end, name, name_end

def gen_pipeline_dryp_files(path_json_template, path_outputs=None, model_name=None, start_year=2000, end_year=2023):

    # Get input file as dictionary
    with open(path_json_template, 'r') as file:
            json_template = json.load(file)
    if path_outputs is None:
            path_outputs = json_template["OUTPUT"]["path_output"]

    # get path, filename and model name
    path, file = os.path.split(path_json_template)

    # identify filename of paramter file
    if model_name is None:
            model_name = file.split('.')[0]

    # avoid updating forcing dataset
    path_pre, path_pet = None, None

    # iterate over the simulatio period
    for i, iyear in enumerate(range(start_year, end_year)):
            mname = model_name + "_" + str(iyear)
            path_new_input = path + "/" + mname + ".json"
            path_new_settings = path + "/"  + mname + "_settings.json"
            # get name of initial conditions
            # skip if the initial year
            if i != 0:
                    mname0 = model_name + "_" + str(iyear-1)
                    path_Qo = path_outputs + '/'+ mname0 +'_avg_Q_ini.asc'
                    path_uz_theta = path_outputs + '/'+ mname0 +'_avg_tht_ini.asc'
                    path_sz_wte = path_outputs + '/'+ mname0 +'_avg_wte_ini.asc'
                    path_rp_theta = path_outputs + '/'+ mname0 + '_avg_tht_rp_ini.asc'
                    path_pnd_Vo = path_outputs + '/' + mname0 + '_avg_V_pnd_ini.asc'
            else:
                    path_Qo=None
                    path_uz_theta=None
                    path_sz_wte=None
                    path_rp_theta=None
                    path_pnd_Vo=None

            # get simulation period
            if i == 0:
                    date_ini, date_end = None, None
            else:
                    date_ini, date_end, name, name_end = gen_inital_end_simulation_dates(iyear-1, 1, 1)

            # write json files
            write_JSON_dryp_files(json_template, mname, path_pre, path_pet, path_new_input,
                                               start_date=date_ini, end_date=date_end, new_setting_file=path_new_settings,
                                               path_Qo=path_Qo,
                                               path_uz_theta=path_uz_theta,
                                               path_sz_wte=path_sz_wte,
                                               path_rp_theta=path_rp_theta,
                                               path_pnd_Vo=path_pnd_Vo
                                               )

Generate pipeline simulation files

[ ]:
# ========================================================
# LOOP FOR CREATING MULTIPLE IMPORT FILES FOR RUNNING IN HPC
fname = [
training_general_path + "model/HAD_IMERG_Tana_input_sim.json",
]

path_outputs = "/home/c1755103/HAD/HAD_output"
# WARNINGS
# The initial file should specify the initial conditions at the begining of
# the simulation of the entire period, therefore, the name of inital model
# must be according to the name required to the next simulation
# subsequent period must conside
for ifname_input in fname:
    gen_pipeline_dryp_files(ifname_input, path_outputs)
[ ]:
import glob
files = glob.glob(regional_path + "inputs/*.asc")
[ ]:
new_path = "D:/"
fnames = []
for ifile in files:
    fnames.append(ifile.split("\\")[-1])
    new_name = new_path + ifile.split("\\")[-1]
    print(new_name)
    #new_name = new_path + fnames[0]
    #process clip raster
[ ]:
#fnames
[ ]:
#new_path = "D:/"
#new_name = new_path + fnames[0]