Sensitivity_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6cd3f109-a5ec-4ce1-b7cf-8309b149e3a1",
   "metadata": {},
   "source": [
    "# Sensitivity analysis\n",
    "## Upper Olivares catchment - january, 2023\n",
    "#### Paul Sandoval Quilodrán - https://github.com/SQPaul/Olivares"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8b4f93a-7e9f-4426-9e44-1c66b681e5cb",
   "metadata": {},
   "source": [
    "#### import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c7fa67e7-b235-405d-a845-756cafb7daba",
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import os\n",
    "import time\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import openpyxl\n",
    "import datetime\n",
    "import hydroeval as he\n",
    "from configparser import ConfigParser\n",
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f6d0365-5048-4f86-8d37-3b6d0eb0650f",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1.- Define the functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f6318ba-7e4c-454f-9aae-14b49da2a4af",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Start time\n",
    "start_time = datetime.now()\n",
    "#Define some functions\n",
    "def run_sphy():\n",
    "    os.system(\"python -m SPHY.main_sw -i P:\\Projects\\Olivares\\SPHY\\cfg\\sw\\sphy_config_sw_SA.cfg\")\n",
    "    \n",
    "def performance(kge_r,corr,gamma,betha,nash_r):\n",
    "    #Import paths\n",
    "    obs_path = r\"P:\\Projects\\Olivares\\Streamflow\\qobs_2011.csv\"\n",
    "    sim_path = r\"P:\\Projects\\Olivares\\SPHY\\output\\QTOTSubBasinTSS.tss\"\n",
    "    \n",
    "    #Read and prepare obs streamflow\n",
    "    obs = pd.read_csv(obs_path)\n",
    "    obs = obs.iloc[:,1]\n",
    "\n",
    "    #Read and prepare sim streamflow\n",
    "    sim_read = pd.read_csv(sim_path,index_col=False,skiprows=4)\n",
    "    sim = []\n",
    "    \n",
    "    for v in range(len(sim_read)):\n",
    "        val = sim_read.iloc[v,0].strip()\n",
    "        val = \" \".join(val.split())\n",
    "        sim.append(val)\n",
    "\n",
    "    sim = pd.DataFrame(sim)\n",
    "    sim.columns = [\"name\"]\n",
    "    sim = sim[\"name\"].str.split(\" \",expand=True)\n",
    "    sim = pd.DataFrame(sim.iloc[:,2]).astype(float)\n",
    "    sim = sim.iloc[:,0]\n",
    "\n",
    "    #Calculate KGE'\n",
    "    kge = pd.DataFrame(he.evaluator(he.kgeprime, sim, obs),index=[\"KGE\",\"r\",\"γ\",\"β\"])\n",
    "    nash = pd.DataFrame(he.evaluator(he.nse, sim, obs),index=[\"NASH\"])\n",
    "    kge_r.append(kge[0][0])\n",
    "    corr.append(kge[0][1])\n",
    "    gamma.append(kge[0][2])\n",
    "    betha.append(kge[0][3])\n",
    "    nash_r.append(nash[0][0])\n",
    "    print(kge)\n",
    "    print(nash)\n",
    "    \n",
    "def soilpars(param,value):\n",
    "    #Read config_file\n",
    "    config_path = r\"P:\\Projects\\Olivares\\SPHY\\cfg\\sw\\sphy_config_sw_SA.cfg\" \n",
    "    config_oliv = ConfigParser()\n",
    "    config_oliv.read(config_path)\n",
    "    print(param)\n",
    "    #1st step 0.5*param\n",
    "    config_oliv[\"SOILPARS\"][str(param)] = str(value*0.5)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "    run_sphy()\n",
    "    performance(kling_gupta,coef_corr,variability,bias,nash_sutt)\n",
    "    #2nd step 1.5*param\n",
    "    config_oliv[\"SOILPARS\"][str(param)] = str(value*1.5)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "    run_sphy()\n",
    "    performance(kling_gupta,coef_corr,variability,bias,nash_sutt)        \n",
    "    #3rd original value\n",
    "    config_oliv[\"SOILPARS\"][str(param)] = str(value)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "\n",
    "def groundwpars(param,value):\n",
    "    #Read config_file\n",
    "    config_path = r\"P:\\Projects\\Olivares\\SPHY\\cfg\\sw\\sphy_config_sw_SA.cfg\" \n",
    "    config_oliv = ConfigParser()\n",
    "    config_oliv.read(config_path)\n",
    "    print(param)\n",
    "    #1st step 0.5*param\n",
    "    config_oliv[\"GROUNDW_PARS\"][str(param)] = str(value*0.5)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "    run_sphy()\n",
    "    performance(kling_gupta,coef_corr,variability,bias,nash_sutt)\n",
    "    #2nd step 1.5*param\n",
    "    config_oliv[\"GROUNDW_PARS\"][str(param)] = str(value*1.5)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "    run_sphy()\n",
    "    performance(kling_gupta,coef_corr,variability,bias,nash_sutt)        \n",
    "    #3rd original value\n",
    "    config_oliv[\"GROUNDW_PARS\"][str(param)] = str(value)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "        \n",
    "def glacier(param,value):\n",
    "    #Read config_file\n",
    "    config_path = r\"P:\\Projects\\Olivares\\SPHY\\cfg\\sw\\sphy_config_sw_SA.cfg\" \n",
    "    config_oliv = ConfigParser()\n",
    "    config_oliv.read(config_path)\n",
    "    print(param)\n",
    "    #1st step 0.5*param\n",
    "    config_oliv[\"GLACIER\"][str(param)] = str(value*0.5)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "    run_sphy()\n",
    "    performance(kling_gupta,coef_corr,variability,bias,nash_sutt)\n",
    "    #2nd step 1.5*param\n",
    "    config_oliv[\"GLACIER\"][str(param)] = str(value*1.5)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "    run_sphy()\n",
    "    performance(kling_gupta,coef_corr,variability,bias,nash_sutt)        \n",
    "    #3rd original value\n",
    "    config_oliv[\"GLACIER\"][str(param)] = str(value)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "\n",
    "def snow(param,value):\n",
    "    #Read config_file\n",
    "    config_path = r\"P:\\Projects\\Olivares\\SPHY\\cfg\\sw\\sphy_config_sw_SA.cfg\" \n",
    "    config_oliv = ConfigParser()\n",
    "    config_oliv.read(config_path)\n",
    "    print(param)\n",
    "    #1st step 0.5*param\n",
    "    config_oliv[\"SNOW\"][str(param)] = str(value*0.5)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "    run_sphy()\n",
    "    performance(kling_gupta,coef_corr,variability,bias,nash_sutt)\n",
    "    #2nd step 1.5*param\n",
    "    config_oliv[\"SNOW\"][str(param)] = str(value*1.5)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "    run_sphy()\n",
    "    performance(kling_gupta,coef_corr,variability,bias,nash_sutt)        \n",
    "    #3rd original value\n",
    "    config_oliv[\"SNOW\"][str(param)] = str(value)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "        \n",
    "def climate(param,value):\n",
    "    #Read config_file\n",
    "    config_path = r\"P:\\Projects\\Olivares\\SPHY\\cfg\\sw\\sphy_config_sw_SA.cfg\" \n",
    "    config_oliv = ConfigParser()\n",
    "    config_oliv.read(config_path)\n",
    "    print(param)\n",
    "    #1st step 0.5*param\n",
    "    config_oliv[\"CLIMATE\"][str(param)] = str(value*0.5)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "    run_sphy()\n",
    "    performance(kling_gupta,coef_corr,variability,bias,nash_sutt)\n",
    "    #2nd step 1.5*param\n",
    "    config_oliv[\"CLIMATE\"][str(param)] = str(value*1.5)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "    run_sphy()\n",
    "    performance(kling_gupta,coef_corr,variability,bias,nash_sutt)        \n",
    "    #3rd original value\n",
    "    config_oliv[\"CLIMATE\"][str(param)] = str(value)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "        \n",
    "def routing(param,value):\n",
    "    #Read config_file\n",
    "    config_path = r\"P:\\Projects\\Olivares\\SPHY\\cfg\\sw\\sphy_config_sw_SA.cfg\" \n",
    "    config_oliv = ConfigParser()\n",
    "    config_oliv.read(config_path)\n",
    "    print(param)\n",
    "    #1st step 0.5*param\n",
    "    config_oliv[\"ROUTING\"][str(param)] = str(value*0.5)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "    run_sphy()\n",
    "    performance(kling_gupta,coef_corr,variability,bias,nash_sutt)\n",
    "    #2nd step 1.5*param\n",
    "    config_oliv[\"ROUTING\"][str(param)] = str(value*1.5)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)\n",
    "    run_sphy()\n",
    "    performance(kling_gupta,coef_corr,variability,bias,nash_sutt)        \n",
    "    #3rd original value\n",
    "    config_oliv[\"ROUTING\"][str(param)] = str(value)\n",
    "    with open(config_path, 'w') as conf: #Save the config file with the params edited\n",
    "        config_oliv.write(conf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51b01df8-d357-48fe-b7ea-71b36e81d50d",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2.- Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87289ad4-a619-45be-9763-b3e17c62485c",
   "metadata": {},
   "outputs": [],
   "source": [
    "kling_gupta = []\n",
    "coef_corr = []\n",
    "variability = []\n",
    "bias = []\n",
    "nash_sutt = []\n",
    "\n",
    "#Base run\n",
    "run_sphy()\n",
    "performance(kling_gupta,coef_corr,variability,bias,nash_sutt)\n",
    "\n",
    "#Vector with original values\n",
    "\n",
    "param_names = [\"RootDepthFlat\",\n",
    "              \"SubDepthFlat\",\n",
    "              \"CapRiseMax\",\n",
    "              \"GwDepth\",\n",
    "              \"GwSat\",\n",
    "              \"deltaGw\",\n",
    "              \"BaseThresh\",\n",
    "              \"alphaGw\",\n",
    "              \"YieldGw\",\n",
    "              \"DDFG\",\n",
    "              \"DDFDG\",\n",
    "              \"Debris_factor\",\n",
    "              \"GlacF\",\n",
    "              \"Tcrit\",\n",
    "              \"SnowSC\",\n",
    "              \"DDFS\",\n",
    "              \"r\",\n",
    "              \"Swfactor\",\n",
    "              \"kx\"]\n",
    "\n",
    "param_values = [300,\n",
    "               915,\n",
    "               5,\n",
    "               4000,\n",
    "               4000,\n",
    "               3,\n",
    "               0.5,\n",
    "               0.5,\n",
    "               0.05,\n",
    "               21.14,\n",
    "               5.8,\n",
    "               0.2,\n",
    "               0.5,\n",
    "               12,\n",
    "               0.5,\n",
    "               12,\n",
    "               0.1,\n",
    "               0.6,\n",
    "               0.5]\n",
    "\n",
    "for p in range(len(param_names)):\n",
    "    if p == 0 or p == 1 or p == 2:\n",
    "        soilpars(param_names[p],param_values[p])\n",
    "    if p == 3 or p == 4 or p == 5 or p == 6 or p == 7 or p == 8:\n",
    "        groundwpars(param_names[p],param_values[p])\n",
    "    if p == 9 or p == 10 or p == 11 or p == 12:\n",
    "        glacier(param_names[p],param_values[p])\n",
    "    if p == 13 or p == 14 or p == 15 or p == 16:\n",
    "        snow(param_names[p],param_values[p])\n",
    "    if p == 17:\n",
    "        climate(param_names[p],param_values[p])\n",
    "    if p == 18:\n",
    "        routing(param_names[p],param_values[p])\n",
    "\n",
    "#Array to DataFrame\n",
    "kling_gupta = pd.DataFrame(kling_gupta)\n",
    "coef_corr = pd.DataFrame(coef_corr)\n",
    "variability = pd.DataFrame(variability)\n",
    "bias = pd.DataFrame(bias)\n",
    "nash_sutt = pd.DataFrame(nash_sutt)\n",
    "\n",
    "#DataFrame to excel\n",
    "kling_gupta.to_excel(r\"P:\\Projects\\Olivares\\Calibration\\Sensitivity_analysis\\kge.xlsx\", index = False)\n",
    "coef_corr.to_excel(r\"P:\\Projects\\Olivares\\Calibration\\Sensitivity_analysis\\correlation.xlsx\", index = False)\n",
    "variability.to_excel(r\"P:\\Projects\\Olivares\\Calibration\\Sensitivity_analysis\\variability.xlsx\", index = False)\n",
    "bias.to_excel(r\"P:\\Projects\\Olivares\\Calibration\\Sensitivity_analysis\\bias.xlsx\", index = False)\n",
    "nash_sutt.to_excel(r\"P:\\Projects\\Olivares\\Calibration\\Sensitivity_analysis\\nash.xlsx\", index = False)\n",
    "\n",
    "end_time = datetime.now()\n",
    "print('Duration: {}'.format(end_time - start_time))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}