python_essentials_III_session_03_colab.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "# Session 3\n",
    "\n",
    "During this final session we will look at some practical hydrological applications of the various Python libraries. There is only so much we can cover during a single session. A great source of information is <A href=\"https://agu-h3s.org/2021/03/29/resources-for-programming-in-hydrology/\">this overview of Python resources for programming in hydrology</A> created by Jerom Aerts. You will find information there that we also already covered during this course, but also many more pointers to resources that you might find useful in your own day-to-day work. \n",
    "\n",
    "## Time series, graphs and spatial data\n",
    "Before doing anything, we must import the libraries we'll be using"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "!wget https://raw.githubusercontent.com/vincentpost/PythonEssentialsIII/main/09246500_time_series_data.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time series statistics\n",
    "\n",
    "### Stream discharge\n",
    "\n",
    "We ended the previous session with fitting a power equation to the stage-discharge data for Elkhead creek. We'll continue by converting the measured stream stages to discharge values. From the time series obtained that way, we can create a flow duration curve, that is, a cumulative density plot of the probability (in percent of time of a record) that the discharge equals or exceeds a certain value.\n",
    "\n",
    "Like we learned last session, we can load the data and parse the dates in the csv file using `Pandas`' `read_csv` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# USGS 09246500 ELKHEAD CREEK NEAR CRAIG, CO\n",
    "# https://waterdata.usgs.gov/nwis/measurements/?site_no=09246500\n",
    "df = pd.read_csv('09246500_time_series_data.txt', \n",
    "                 sep='\\t',\n",
    "                 index_col=2,\n",
    "                 parse_dates=True,\n",
    "                 skiprows=[i for i in range(32) if i != 30])\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file header explains that the column named '279703_00065' contains the gauge heights, and the column '212023_00060' has the discharge (in cubic feet per second, CFS). Printing the first five rows of data shows there are lots of NaN values. Closer inspection of the file (in a text editor or a spreadsheet program) reveals that most of the stage data fields are empty. The year 2021 has a nearly complete set of measurements though, so let's select that for the remainder of the exercise. This is easily done using the `year` attribute of the date/time field of the index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df = df.loc[df.index.year == 2021]\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's quickly plot the time series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(df.index, df['279703_00065'], '-')\n",
    "plt.ylabel('Stage [feet]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the rating curve from session 2 ($Q = 0.547  h ^ {5.15}$) we can convert the stream stage to a discharge. We'll add it to the dataframe `df` and give it column label Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# From session 2:\n",
    "b, a = 5.1528685333467035, 0.5467695289254147\n",
    "df['Q'] = a * df['279703_00065'] ** b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the calculated discharge and compare it to the reported discharge from the downloaded file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(df.index, df['Q'], '-', label='Session 2 rating curve')\n",
    "plt.plot(df.index, df['212023_00060'], '-', label='USGS reported value')\n",
    "plt.legend()\n",
    "plt.ylabel('Discharge [CFS]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The curves don't exactly line up, but since this is an exercise only, let's continue with our calculated discharge data. If this were a real investigation, we'd be scrutinising the data a lot more of course and try to understand the causes for these differences. Note for example that the file hardly contains any discharge figures for the winter season (even though the stream stage was measured), presumably because the stream was frozen. The results of what follows therefore make no hydrological sense, but the purpose is not to analyse this stream, the purpose is to demonstrate the usage of `Pandas` to calculate the probabilities of exceedance. \n",
    "\n",
    "The first step in the process is to resample the 15-minute data to daily values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dfd = df['Q'].resample('1D').mean()\n",
    "dfd.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because there is only one column of data, `Pandas` actualy just created a `Series` object instead of a `DataFrame`. Although they are similar, they are different at the same time. For example, both a `Series` and a `DataFrame` have an index, unlike a `DataFrame`, a `Series` does not have column names. Also, we can't easily add a column to `dfd` like we did in the above example for the discharge Q. To overcome this, we can simply use `to_frame` to convert the `Series` back into a `DataFrame`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dfd = dfd.to_frame()\n",
    "dfd.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `rank` method in `Pandas` computes the numerical data ranks. If there are duplicate values then equal values are assigned a rank that is the average of the ranks of those values. Using `pct=True` converts the ranks into percentiles. For plotting purposes, the data need to be sorted, which is done using the `sort_values` function. The first argument determines which column is used for sorting, in this case it is 'cdf'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dfd['cdf'] = dfd.rank(method='average', pct=True)\n",
    "dfd = dfd.sort_values('cdf', ascending=False)\n",
    "dfd.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the cumulative density function is easy. For flow duration curves it is common to show the discharge values in order of decreasing probability, which is accomplished by using the `xlim` function to reverse the min and max values of the $x$ axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.semilogy(dfd['cdf'] * 100, dfd['Q'])\n",
    "plt.xlim(100, 0)\n",
    "plt.xlabel('Exceedance probability [%]')\n",
    "plt.ylabel('Discharge [CFS]')\n",
    "plt.grid(True, which=\"both\", ls=':', color='0.65')\n",
    "#plt.savefig('cdf.png', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that there are different ways to calculate a flow duration curve. This is just one illustrative example. More sophisticated methods can also be implemented, but that is beyond an essentials course.\n",
    "\n",
    "### Sydney rainfall\n",
    "\n",
    "The code cell below loads the data from a comma separated values file that contains the daily rainfall amounts for the Randwick BOM station in  Sydney. We will use the data to calculate the yearly totals and plot the probality of exceedance graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# http://www.bom.gov.au/jsp/ncc/cdio/weatherData/av?p_nccObsCode=136&p_display_type=dailyDataFile&p_stn_num=66037\n",
    "!wget https://raw.githubusercontent.com/vincentpost/PythonEssentialsIII/main/IDCJAC0009_66037_1800_Data.csv\n",
    "    \n",
    "df = pd.read_csv('IDCJAC0009_66037_1800_Data.csv',\n",
    "                 parse_dates={'dt': [2, 3, 4]})\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the `parse_dates` argument is used to make sure that the year, month and day values, which are all stored in a separated column, are joined together and converted to a date format. The only thing is that the time is set to midnight, while rainfall amouns are read at 9am every morning (see also the file header). Not to worry, we'll simply add 9 hours to the timestamp using the `TimeDelta` function. Once we've done that, we'll make the 'dt' column the index for the DataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df['dt'] = df['dt'].add(pd.Timedelta(9, unit='h'))\n",
    "df.set_index('dt', inplace=True, drop=True)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The column with the rainfall amounts has a really long name, so to avoid having to type it all the time, we store it as a string variable, which we'll call `cname`. We'll then get rid of all the other columns by selecting only the column with the measured rainfall amounts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cname = df.columns[2]\n",
    "df = df[cname]\n",
    "print(df.tail(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating the yearly sums is easy with the `resample` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dfy = df.resample('1Y').sum()\n",
    "dfy.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A bar graph gives a good first impression of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.bar(dfy.index, dfy, width=365)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's calculate the 1961 - 1990 average of the annual rainfall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "idx = (dfy.index.year >= 1961) & (dfy.index.year <= 1990)\n",
    "rmean = dfy.loc[idx].mean()\n",
    "print(f'Mean annual rainfall (1961 - 1990): {rmean:0.1f} mm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The value can also be displayed as a line in the plot above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "x0 = dfy.index[0]\n",
    "x1 = dfy.index[-1]\n",
    "plt.figure()\n",
    "plt.bar(dfy.index, dfy.array, width=365)\n",
    "plt.plot([x0, x1], [rmean, rmean], 'r--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once more we can use the `rank` method in `Pandas` to plot the cumulative probabilty function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "dfy = dfy.to_frame()\n",
    "dfy['cdf'] = dfy.rank(method='average', pct=True)\n",
    "dfy = dfy.sort_values('cdf')\n",
    "dfy.head()\n",
    "plt.figure()\n",
    "plt.plot(dfy[cname], dfy['cdf'] * 100)\n",
    "plt.xlabel('Annual rainfall')\n",
    "plt.ylabel('Probability of exceedance [%]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise 3.1***: Use the `Pandas` functions `min`, `max`, `median` on the 1961 - 1990 portion of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Type your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise 3.2***: Use the `count` function to determine the number of years in which the rainfall exceeded 2000 mm (use the entire dataset, not just the 1961 - 1990 period). Check your answer by comparing it to the bar chart above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Type your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More on graphs\n",
    "\n",
    "You may have seen them: The warming stripes that show how annual temperatures have increased over the past few decades. They were first created by  Ed Hawkings (see <A href=\"https://www.climate-lab-book.ac.uk/2018/warming-stripes/\">Climate Lab Book</A>) and were adopted as the logo of the <A href=\"https://scientistsforfuture.org/\">Scientists for Future</A> movement. The `Matplotlib` homepage contains an <A href=\"https://matplotlib.org/matplotblog/posts/warming-stripes/\">explanation</A> on how to create them, but it is a rather complex example and there is in fact a much easier way, which is included as a bonus exercise at the end of this notebook. \n",
    "\n",
    "For the next exercise we will be using the HadCRUT5 near surface temperature data set to demonstrate some of `Matplotlib`'s capabilities. The aim is to plot three graphs in a single figure and add an uncertainty band around the data. The HadCRUT5 data set is created by the British Met Office Hadley Centre. The files for this exercise were downloaded on 19 April 2022 from <A href=\"https://www.metoffice.gov.uk/hadobs/hadcrut5/data/current/download.html\">https://www.metoffice.gov.uk/hadobs/hadcrut5/data/current/download.html</A>. They contain the annual mean temperature anomalies for the northern- and southern hemisphere, as well as the entire globe since 1850. The temperature anomaly is the difference of the annual temperature and the mean annual temperature for the reference period 1961 - 1990.\n",
    "\n",
    "Let's have a look at the data. First, we'll store the filenames in a `list`, so that we can automate the creation of the graphs with a `for` loop later on"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!wget https://raw.githubusercontent.com/vincentpost/PythonEssentialsIII/main/HadCRUT.5.0.1.0.analysis.summary_series.global.annual.csv\n",
    "!wget https://raw.githubusercontent.com/vincentpost/PythonEssentialsIII/main/HadCRUT.5.0.1.0.analysis.summary_series.northern_hemisphere.annual.csv\n",
    "!wget https://raw.githubusercontent.com/vincentpost/PythonEssentialsIII/main/HadCRUT.5.0.1.0.analysis.summary_series.southern_hemisphere.annual.csv\n",
    "\n",
    "# https://www.metoffice.gov.uk/hadobs/hadcrut5/data/current/download.html\n",
    "fnames = ['HadCRUT.5.0.1.0.analysis.summary_series.global.annual.csv',\n",
    "          'HadCRUT.5.0.1.0.analysis.summary_series.northern_hemisphere.annual.csv',\n",
    "          'HadCRUT.5.0.1.0.analysis.summary_series.southern_hemisphere.annual.csv',\n",
    "          ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, later on we'll want a title above each graph, so let's store those in a `list` as well (making sure that the order of the titles is consistent with the order in `fnames`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "titles = ['Global', 'Northern hemisphere', 'Southern hemisphere']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As there are three sets of data, we'll need three individual graphs. For this, we'll use `Matplotlib`'s `subplots` function.\n",
    "\n",
    "The `for` loop executes a number of steps. First, it uses `Pandas`' `read_csv` function to read the data files. It then uses `Matplotlib`'s `plot`function to plot the temperature anomaly as a function of time. Because the data also contain information about the uncertainty, we can plot an envelope around the line to indicate the uncertainty range. Note that the envelope is plotted in black (hence the 'k' for the color) but with a transparency `alpha` of 20% to make it appear grey. Each graph is also give a label on the y axis and a title. The final command (`plt.tight_layout()`) optimizes the layout of the figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(nrows=3, figsize=(9,6))\n",
    "\n",
    "for fn, title, ax in zip(fnames, titles, axs):\n",
    "    df = pd.read_csv(fn, index_col=0)\n",
    "    \n",
    "    ax.plot(df.index, df['Anomaly (deg C)'])\n",
    "    ax.fill_between(df.index, df.iloc[:, 1], df.iloc[:, 2], color='k', alpha=0.2 )\n",
    "\n",
    "    ax.set_ylabel('Temp. anom. [\\u2103]')\n",
    "    ax.set_title(title)\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise 3.3***: Run the code cell above but comment out the `plt.tight_layout()` line. Observe what happens.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise 3.4***: The range of the horizontal axis can be set using `ax.set_xlim(xmin, xmax)`, in which `xmin` and `xmax` are the axis limits. Add this function to the code cell above the figure to plot only the data between 1920 and 2020."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Animation\n",
    "\n",
    "In the following example it will be demonstrated how the drawdown (i.e. the fall of the groundwater level) of an abstraction well as a function of time can be visualised in an animation. The first step is to define the function that calculates the drawdown of the well as a function of time $t$ and radial distance $r$\n",
    "\n",
    "$s = \\frac{Q}{4 \\pi T} W(u)$\n",
    "\n",
    "where $Q$ is the well abstraction rate, $T$ the aquifer transmissivity and $W(u)$ is the Theis well function, known as the exponential integral outside hydrogeology\n",
    "\n",
    "$\\int_u^\\infty \\frac{e^u}{u} du$\n",
    "\n",
    "in which \n",
    "\n",
    "$u = \\frac{r^2S}{4Tt}$\n",
    "\n",
    "where $S$ is the aquifer storativity.\n",
    "\n",
    "As luck would have it, the exponential integral is part of `scipy`'s special functions, so all we have to do is import it, and we can code up a function to calculate the drawdown according to the Theis equation. You'll notice that a function is defined by typing `def` followed by the function name (`theis` in this example). The arguments that can be passed to the function are enclosed between parentheses. The arguments `r` (radial distance from the well), `t` (time), `S` (aquifer storativity) and `T` (aquifer transmissivity) are compulsory. The argument `Q` (the well discharge) is a so-called keyword argument or kwarg, which is an argument that is assigned a default value (in this case 100) in the function definition. The user does not have to specify a value for `Q` unless they want to override the default value. The drawdown is stored in array `s` and in the final line of the function, this array is 'returned' to the part of the program that called the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from scipy.special import exp1\n",
    "\n",
    "def theis(r, t, S, T, Q=100): # Q is a kwargt (keyword argument), it is optional and has a default value\n",
    "    u = S * r ** 2 / (4.0 * T * t)\n",
    "    s = Q / (4 * np.pi * T) * exp1(u)\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the function defined, we can create a graph that shows the drawdown for a given set of aquifer parameters after 0.1, 1, and 10 days since the well started pumping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "S = .01 # Storativity\n",
    "T = 10 # Transmissivity\n",
    "r = np.linspace(0.1, 100, 200) # Radial distance from well\n",
    "for t in [0.1, 1, 10]: # Times\n",
    "    ax.plot(r, theis(r, t, S, T), label = f'time = {t} d')\n",
    "\n",
    "ax.set_xlabel('radial distance (m)')\n",
    "ax.set_ylabel('drawdown (m) ')\n",
    "ax.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This does not look as nice as it could. A different way of visualizing is to show a single line and plot the different times one after another. This can be done interactively in `matplotlib` using the `FuncAnimation` function. In the following example, an empty graph is created first. The animation is created by `FuncAnimation`, which requires information about the figure that contains the graph to be animated, a function that does the actual updating, and the times to be passed to the `update` function. An `init` function is optional, it can be used to create the intial appearance of the graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from matplotlib.animation import FuncAnimation\n",
    "r = np.linspace(0.1, 1000, 200)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ln, = ax.plot(r, np.zeros_like(r))\n",
    "t_str = ax.text(800, 8, f\"t = {0:0.1f} d\",\n",
    "                ha='center', va='center')\n",
    "\n",
    "def init():\n",
    "    ax.set_xlim(0, 1000)\n",
    "    ax.set_ylim(10, 0)\n",
    "    ax.grid(True, ls=':')\n",
    "\n",
    "def update(t):\n",
    "    ln.set_ydata(theis(r, t, S, T))\n",
    "    t_str.set_text(f\"t = {t:0.2f} d\")\n",
    "\n",
    "ani = FuncAnimation(fig, update, frames=[0.1, 1, 10], init_func=init)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It works, but there is room for improvement! Before continuing, hit the blue/white on-off button in the upper right corner of the graph above to stop interaction. This will stop the animation and avoid any interference with  the next example, which is actually an exercise.\n",
    "\n",
    "By playing around with the `frames` argument of the `FuncAnimation` function, a much smoother looking animation can be created. Here it is done using the `logspace` function in `numpy`. This can be used to create an array with time intervals that are small initially and get bigger with time. This is useful because drawdown near the well changes rapidly in the beginning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from matplotlib.animation import FuncAnimation\n",
    "r = np.linspace(0.1, 1000, 200)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ln, = ax.plot(r, np.zeros_like(r))\n",
    "t_str = ax.text(800, 8, f\"t = {0:0.1f} d\",\n",
    "                ha='center', va='center')\n",
    "\n",
    "def init():\n",
    "    ax.set_xlim(0, 1000)\n",
    "    ax.set_ylim(10, 0)\n",
    "    ax.grid(True, ls=':')\n",
    "\n",
    "def update(t):\n",
    "    ln.set_ydata(theis(r, t, S, T))\n",
    "    t_str.set_text(f\"t = {t:0.2f} d\")\n",
    "\n",
    "ani = FuncAnimation(fig, update, frames=np.logspace(-1, 2, 100), init_func=init,\n",
    "                    blit=True, interval=1, repeat=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise 3.5***: Try to make the animation run twice as slow. Also move the text that shows the time to a different position and increase the fontsize."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geopandas\n",
    "\n",
    "Before you can run the next code cell, you must install `geopandas`. The installation guide can be found <A href=\"https://geopandas.org/getting_started/install.html\">here</A>. Once `geopandas` is installed you can import it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import geopandas as gpd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the function `read_file` it is possible to open a shapefile. The result will be a `GeoDataFrame`, which is like a `pandas` DataFrame, but with the spatial attributes of the features in the shapefile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "gdf = gpd.read_file(\"shp/river_vertices.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the file is loaded, we can loop over all the river nodes, and get their coordinates. These are used to calculate the distance along the river course"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "xvalues = np.array([pt.x for pt in gdf[\"geometry\"]])\n",
    "yvalues = np.array([pt.y for pt in gdf[\"geometry\"]])\n",
    "dx = np.diff(xvalues, prepend=0)\n",
    "dy = np.diff(yvalues, prepend=0)\n",
    "l = np.sqrt(dx ** 2 + dy ** 2)\n",
    "lsum = np.cumsum(l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the distance of each node known, `numpy`'s `interp` function can be used to do the interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "l0 = lsum[0]\n",
    "l1 = lsum[-1]\n",
    "stage0 = gdf[\"stage\"].iloc[0]\n",
    "stage1 = gdf[\"stage\"].iloc[-1]\n",
    "\n",
    "lint = lsum[1:-1]\n",
    "\n",
    "stageint = np.interp(lint, [l0, l1], [stage0, stage1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the interpolated values are stored back into the `GeoDataFrame`, which is then saved as a shapefile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "gdf[\"stageint\"] = np.hstack((stage0, stageint, stage1))\n",
    "\n",
    "gdf.to_file(\"shp/river_vertices_int.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bonus exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discharge calculation using the salt dilution method\n",
    "\n",
    "The principe of the salt dilution method is that when a mass of salt is inserted to the flow of a stream, measuring the concentration change with time will allow calculating the stream discharge. When the volumetric flow rate of the stream (for example in l/s) is $Q$ and the concentration of the salt is $C$ (in g/l), the salt mass that passes at a point downstream of the injection point per unit of time (i.e., g/s) is $QC$. The total mass that passes that point has to be equal to the mass of salt injected $M$ (in g), which can be expressed mathematically as\n",
    "\n",
    "$$M = \\int QC dt$$\n",
    "\n",
    "which, when $Q$ remained constant during the measurement, is\n",
    "\n",
    "$$M = Q \\int C dt$$\n",
    "\n",
    "The term $\\int C dt$ is equal to the area below the curve in a salt concentration versus time plot, as shown in the figure below\n",
    "\n",
    "<IMG src=\"https://github.com/vincentpost/PythonEssentialsIII/blob/main/salt_diluton.png\" width=450>\n",
    "    \n",
    "So when $\\int C dt$ can be determined, $Q$ is calculated as $M / \\int C dt$. Back in the old days, the area under the curve was determined by plotting the data on milimeter paper and counting the squares. Nowadays, the integral can be determined with one line of Python code. Let's use some real data to show how it is done. The data are stored in a an Excel file called `stream_ec_data.xlsx`. You may want to inspect the contents of this file in a spreadsheet program  in order to better understand the arguments passed with the `read_excel` function in the code cell below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!wget https://github.com/vincentpost/PythonEssentialsIII/raw/main/stream_ec_data.xlsx \n",
    "\n",
    "df = pd.read_excel('stream_ec_data.xlsx',\n",
    "                   #index_col = [0, 1],\n",
    "                   parse_dates={'dt': [0, 1]},\n",
    "                   skiprows=[1])\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the file contains timestamps and measured electrical conductivity (or rather: specific conductance) values. For the integration procedure it is easiest to convert the timestamps to integer values. We'll convert them to seconds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "xvalues = (df['dt'] - df['dt'].values[0])\n",
    "xvalues = xvalues.dt.total_seconds().array\n",
    "print(xvalues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll also convert the electrical conductivities to salt (grams of NaCl) using a conversion factor. Note that the background conductivity is subtracted first. If the value is less than zero, we'll force it to be zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sc_bg = 2.495\n",
    "df['SpCond'] = df['SpCond'] - sc_bg\n",
    "idx = df['SpCond'] < 0\n",
    "df['SpCond'].loc[idx] = 0\n",
    "df['NaCl'] = df['SpCond'] * 0.6282"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we use the `trapz` function from `scipy.integrate` to determine the area under the curve. It is stored in the variable `V`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import trapz\n",
    "V = trapz(df['NaCl'], x=xvalues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The injected salt mass $M$ was 2000 grams for this example, so Q becomes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "M = 2000\n",
    "Q = M / V\n",
    "print(f'The stream discharge is {Q:0.1f} l/s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "## Climate stripes\n",
    "\n",
    "This next example will use the HeadCRUT5 data set that we visualized earlier to recreate the famouse climate stripes (see <A href=\"https://www.climate-lab-book.ac.uk/2018/warming-stripes/\">Climate Lab Book</A>). The graph will look slightly different, but definitively close enough. You will notice how compact the code is, demonstrating the power of Python to get rather complex tasks done with a minimal amount of code.\n",
    "\n",
    "Any temperature data set could be used for this purpose, so you may want to try your own after completing this exercise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's use the global temperature data set to create the warming stripes. We'll need the first filename from the list and we're only interested in the anomalies, so we'll just store this column in a DataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv(fnames[0], index_col=0)['Anomaly (deg C)']\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note once more that `df` is not a DataFrame but a Series. For the purposes of this exercise this does not matter though.\n",
    "\n",
    "To control the colour of the warming stripes, we will use one of `Matplotlib`'s colormaps (see <A href=\"https://matplotlib.org/stable/tutorials/colors/colormaps.html\">https://matplotlib.org/stable/tutorials/colors/colormaps.html</A> for more information). The one that gives the closest resemblance to the original warming stripes is the colormap called bwr (short for blue, white, red). In order to be able to use it, we must first import the `cm` package from `Matplotlib`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from matplotlib import cm\n",
    "\n",
    "cm.bwr(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first line does the importing, the second line calls the `bwr` colormap with 0 as a function argument. The return value is a tuple with four values. The first three are the RGB (red, green and blue) values, the fourth is the alpha value, which we have seen in the example above. Since the color blue is on left side of the colormap, it makes sense that the RGB value is 1.0 (i.e., 100%) for blue. The representation of colours in this way is called the RGBA model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise 3.6***: Execute the code cell below. Does the result make sense?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(cm.bwr(0.5))\n",
    "print(cm.bwr(1.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So how can we use this to create the warming stripes? The first step is to scale all temperature anomalies between 0 and 1. Then we can use `cm.bwr` to get the corresponding colours from the colormap. This is where list comprehension comes in handy again: For each year we can store the colour in a list by looping over all values, scale it between 0 and 1 and convert it to a RGBA colour representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ymin = df.min()\n",
    "ymax = df.max()\n",
    "colors =  [cm.bwr((y - ymin) / (ymax - ymin)) for y in df]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final step is to create a bar graph of the data using `Matplotib`'s `bar` function. Note that all bars have a height of 1 and a widht of 1 (year). The color of each bar is passed using the `color` argument of the `bar` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 1))\n",
    "ax.bar(df.index, np.ones(len(df)), width=1, color=colors)\n",
    "ax.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise 3.7***: Create the same warming stripes but use another colormap. You can find their names on <A href=\"https://matplotlib.org/stable/tutorials/colors/colormaps.html\">https://matplotlib.org/stable/tutorials/colors/colormaps.html</A>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Type your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if you wanted to also plot the actual data in the same graph as the warming stripes? Not a problem: You can use `twinx` to define a secondary y-axis. The code cell below shows how it is done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig, ax0 = plt.subplots(figsize=(8, 1))\n",
    "ax0.bar(df.index, np.ones(len(df)), width=1, color=colors)\n",
    "ax0.axis('off')\n",
    "ax1 = ax0.twinx()\n",
    "ax1.plot(df.index, df.values, 'k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise 3.7***: Add lines of code to the code cell above to (i) set the x axis range from 1850 to 2021 and (ii) add a label to the secondary y axis. *Hint: You must set the x axis range of both ax0 and ax1. Use a `for` loop to do this.*"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "hide_code_all_hidden": false,
  "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.8.8"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "name": "_merged_merged_merged"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}