session_07_plotting.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## From the previous session\n",
    "\n",
    "***Exercise (homework)***: Use the code examples above to create a series of 100 model runs in which the pan factor, rainfall chloride concentration as well as the measured rainfall, evaporation and the dam water volume are all stochastic parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the required packages\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like last time, before we can do any calculations we need to import the file with the measured data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_excel('data/water_balance_data.xlsx',\n",
    "                   index_col=0,\n",
    "                   parse_dates=True)\n",
    "df.head()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the assignment, all water balance components must have a random error. This means that the function `get_conc_cl` must be updated so tha a random error is calculated, and added to the original data. In the code cell below, the rainfall, evaporation and volume measurements are all modulated by adding a (circa) 10% measurement error, which is assumed to be normally distributed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Cl_0 = 20 # g/m^3 = mg/l\n",
    "\n",
    "def get_conc_cl(df, pan_factor=1.3, Cl_rain=2, sd1=0.033):\n",
    "    # Calculate a normally distributed error that is roughly 10% (3 x a standard deviation \n",
    "    # of 3.3%) of the measured value\n",
    "    e_rain = np.random.normal(0, sd1, len(df['rain'])) * df['rain']\n",
    "    e_evap = np.random.normal(0, sd1, len(df['evaporation'])) * df['evaporation']\n",
    "    e_volume = np.random.normal(0, sd1, len(df['volume'])) * df['volume']\n",
    "\n",
    "    # Add the error to the orginal data\n",
    "    df['rain'] = df['rain'] + e_rain\n",
    "    df['evaporation'] = df['evaporation'] + e_evap\n",
    "    df['volume'] = df['volume'] + e_volume\n",
    "\n",
    "    # Proceed with the calculations as before\n",
    "    df['P'] = df['area'] * df['rain'] / 1000.\n",
    "    df['E'] = df['area'] * df['evaporation'] / (1000. * pan_factor)\n",
    "    df['dV'] = -df['volume'].diff(periods=-1)\n",
    "    df['I'] = df['P'] - df['E'] - df['dV']\n",
    "\n",
    "    M_Cl_g = np.empty(len(df))\n",
    "    conc_Cl = np.empty(len(df))\n",
    "\n",
    "    P = df[\"P\"].to_numpy()\n",
    "    I = df[\"I\"].to_numpy()\n",
    "    V = df[\"volume\"].to_numpy()\n",
    "    for i, (Vi, Pi, Ii) in enumerate(zip(V, P, I)):\n",
    "        if i == 0: # First day\n",
    "            M_Cl_g[0] = Vi * Cl_0\n",
    "            conc_Cl[0] = M_Cl_g[0] / Vi # Gives Cl_0 of course!\n",
    "        else:\n",
    "            M_Cl_g[i] = M_Cl_g[i - 1] + dM_P - dM_I\n",
    "            conc_Cl[i] = M_Cl_g[i] / Vi\n",
    "\n",
    "        dM_P = Cl_rain * Pi\n",
    "        dM_I = conc_Cl[i] * Ii\n",
    "\n",
    "    return conc_Cl"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even though the homework exercise asked to run the model 100 times, a number of 99 runs is done. The reason is that this allows drawing realizations for the pan factor and rainfall chloride concentration at exactly 1% intervals between the 1 and 99% percentiles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 99\n",
    "sample_percentiles = np.linspace(0.01, 0.99, N)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that we sliced the array with randomly ordered indices to shuffle the selected values for the pan factors and rainfall chloride concentrations. Otherwise they'd both be increasing in tandem, which would bias the model outcomes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "indices = np.arange(N, dtype=int)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw values from the normally distributed pan factors with a mean of 1.3 and a standard deviation of 0.1 (using a relatively large value to emphasize the differences in model outcomes, the real range is likely to be smaller)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_p = 1.3\n",
    "sigma_p = 0.1\n",
    "\n",
    "dist_p = norm(loc=mu_p, scale=sigma_p)\n",
    "pan_factors = dist_p.ppf(sample_percentiles)\n",
    "\n",
    "np.random.shuffle(indices)\n",
    "pan_factors = pan_factors[indices]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also draw the realizations for the rainfall chloride concentration. Note that the mean and standard deviation were chosen such that the concentrations do not become negative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_c = 4\n",
    "sigma_c = 1\n",
    "\n",
    "dist_c = norm(loc=mu_c, scale=sigma_c)\n",
    "Cl_rains = dist_c.ppf(sample_percentiles)\n",
    "\n",
    "np.random.shuffle(indices)\n",
    "Cl_rains = Cl_rains[indices]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running the model N times is extremely easy if we loop over the elements of the `pan_factors` and `Cl_rains` arrays simultaneously (using `zip`). Note that there is no need to add a counter to keep track of the number of model runs, it simply follows from the number of elements in the arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,4))\n",
    "for pan_factor, Cl_rain in zip(pan_factors, Cl_rains):\n",
    "    Cl = get_conc_cl(df.copy(), pan_factor=pan_factor, Cl_rain=Cl_rain, sd1=0.033)\n",
    "    ax.plot(df.index, Cl)\n",
    "ax.plot(df.index, df[\"Cl_sample\"], 'o')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes though, it is convenient to add a counter. For example, if we wish to store the outcome of each model run in a DataFrame. The script below uses a counter to add a unique column name to each column of `dfnew`, the DataFrame with the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfnew = pd.DataFrame(index=df.index)\n",
    "for n, (pan_factor, Cl_rain) in enumerate(zip(pan_factors, Cl_rains)):\n",
    "    Cl = get_conc_cl(df.copy(), pan_factor=pan_factor, Cl_rain=Cl_rain, sd1=0.033)\n",
    "    col_name = f\"simulation_{n}\"\n",
    "    dfnew[col_name] = Cl"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's try to draw the mean value for each time step, and add some uncertainty bands (rather than plotting 99 colourful lines). One way to do this is to calculate the descriptive statistics using the `describe` function of Pandas. This function works on columns rather than rows so we have to transpose (that is, swap the rows for the columns) the DataFrame first using `dfnew.T`. The result of calling the `describe` method is stored in the DataFrame `dfstats`, which we also transpose so that the statistics for each time step appear per row for each time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfstats = dfnew.T.describe()\n",
    "dfstats = dfstats.T\n",
    "\n",
    "print(dfstats.head)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these numbers we can now plot the median (50% percentile) model outcome as a blue line, and include shaded areas that represent the min-max and 25%-75% percentile ranges. Note the use of the `alpha` value, which controls the transparency of an object in the figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,4))\n",
    "ax.fill_between(dfstats.index, dfstats[\"min\"], dfstats[\"max\"], color='k', alpha=0.1)\n",
    "ax.fill_between(dfstats.index, dfstats[\"25%\"], dfstats[\"75%\"], color='k', alpha=0.1)\n",
    "ax.plot(dfstats.index, dfstats[\"50%\"], color='C0');"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise***: Use the `describe` method to calculate the 16% and 84% percentiles instead of the default values."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise***: By default, Matplotlib draws a line around the filled areas. Look up the documentation and find a way to plot the filled area without a visible line around its borders. Note that there are multiple ways of achieving this."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing data\n",
    "\n",
    "## Matplotlib\n",
    "\n",
    "We have already seen that Matplotlib is heavily object-oriented. It offers many options for drawing, but its vocabulary can be difficult to master. The links below provide some insighs.\n",
    "\n",
    "<A href=\"https://matplotlib.org/stable/tutorials/introductory/quick_start.html\">https://matplotlib.org/stable/tutorials/introductory/quick_start.html</A>\n",
    "\n",
    "<A href=\"https://stackoverflow.com/questions/37970424/what-is-the-difference-between-drawing-plots-using-plot-axes-or-figure-in-matpl\">https://stackoverflow.com/questions/37970424/what-is-the-difference-between-drawing-plots-using-plot-axes-or-figure-in-matpl</A>\n",
    "\n",
    "It is useful to know that an `Artist` is basically everything visible on a `Figure` (the area being drawn on). An `Axes` represents one (sub-)plot in a `Figure`. It contains the plotted data, axis ticks, labels, title, legend, etc. Its methods are the main interface for manipulating the plot. An `Axis` is one of the axes in a graph (for example the x or y axis).\n",
    "\n",
    "This figure from the matplotlib documentation pages is a great help in understanding the matplotlib terminology: <A href=\"https://matplotlib.org/stable/tutorials/introductory/quick_start.html\">https://matplotlib.org/stable/tutorials/introductory/quick_start.html</A>\n",
    "\n",
    "<IMG src=\"https://matplotlib.org/stable/_images/anatomy.png\" width=800>\n",
    "\n",
    "The difference between a `Figure` and an `Axes` may be illustrated by creating a handle to both types of objects using `subplots` (confusingly there is also a function called `subplot`). We color the `Figure` blue, and the `Axes` red."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "fig.set_facecolor('b')\n",
    "ax.set_facecolor('r')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the graph with the model outcomes to explore some ways of prettifying the graph. The code cell below adds a y-axis label, grid lines and formats the dates on the x-axis. For conveniences, here is a summary list of formatting codes for common date and time elements.\n",
    "\n",
    "  * %Y - 4 digit year with upper case Y\n",
    "  * %y - 2 digit year with lower case y\n",
    "  * %m - month as a number with lower case m\n",
    "  * %b - month as abbreviated name with lower case b\n",
    "  * %d - day with lower case d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fig, ax = plt.subplots(figsize=(8,4))\n",
    "ax.fill_between(dfstats.index, dfstats[\"min\"], dfstats[\"max\"], color='k', alpha=0.1)\n",
    "ax.fill_between(dfstats.index, dfstats[\"25%\"], dfstats[\"75%\"], color='k', alpha=0.1)\n",
    "ax.plot(dfstats.index, dfstats[\"mean\"], color='C0')\n",
    "\n",
    "# Add a label on the y-axis\n",
    "ax.set_ylabel(\"Chloride concentration (mg/l)\")\n",
    "\n",
    "# Add grid lines\n",
    "ax.grid(linestyle=':')\n",
    "ax.set_axisbelow(True) # Makes sure that grid lines do not cover graph elements\n",
    "\n",
    "# Format the dates on the x-axis\n",
    "from matplotlib.dates import DateFormatter\n",
    "date_format = DateFormatter(\"%d-%b\")\n",
    "ax.xaxis.set_major_formatter(date_format)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's explore a bit more by plotting the data from the `df` DataFrame. Creating a bar graph is done using `bar`. Rainfall was exceptionally high on one day, so the vertical axis is truncated at 8 mm, and an arrow is added to indicate the value for the bar that extends above this maximum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin = pd.to_datetime(\"2013/08/24\")\n",
    "xmax = pd.to_datetime(\"2013/12/16\")\n",
    "\n",
    "idx = df['rain'] > 8\n",
    "print(df.loc[idx, 'rain'])\n",
    "fig, ax = plt.subplots(figsize=(16, 9))\n",
    "\n",
    "ax.bar(df.index, df[\"rain\"], ec='w', width=1, label=\"Rainfall\")\n",
    "ax.set_ylim(0, 8)\n",
    "ax.set_ylabel(\"Rainfall (mm)\")\n",
    "\n",
    "x0 = pd.to_datetime(\"2013/09/12 15:00\")\n",
    "x1 = pd.to_datetime(\"2013/09/16\")\n",
    "ax.annotate(\n",
    "    text='15 mm', \n",
    "    xy=(x0, 8), \n",
    "    xytext=(x1, 7),\n",
    "    arrowprops=dict(\n",
    "        arrowstyle=\"->\",\n",
    "        connectionstyle=\"angle3,angleA=0,angleB=-90\",\n",
    "    )\n",
    ");"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next code cell plots the measured evaporation. It demonstrates how each bar in the chart can be given a color depending on its value. It is not very useful in this case, but it is a demonstration of the use of a color map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.cm as cm\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(16, 9))\n",
    "\n",
    "scaled_values = (df[\"evaporation\"] - df[\"evaporation\"].min()) / (df[\"evaporation\"].max() - df[\"evaporation\"].min())\n",
    "colors = cm.viridis(scaled_values)\n",
    "ax.bar(df.index, df[\"evaporation\"], color=colors, ec='w', width=1)\n",
    "ax.set_ylabel(\"Pan evaporation (mm)\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise***: Have a look at the available color maps on <A href=\"https://matplotlib.org/stable/tutorials/colors/colormaps.html\">https://matplotlib.org/stable/tutorials/colors/colormaps.html</A>. Modify the previous code example to have the bars plotted using the color scheme of your choice."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The (lengthy) code below draws a three-axes figure. It shows, amongst others, how a horizontal line can be added using `axhline` and how the colors of the vertical axis can be changed. Note how `for` loop is used to set the properties of the `Axes` that should be the same or different for each of the three graphs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (axr, axe, axwl) = plt.subplots(nrows=3, figsize=(16, 9))\n",
    "axt = axr.twinx()\n",
    "\n",
    "axr.bar(df.index, df[\"rain\"], ec='w', width=1, label=\"Rainfall\")\n",
    "axr.set_yticks([0, 4, 8, 12, 16, 20])\n",
    "axr.set_ylim(0, 27 / 25 * 20)\n",
    "axr.set_ylabel(\"Rainfall (mm)\")\n",
    "axt.plot(df[\"temperature\"], 'C3-', label=\"Temperature\")\n",
    "axt.set_yticks([5, 10, 15, 20, 25, 30])\n",
    "axt.set_ylim(5, 32)\n",
    "axt.set_ylabel(\"Temperature (\\u00B0C)\")\n",
    "\n",
    "axe.bar(df.index, df[\"evaporation\"], color=\"C1\", ec='w', width=1)\n",
    "axe.set_ylabel(\"Pan evaporation (mm)\")\n",
    "\n",
    "axwl.plot(df.index, df[\"wl\"], label=\"Water level\")\n",
    "axwl.axhline(-0.5, color='gray', ls='--')\n",
    "axwl.set_ylabel(\"Water level (m)\")\n",
    "axwl.legend(loc='lower center', ncol=2)\n",
    "axwl.text(xmax, -0.5, \"Drainage level \", ha='right', va='bottom')\n",
    "\n",
    "for ax in (axr, axe, axt, axwl):\n",
    "    if ax != axt:\n",
    "        ax.grid(ls=\":\")\n",
    "    ax.set_xlim(xmin, xmax)\n",
    "    if ax != axwl:\n",
    "        ax.set_xticklabels([])\n",
    "\n",
    "axwl.xaxis.set_major_formatter(DateFormatter('%d-%b'))\n",
    "\n",
    "axt.spines['right'].set_color('C3')\n",
    "axt.tick_params(axis='y', colors='C3')\n",
    "axt.yaxis.label.set_color('C3')\n",
    "axt.spines['left'].set_color('C0')\n",
    "axr.spines['left'].set_color('C0')\n",
    "axr.tick_params(axis='y', colors='C0')\n",
    "axr.yaxis.label.set_color('C0')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise***: Create a shared x-axis just like in the upper graph of the previous figure to also include a line that shows the water volume in the farm dam (df[\"volume\"]). Modify the code cell below to create this graph."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now see two lines but without a legend, it is impossible to tell which line is which. Adding a legend to a graph with dual axes is a bit involved, but it can be done."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(16, 3))\n",
    "\n",
    "ax.plot(df.index, df[\"wl\"], label=\"Water level\")\n",
    "ax.axhline(-0.5, color='gray', ls='--')\n",
    "ax.set_ylabel(\"Water level (m)\")\n",
    "ax.legend(loc='lower center', ncol=2)\n",
    "ax.text(xmax, -0.5, \"Drainage level \", ha='right', va='bottom')\n",
    "\n",
    "ax.grid(ls=\":\")\n",
    "ax.set_xlim(xmin, xmax)\n",
    "\n",
    "ax.xaxis.set_major_formatter(DateFormatter('%d-%b'))\n",
    "\n",
    "axv = ax.twinx()\n",
    "l1 = axv.plot(df.index, df[\"volume\"], 'C3', label=\"Volume\")\n",
    "\n",
    "# See: https://matplotlib.org/stable/tutorials/intermediate/legend_guide.html\n",
    "handles0, labels0 = ax.get_legend_handles_labels()\n",
    "handles1, labels1 = axv.get_legend_handles_labels()\n",
    "ax.legend(handles0+handles1, labels0+labels1)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want a consistent font size for all text in your figure, you can set the default using the line of code below. Note that you must include this a the beginning of a script in order to take effect. For help on setting other defaults, see <A href=\"https://matplotlib.org/stable/tutorials/introductory/customizing.html\">https://matplotlib.org/stable/tutorials/introductory/customizing.html</A>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'font.size': 14})"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geopandas_env",
   "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.16"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}