Assignment_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CEE 6200 Assignment 2\n",
    "Understanding hydrologic processes is a key to water resources systems managment. This assignment will walk you through some foundational concepts of hydrology and hydrologic modeling.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from find_droughts import find_droughts\n",
    "from scipy.stats import genextreme as gev\n",
    "from Hymod_functions import hymod_main"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Visualizing Hydrologic Data\n",
    "### (1a)\n",
    "Streamflow and precipitation data are usually plotted together as a combined hydrograph and hyetograph (from Greek hyetos, “rain”). Streamflow is plotted as a time series, while rainfall is shown as an inverted bar plot along the top of the graph. Streamflow labels are shown on the left y-axis, while rainfall labels are shown on the right y-axis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a combined hydrograph/hyetograph for this one year of Leaf River data. It should look like this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](Hydro-hyetograph_sol.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a first step, we'll load the data as a pandas data frame. Pandas is a helpful library for managing data in python: https://pandas.pydata.org/pandas-docs/stable/index.html. \n",
    "\n",
    "We'll the data using the function \"read_csv\", the data file is a csv file, so the \"separator\" is a ','. We'll then extract the first year of data using the \"iloc\" function (extract the first 365 days which correspond to rows 0 to 364)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the data\n",
    "leaf_data = pd.read_csv('LeafCatch.csv', sep=',')\n",
    "\n",
    "# extract the first year of data\n",
    "leaf_data = leaf_data.iloc[0:364]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can access a single column like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "precip = leaf_data['Precip']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print it to see what it looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(precip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot, we'll set it up for you. Use the basic plot function to create the hydrograph (https://matplotlib.org/api/_as_gen/matplotlib.pyplot.plot.html) and the bar function to create the hyetograph (https://matplotlib.org/api/_as_gen/matplotlib.pyplot.bar.html) **(5 points)**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first we'll make an axis for the hydrograph\n",
    "fig, strmflw_ax = plt.subplots(figsize=[12,6])\n",
    "strmflw_ax.set_ylim([0, 40])\n",
    "\n",
    "# then we'll make a second y-axis for the hyetograph\n",
    "precip_ax = strmflw_ax.twinx()\n",
    "precip_ax.set_ylim([0, 200])\n",
    "precip_ax.invert_yaxis()\n",
    "\n",
    "# to plot streamflow use strmflow_ax.plot()\n",
    "# to plot precip use precip_ax.bar()\n",
    "\n",
    "############# add your code here ###################\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "####################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (1b)\n",
    "Next, let's take a look at how the magnitude of streamflows vary. To visualize the variability of streamflow in a particular watershed, hydrologists often use a Flow Duration Curve. This is essentially a cumulative distribution function (CDF) that, for historical reasons, is flipped sideways. Here is an example of the flow duration curve that you’ll make:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](flow_duration_curve.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Streamflow is shown on the y-axis using a log scale, to emphasize differences in extreme values. The x-axis shows the fraction of time that a given flow value is exceeded. For example, the highest recorded flow value is exceeded 0% of the time, while the smallest recorded flow value is exceeded 100% of the time (at least, as far as historical data can tell us). \n",
    "Reproduce this plot for the same Leaf River data we used in Problem 1a. There are functions in python to create empirical CDFs, but it’s fairly straightforward to do yourself. Think about it this way: you already have the streamflow data (y-axis), and if you sort the streamflow values in order, you know exactly what percentage of days are greater or less than each flow value. **(5 points)**.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To set the y-axis to log scale, use the matplotlib command 'semilogy': https://kite.com/python/docs/matplotlib.pyplot.semilogy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "############# add your code here ###################\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "####################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Quantifying Extremes\n",
    "\n",
    "Water resources analysts are often interested in the extremes of data sets (floods/droughts etc.); in Part 2, we will examine two metrics for quantifying extremes. To complete this part of the assignment download the you will need the Leaf River data set from Part 1 as well as the Falls-Lake-Log-Streamflows.csv data set and the findDroughts function. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (2a) \n",
    "Return periods are a commonly used metric for quantifying the frequency of extreme events. The return period (RP) can be thought of as the expected waiting time between extreme events. The return period of an event with magnitude x is defined as:\n",
    "\n",
    "$$RP = \\frac{1}{1-F(x)}$$\n",
    "\n",
    "Where $F(x)$ is the cumulative distribution function of all events X. $1-F(x)$ is known as the exceedance probability. A ten year 24-hour rainfall event, for example, will have a rainfall depth that will only be exceeded on average once every ten years, and therefore its exceedance probability is: $1-0.90=0.10$\n",
    "\n",
    "$$10 years = \\frac{1}{1-0.90}$$\n",
    "\n",
    "To solve for the magnitude of the 10-year storm, one would solve for x in equation 1 by finding the inverse of CDF of the corresponding return period. The magnitude of the 10-year 24-hour rainfall event from the Leaf Catchment data set can be calculated as follows:\n",
    "\n",
    "$$10 years = \\frac{1}{1-F(x)} \\rightarrow F(x)=0.90 \\rightarrow x = F(x)^{-1} = 126.59mm$$\n",
    "\n",
    "Where $F(x)^{-1}$ is the inverse cdf of the GEV fit to the historical data set.\n",
    "\n",
    "In this exercise, we will calculate the magnitude, x, of the 20, 50, 100, 200 and 500-year return period 24 hour rainfall events based of the most current 30 years of the Leaf Catchment historical data set.  In Matlab or your language of choice, use the following procedure to find the magnitude of each return period **(10 points)**:\n",
    "\n",
    "\n",
    "1.\tReload the leaf river data and extract the most current 30 years of precipitation. You can assume all years have 365 days (i.e. ignore leap years). \n",
    "2.\tFind the annual maximum value for each of the 30 years.\n",
    "3.\tFit a GEV to the annual maximum data set. In python, use the command gev.fit(), we've already loaded this from the scipy library: https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.stats.genextreme.html\n",
    "4.\tCalculate the magnitude of each return period. A helpful command to use is .ppf from the GEV distribution made above\n",
    "5.\tSummarize your findings in the markdown table below the code (just double click to edit)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "############# add your code here ###################\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "####################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "---\n",
    "\n",
    "\n",
    "*...Edit the table below-->(double click to edit)...*\n",
    "\n",
    "| RP            | Magnitude (mm)                     |\n",
    "| ------------- |:----------------------------------:|\n",
    "| 10 year       | 126.59 (good sanity check)         |\n",
    "| 20 years      | (edit)                             |\n",
    "| 50 year       | (edit)                             |\n",
    "| 100 year      | (edit)                             |\n",
    "| 200 year      | (edit)                             |\n",
    "| 500 year      | (edit)                             |\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (2b) \n",
    "Defining what constitutes a period of drought is a difficult question that has been a source of contention in the water resources community for decades. One method of defining drought is through a six month Standardized Streamflow Indicator (SSI6), also known as SRI, Standardized Runoff Indicator (Mckee et. al 1993). The SSI6 is a measure of the deviation from the historical mean of the 6 month moving average of a given streamflow record. Using this statistic, a drougth can be defined as a time interval where the SSI6 < 0 for at least 6 consecutive months and SSI6 < -1 for at least one week within that period (Herman et al 2016). In this exercise, we will use this criterion to search through a historical data set, identify periods of drought and visualize them over the historical record.  Streamflows in the data set of interest are known to follow a log-normal distribution. Log(streamflow) data has been uploaded to blackboard. \n",
    "\n",
    "To find the periods of historical drought, follow the procedure below **(5 points)**:\n",
    "\n",
    "\n",
    "1.\tLoad the log-transformed streamflow data from the Falls-Lake-Log-Streamflows.csv file on blackboard. This file contains a time series of $Y_t=ln⁡(Q_t )$, where Qt is the observed inflows and Yt is the log-transformed stream flows.\n",
    "2.\tNormalize the weekly log(streamflow) by subtracting the mean of the entire record and dividing by the standard deviation:\n",
    "$Z_t=(Y_t-μ ̂)/\\hat{\\sigma} ̂$\n",
    "This value represents how many standard deviations away from the sample mean a given week’s inflow was.\n",
    "3.\tCalculate the SSI6, the 6-month moving average of Zt for each week. Save these values in a new pandas dataframe with a column named \"ssi\" (https://kite.com/python/answers/how-to-create-pandas-dataframe-from-a-numpy-array-in-python).\n",
    "4.\tUse the findDroughts function to search through the historical record and identify the periods of drought, defined as at least 6 months (24 consecutive weeks) where the SSI6 < 0 and for at least 1 week $SSI6 < -1$. Make sure to examine the function to until you understand how it works. The output of this function will be a new Pandas dataframe with three columns \"start\" containing the start week of the drought, \"end\" containing the end weeks of each drought and \"severity\" containing the minimum SSI6 value during the drought.  \n",
    "5.\tPlot the time series of log(inflows) and shade the regions that have been identified as droughts. The final plot should resemble the one shown below. \n",
    "\n",
    "!['example'](fallsLakeDroughts.png)\n",
    "\n",
    "Hint: to make this plot, use the python function \"axvspan\" to shade the droughts, make sure to adjust the transparency using the alpha argument (https://matplotlib.org/gallery/subplots_axes_and_figures/axhspan_demo.html#sphx-glr-gallery-subplots-axes-and-figures-axhspan-demo-py)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the data (we're doing this for you because of the date format)\n",
    "df = pd.read_csv('Falls-Lake-Log-Streamflows.csv', parse_dates=[0], index_col=0)\n",
    "\n",
    "# calculate SSI6 (hint: Pandas function mean, std and rolling will be helpful)\n",
    "\n",
    "############# add your code here ###################\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "####################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add plotting code below "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "############# add your code here ###################\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "####################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 3: Rainfall-Runoff Modeling\n",
    "The objective of Part 3 is to familiarize you with conceptual rainfall-runoff modeling, including explorations of parameter uncertainty and basic sensitivity analysis. This directory contains Python code for the HyMod rainfall-runoff model. Take some time to review both the model diagram below, and the code that makes it work in the Hymod_functions folder. Then answer the following questions:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](hymod.png)\n",
    "\n",
    "---\n",
    "\n",
    "### Model Parameters\n",
    "**Huz**: Maximum storage size\n",
    "\n",
    "**b**:      Distribution of storage sizes between 0 and Huz\n",
    "\n",
    "**Alp**:    Fraction of runoff contributing to quick flow \n",
    "\n",
    "**Kq**:     Quick flow residence time of linear infinite reservoir (the Kq values of all three linear reservoirs are the same)\n",
    "\n",
    "**Ks**:Slow flow residence time of linear infinite reservoir\n",
    "\n",
    "### Variables\n",
    "\n",
    "**PP**:    Precipitation\n",
    "\n",
    "**ET**:    Evapotranspiration\n",
    "\n",
    "**OV**:    Runoff\n",
    "\n",
    "**Qq**:    Quick Flow\n",
    "\n",
    "**Qs**:    Slow Flow\n",
    "\n",
    "**QQ**:    Streamflow (Quick Flow + Slow Flow)\n",
    "\n",
    "**XHuz** and **XCuz**: Current moisture state of soil moisture accounting component (as depth XH or volume XC)\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (3a)\n",
    "In the context of the diagram above, what does the Nash function do? Enter your answer in the cell below **(2 points)**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...*Add markdown explanation here-->(double click to edit)...*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (3b) \n",
    "In the context of the diagram above, what does the Pdm01 function do? **(2 points)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...*Add markdown explanation here-->(double click to edit)...*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (3c)\n",
    "Use Hymod to perform single model runs, using the same Leaf River data as in Part 1. There are five parameters that must be chosen by the user: Kq, Ks, Alp, Huz, and B. Alter these parameters manually until you observe an acceptable match between the observed and simulated streamflow. (This can be done visually by plotting LeafCatch['Strmflw'].iloc[0:365] and Model['Q'] vs. time). Report the final parameter values you used. Provide the plot comparing the observed and simulated streamflows on a hydrograph (no need to include the hyetograph this time). Also provide a plot comparing the observed and simulated flow duration curves. How well do they match up? **(11 points)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### assign your parameters here #####\n",
    "\n",
    "Nq =    # number of quickflow routing tanks\n",
    "Kq =     # quickflow routing tanks parameters \t\t\t\t- Range [0.1, 1]\n",
    "Ks =     # slowflow routing tanks rate parameter \t\t\t- Range [0, 0.1]\n",
    "Alp =    # Quick-slow split parameters \t\t\t\t\t\t- Range [0, 1]\n",
    "Huz =    # Max height of soil moisture accounting tanks \t- Range [0, 500]\n",
    "B =     # Distribution function shape parameter \t\t\t\t- Range [0, 2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the Model in the following cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Model = hymod_main.main(Nq, Kq, Ks, Alp, Huz, B, 'LeafCatch.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the results below, you'll likely have to do this several times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "############# add your code here ###################\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "####################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (3d)\n",
    "Create your own function(s) to compute two error metrics: the root-mean squared error (RMSE) and runoff coefficient error (ROCE, which describes the long-term average ratio of streamflow to precipitation). Use these to quantitatively compare the observed and simulated time series you found in Part 3c. How well does your “best” parameter set perform?  **(5 points)**\n",
    "\n",
    "$$RMSE = \\sqrt{\\frac{1}{N}\\sum_{t=1}{N}(Q_{sim,t}-Q_{obs,t})^2}$$\n",
    "\n",
    "$$ROCE = \\left\\lvert\\frac{\\bar{Q}_{sim}-\\bar{Q}_{obs}}{\\bar{P}}\\right\\rvert$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "############# add your code here ###################\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "####################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...*Add markdown explanation here-->(double click to edit)...*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (3e) \n",
    "Perform a Monte Carlo analysis of the feasible parameter space using uniform random sampling (assuming the model parameters are independent). The numpy library’s rand function can help you with this. For example, to create 20 samples of the parameter Kq, using its lower and upper bounds, you could do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "Kq_max = 1 # upper bound\n",
    "Kq_min = 0.1 # lower bound\n",
    "N = 20 # number of samples\n",
    "Kq_samples = (Kq_max-Kq_min)*np.random.rand(N)+Kq_min"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take all N parameter sets that you generate and run them through the model. The value of N is up to you, depending on your computing resources, but you should be able to run at least 1,000 evaluations without issue. Use the functions you wrote in (3d) to calculate the RMSE and ROCE values for each model evaluation. When these model runs are complete, please do the following:\n",
    "\n",
    "i.\tYou can think of this uniform random sample as a crude form of optimization. Find the parameter set with the best RMSE value and report its parameters. Then, find the parameter set with the best ROCE value and report its parameters. Are these two parameter sets approximately the same? Why or why not? **(5 points)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "############# add your code here ###################\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "####################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...*Add markdown explanation here-->(double click to edit)...*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ii.\tCreate scatter plots of each parameter vs. the RMSE and ROCE metrics. Since you have 5 parameters and 2 error metrics, you should have 10 total plots (or subplots). Please put the error metric on the y-axis, and the parameter value on the x-axis. Which parameters appear to be most responsible for changes in the RMSE value? What about the ROCE value? How do you know? This is a simple form of sensitivity analysis, a topic we’ll return to later in the class. **(5 points)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "############# add your code here ###################\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "####################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...*Add markdown explanation here-->(double click to edit)...*"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}