Rainfall-gauge-analysis.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Untitled0.ipynb",
      "provenance": [],
      "authorship_tag": "ABX9TyMEcqVqAHkhNEo4Sjw5Apks",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/jifarquharson/Farquharson_Amelung_2020_Kilauea-Nature/blob/master/Rainfall-gauge-analysis.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0nLnckFKZI2W"
      },
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import os\n",
        "import glob\n",
        "import datetime as dt\n",
        "import warnings\n",
        "warnings.filterwarnings('ignore')"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6oA2Q6JJZJ_Q"
      },
      "source": [
        "'''\n",
        "Function converts mm input to inches (for plotting figures the correct size).\n",
        "'''\n",
        "\n",
        "def mm2inch(*tupl):\n",
        "    if isinstance(tupl[0], tuple):\n",
        "        return tuple(k*0.0393701 for k in tupl[0])\n",
        "    else:\n",
        "        return tuple(k*0.0393701 for k in tupl)"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 204
        },
        "id": "mcZ72fSvZLjp",
        "outputId": "948bcb04-466f-4e37-81e7-8093184a05c7"
      },
      "source": [
        "'''\n",
        "Read in NCDC gauge data (accesses the GitHub project where we have uploaded the gauge data)\n",
        "'''\n",
        "all_gauges = pd.read_csv(\"https://raw.githubusercontent.com/jifarquharson/Farquharson_Amelung_2020_Kilauea-Nature/master/gauge_data/all_gauges_.csv\")\n",
        "all_gauges.head()"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>STATION</th>\n",
              "      <th>NAME</th>\n",
              "      <th>LATITUDE</th>\n",
              "      <th>LONGITUDE</th>\n",
              "      <th>ELEVATION</th>\n",
              "      <th>DATE</th>\n",
              "      <th>DAPR</th>\n",
              "      <th>MDPR</th>\n",
              "      <th>PRCP</th>\n",
              "      <th>SNOW</th>\n",
              "      <th>SNWD</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>USC00517023</td>\n",
              "      <td>KEAAU 92, HI US</td>\n",
              "      <td>19.63333</td>\n",
              "      <td>-155.03333</td>\n",
              "      <td>85.0</td>\n",
              "      <td>1930-01-01</td>\n",
              "      <td>NaN</td>\n",
              "      <td>NaN</td>\n",
              "      <td>0.0</td>\n",
              "      <td>NaN</td>\n",
              "      <td>NaN</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>USC00517023</td>\n",
              "      <td>KEAAU 92, HI US</td>\n",
              "      <td>19.63333</td>\n",
              "      <td>-155.03333</td>\n",
              "      <td>85.0</td>\n",
              "      <td>1930-01-02</td>\n",
              "      <td>NaN</td>\n",
              "      <td>NaN</td>\n",
              "      <td>3.0</td>\n",
              "      <td>NaN</td>\n",
              "      <td>NaN</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>USC00517023</td>\n",
              "      <td>KEAAU 92, HI US</td>\n",
              "      <td>19.63333</td>\n",
              "      <td>-155.03333</td>\n",
              "      <td>85.0</td>\n",
              "      <td>1930-01-03</td>\n",
              "      <td>NaN</td>\n",
              "      <td>NaN</td>\n",
              "      <td>5.8</td>\n",
              "      <td>NaN</td>\n",
              "      <td>NaN</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>USC00517023</td>\n",
              "      <td>KEAAU 92, HI US</td>\n",
              "      <td>19.63333</td>\n",
              "      <td>-155.03333</td>\n",
              "      <td>85.0</td>\n",
              "      <td>1930-01-04</td>\n",
              "      <td>NaN</td>\n",
              "      <td>NaN</td>\n",
              "      <td>2.3</td>\n",
              "      <td>NaN</td>\n",
              "      <td>NaN</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>USC00517023</td>\n",
              "      <td>KEAAU 92, HI US</td>\n",
              "      <td>19.63333</td>\n",
              "      <td>-155.03333</td>\n",
              "      <td>85.0</td>\n",
              "      <td>1930-01-05</td>\n",
              "      <td>NaN</td>\n",
              "      <td>NaN</td>\n",
              "      <td>0.0</td>\n",
              "      <td>NaN</td>\n",
              "      <td>NaN</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "       STATION             NAME  LATITUDE  LONGITUDE  ...  MDPR PRCP  SNOW  SNWD\n",
              "0  USC00517023  KEAAU 92, HI US  19.63333 -155.03333  ...   NaN  0.0   NaN   NaN\n",
              "1  USC00517023  KEAAU 92, HI US  19.63333 -155.03333  ...   NaN  3.0   NaN   NaN\n",
              "2  USC00517023  KEAAU 92, HI US  19.63333 -155.03333  ...   NaN  5.8   NaN   NaN\n",
              "3  USC00517023  KEAAU 92, HI US  19.63333 -155.03333  ...   NaN  2.3   NaN   NaN\n",
              "4  USC00517023  KEAAU 92, HI US  19.63333 -155.03333  ...   NaN  0.0   NaN   NaN\n",
              "\n",
              "[5 rows x 11 columns]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 3
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0aeg-ZhoZNd3"
      },
      "source": [
        "'''\n",
        "One of the gauges is in a separate file. Read it in and add to the dataframe.\n",
        "'''\n",
        "keaau_2 = pd.read_csv(\"https://raw.githubusercontent.com/jifarquharson/Farquharson_Amelung_2020_Kilauea-Nature/master/gauge_data/keaau_2.csv\")\n",
        "all_gauges = all_gauges.append(keaau_2)"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "A9ztRhz8ZPdT",
        "outputId": "f253eb44-0af1-48e3-c2c3-53350111bff3"
      },
      "source": [
        "datetimes = []\n",
        "for j in all_gauges.DATE:\n",
        "    datetimes.append(dt.datetime.strptime(j, \"%Y-%m-%d\"))\n",
        "all_gauges['dates'] = datetimes\n",
        "\n",
        "gauge_list = [name for name in list(set(all_gauges.NAME))]\n",
        "stn_list = [stn for stn in list(set(all_gauges.STATION))]\n",
        "print(\"List of stations covered: {}\".format(stn_list))"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "List of stations covered: ['USC00516552', 'US1HIHI0060', 'USC00517457', 'USW00021504', 'USC00517023', 'US1HIHI0055', 'USC00511303', 'US1HIHI0008', 'US1HIHI0051', 'USC00513872', 'US1HIHI0003']\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SU4u7FKYZQ61"
      },
      "source": [
        "'''\n",
        "Separate gauges into individual dataframes based on name or number\n",
        "'''\n",
        "\n",
        "kuristown = all_gauges[all_gauges.NAME == 'KURTISTOWN 6.6 SSE, HI US']\n",
        "pahoa_1 = all_gauges[all_gauges.NAME == 'PAHOA 7.5 S, HI US']\n",
        "pahoa_2 = all_gauges[all_gauges.NAME == 'PAHOA 2.7 SSE, HI US']\n",
        "pahoa_3 = all_gauges[all_gauges.NAME == 'PAHOA 2.1 E, HI US']\n",
        "pahoa_4 = all_gauges[all_gauges.NAME == 'PAHOA 65, HI US']\n",
        "hilo = all_gauges[all_gauges.NAME == 'HILO INTERNATIONAL AIRPORT 87, HI US']\n",
        "keaau_1 = all_gauges[all_gauges.STATION == 'USC00517023']\n",
        "keaau_2 = all_gauges[all_gauges.STATION == 'USC00513872']\n",
        "mt_view_1 = all_gauges[all_gauges.NAME == 'MOUNTAIN VIEW 4.5 NNE, HI US']\n",
        "mt_view_2 = all_gauges[all_gauges.NAME == 'MOUNTAIN VIEW 91, HI US']\n",
        "hvnp_1 = all_gauges[all_gauges.NAME == 'HAWAII VOL. NATIONAL PARK HQ. 54, HI US']\n",
        "keaau = keaau_1.append(keaau_2, ignore_index=True) ## append k1 and k2"
      ],
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZAKeRoVbaZsx"
      },
      "source": [
        "'''\n",
        "List of all dataframes for iteration\n",
        "'''\n",
        "df_list = [kuristown, pahoa_1, pahoa_2, pahoa_3, pahoa_4, hilo, keaau_1,keaau_2, mt_view_1, mt_view_2, hvnp_1, keaau]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "60ZvaDBIaZ3n"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RwTX1O8iaaBX"
      },
      "source": [
        "'''\n",
        "Plot daily rainfall. Note that this may take ten minutes or so to execute.\n",
        "'''\n",
        "\n",
        "fig = plt.figure(figsize=mm2inch(180, 300), dpi = 100) \n",
        "for ix, df in enumerate(df_list):\n",
        "    ax = fig.add_subplot(len(df_list), 1, ix+1)\n",
        "    ax.bar(df.dates, df.PRCP, ec = \"k\")\n",
        "    label = [df.STATION.values[0]+\" ({})\".format(df.NAME.values[0]), \"\"]\n",
        "    ax.set_title(label, loc = 'left', fontsize = \"small\")\n",
        "plt.tight_layout()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YOxBNEdBaido"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Uhap_UOnaigc"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "of2xI-2UaijL"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AiZ1MktAail5"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "VQpcezY0ZTJu"
      },
      "source": [
        "def plot_rainfall_timeseries(axis, df, interval = \"30d\", label = False, notes = False, unit_scaler = 1):\n",
        "    ax = axis\n",
        "    index = df.dates ## Use the dates as the index\n",
        "    series = pd.Series(df.PRCP.fillna(0).values, index = index) ## Replace NaN (\"not a number\") values with 0\n",
        "    interval = interval\n",
        "    data = [x/unit_scaler for x in series.rolling(interval).sum()] ## Unit scaler can be used to convert between daily/weekly/monthly average and mm/inches/m etc.\n",
        "    ax.plot(index, data,  alpha= 1, color = \"w\", lw = 1.5, zorder = np.inf) ## Plot timeseries data\n",
        "    ax.plot(index, data,  alpha= 1, color = blue, lw = 0.8, zorder = np.inf)\n",
        "\n",
        "    ## Prettify plot ##\n",
        "    ax.spines[\"top\"].set_visible(False) \n",
        "    ax.spines[\"right\"].set_visible(False)\n",
        "    ax.yaxis.set_tick_params(labelsize=\"x-small\")\n",
        "    ax.tick_params(axis='both', which='major', labelsize=\"x-small\", direction=\"out\")\n",
        "    ax.set_xlim(xmin = \"2010-08-01\")\n",
        "\n",
        "    ## Retrieve mean and standard deviations of data ##\n",
        "    m, std_1, std_1_, std_2, std_2_ = descriptive_stats(df, l = 30)\n",
        "\n",
        "    x_val = \"2018-08-01\"\n",
        "\n",
        "    ax.vlines(x = x_val, ymin = std_2_/unit_scaler, ymax=std_2/unit_scaler,color = blue, lw =.8)\n",
        "\n",
        "    for y in [std_1, std_1_, std_2, std_2_]:\n",
        "        ax.hlines(xmin = \"2018-07-15\", xmax = \"2018-08-15\", y = y/unit_scaler,\n",
        "                 color = blue, lw =.8)\n",
        "    ax.hlines(xmin = \"2018-07-01\", xmax = \"2018-09-01\", y = m/unit_scaler,\n",
        "                 color = blue, lw =.8)\n",
        "   \n",
        "    ax.axhspan(ymin = std_1_/unit_scaler, ymax =std_1/unit_scaler, color= red,alpha = .15, lw =0)\n",
        "    ax.axhspan(ymin = std_2_/unit_scaler, ymax =std_2/unit_scaler, color= red,alpha = .15, lw =0)\n",
        "    ax.axhline(m/unit_scaler,color= red,alpha = .15, lw =.8)\n",
        "    ticks_y = ticker.FuncFormatter(lambda x, pos: '{:.0f}'.format(x))\n",
        "    ax.yaxis.set_major_formatter(ticks_y)\n",
        "    if label == True:\n",
        "        label = df.STATION.values[0]+\" ({})\".format(df.NAME.values[0].split()[0])\n",
        "        ax.annotate(\"{}\".format(label),\n",
        "                    xy = (0,1), xytext = (12,-8),\n",
        "                    xycoords = 'axes fraction', textcoords='offset points',\n",
        "                   fontsize = 8, va = \"center\", ha = \"left\")\n",
        "    if notes == True:\n",
        "#         for i, y in enumerate([m, std_1, std_1_, std_2, std_2_]):\n",
        "#             ax.annotate(r\"${}$\".format(\n",
        "#                 [\"\\overline{x}\",r\"+1\\varsigma\",r\"-1\\varsigma\",r\"+2\\varsigma\",r\"-2\\varsigma\" ][i]\n",
        "#             ),\n",
        "#                        xy = (\"2018-10-01\", y/300),xytext = (10,0),textcoords=\"offset points\",\n",
        "#                         va = \"center\",ha = \"right\", fontsize = \"xx-small\")\n",
        "     for i, y in enumerate([m, std_1, std_2,]):\n",
        "            ax.annotate(r\"${}$\".format(\n",
        "                [\"\\overline{x}\",r\"+1\\varsigma\",r\"+2\\varsigma\",][i]\n",
        "            ),\n",
        "                       xy = (\"2018-10-01\", y/unit_scaler),xytext = (10,0),textcoords=\"offset points\",\n",
        "                        va = \"center\",ha = \"center\", fontsize = \"xx-small\")"
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DdPYgS4VZccH"
      },
      "source": [
        "plot_rainfall_timeseries(axis=axs[0], df=hilo, interval = \"30d\", label = True, notes = True, unit_scaler=300)\n",
        "\n",
        "plot_rainfall_timeseries(axis=axs[1], df=pahoa_3, interval = \"30d\", label = True, notes = True, unit_scaler=300)"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}