Global Reporting on Groundwater Levels 2022 - Algorithm.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# INTRODUCTION\n",
        "# This program calculates 2 indicators of the quantitative satte of aquifers, based on monitoring data:\n",
        "#  - the rank of the mean water level in the latest year of record, as compared with previous years\n",
        "#  - the pluri-annual trend of mean water levels \n",
        "\n",
        "# It takes in input a set of CSV files containing the water level monitoring data. Each CSV file corresponds to an individual \n",
        "# aquifer or reporting unit. CSV files are formatted as follow: site / date / level\n",
        "# The title of CSV files is the name of the aquifer. For example \"Aquifer_X.csv\"\n",
        "# NB: Site names must start with # if in numerical format\n",
        "\n",
        "# It creates in output:\n",
        "#  - A set of PDF files containing the charts of the water level trends (one PDF file per aquifer). \n",
        "#  - A CSV summary file containing the list of aquifers and for each aquifer: the rank of the qter level in the latest year of \n",
        "#    record, the slope of the trend, as well as the number of observation wells used to calculate the trend.\n",
        "\n",
        "# Acknowledgement: This program was built based on the following code, developed in an earlier IGRAC project:\n",
        "# https://github.com/alexaurgilez/G3P_IGRAC/blob/main/insitu_GWLA_temporal_demo.ipynb\n",
        "\n",
        "# IMPORTING LIBRARIES\n",
        "import numpy as np\n",
        "import matplotlib\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "import datetime\n",
        "import os, glob\n",
        "import warnings; warnings.filterwarnings(action='ignore')\n",
        "\n",
        "# INITIALIZE THE OUTPUT FILE\n",
        "list_aquifers = []\n",
        "\n",
        "# OPENING THE INPUT FILES ONE BY ONE\n",
        "path = \"\"\n",
        "for filename in glob.glob(os.path.join(path, '*.csv')):   # NB: only input files must be in CSV format\n",
        "    \n",
        "    aquifer_name_ext = os.path.basename(filename)\n",
        "    aquifer_name, ext = os.path.splitext(aquifer_name_ext)\n",
        "    print(aquifer_name)\n",
        "    \n",
        "    # READ INPUT FILE AS A DATAFRAME\n",
        "    #gdf = pd.DataFrame(pd.read_csv(filename, index_col='date', parse_dates=True, infer_datetime_format=True, encoding='utf8')) # automatic date parsing created issues\n",
        "    gdf = pd.DataFrame(pd.read_csv(filename)) #without automatic parsing of dates\n",
        "    gdf['date'] = pd.to_datetime(gdf['date'], format = '%d/%m/%Y') # update the date format accordingly\n",
        "    gdf.set_index('date', inplace=True)\n",
        "    \n",
        "    #gdf['level'] = gdf['level'].apply(lambda x: x * 0.3048) # use this line if water levels are expressed in feet\n",
        "    gdf['level'] = gdf['level'].apply(lambda x: -x) # use this line if the water levels are expressed as depth\n",
        "\n",
        "    # REORGANIZE THE DATA WITH ONE COLUMN PER SITE, AND RESAMPLE DATA AS MEAN MONTHLY VALUES \n",
        "    # OVER THE PERIOD 01-01-2013 --> 31/12/2022\n",
        "    gdf_pivot = gdf.pivot_table(index='date',columns='site',values='level')\n",
        "    rng = pd.period_range('2013-01-01',periods=120, freq='M').to_timestamp() - pd.DateOffset(day=31)\n",
        "    gdf_pivot = gdf_pivot.resample('1m').mean().reindex(index=rng)\n",
        "    gdf_pivot.index.name = \"date\"\n",
        "\n",
        "    # RESAMPLE THE DATA OVER THE PERIOD 01-01-2013 --> 31/12/2022. \n",
        "    # In principle the data in input should be within this period, but we resample the data to be sure they are.\n",
        "    gdf_pivot_1 = gdf_pivot.copy()['01-01-2013':'31-12-2022']\n",
        "    gdf_pivot_1.head()\n",
        "\n",
        "    # First glance of the data\n",
        "    # gdf_pivot_1.plot(figsize=(15,5), legend=False, xlim=[datetime.date(2013, 1, 1), datetime.date(2022, 12, 31)])\n",
        "    \n",
        "    # SAVING THE NUMBER OF TIME-SERIES CONTAINED IN THE INPUT FILE\n",
        "    nb_sites_input = gdf_pivot_1.shape[1]\n",
        "\n",
        "    # DATA SELECTION\n",
        "    # This step consists in removing observation wells from the analysis, for instance if they are incomplete. \n",
        "\n",
        "    # Removing time series if the percentage of missing values is larger than some percentage\n",
        "    perc_missing = 8/12 # this is 1 data every 3 months (quarterly data), on average\n",
        "\n",
        "    gdf_1 = gdf_pivot_1.copy()\n",
        "    \n",
        "    max_number_of_nas = perc_missing*gdf_1.shape[0]\n",
        "    gdf_1 = gdf_1.loc[:, (gdf_1.isnull().sum(axis=0) <= max_number_of_nas)]\n",
        "    \n",
        "    # Removing time series if data are missing for more than 1 year\n",
        "    gdf_1_year = gdf_1.resample('1y').mean()\n",
        "    gdf_1_year.index.name = \"date\"\n",
        "    \n",
        "    nb_years_missing = 1\n",
        "\n",
        "    col_to_drop = [i for i in gdf_1_year.columns if gdf_1_year[i].isnull().sum() > nb_years_missing]\n",
        "    gdf_1.drop(col_to_drop, inplace = True, axis=1)\n",
        "\n",
        "    # SAVING THE NUMBER OF TIME-SERIES THAT HAVE BEEN SELECTED\n",
        "    nb_sites = gdf_1.shape[1]\n",
        "\n",
        "    # WE CONTINUE IF AT LEAST ONE SITE HAS BEEN SELECTED\n",
        "    if nb_sites:\n",
        "\n",
        "        # REORGANIZE THE DATAFRAME AS SITE/ DATE / LEVEL\n",
        "        gdf_melt = gdf_1.copy() #dataframe from step 3\n",
        "        gdf_melt = gdf_melt.unstack().reset_index(name='level')\n",
        "        gdf_melt.set_index('site', inplace=True)\n",
        "\n",
        "        # CALCULATE MEAN WATER LEVEL PER SITE\n",
        "        waterlevel_means = gdf_melt.reset_index().groupby('site', as_index=False)['level'].mean()\n",
        "        waterlevel_means.set_index('site', inplace=True)\n",
        "\n",
        "        # ASSIGN MEAN WATER LEVEL TO EACH DATA ROW\n",
        "        gdf_melt['mean_level'] = gdf_melt.index.map(waterlevel_means['level'])\n",
        "\n",
        "        # Calculate normalized water level \"norm_level\"\n",
        "        gdf_melt['norm_level'] = gdf_melt['level'] - gdf_melt['mean_level'] \n",
        "    \n",
        "        # aquifer-averaged dataframe\n",
        "        final_df1 = gdf_melt[['date', 'norm_level']].groupby('date', as_index=False).mean()\n",
        "\n",
        "        # CALCULATE LINEAR TREND\n",
        "        final_df1.loc[:, \"date\"] = pd.to_datetime(final_df1.loc[:, \"date\"], format=\"%d-%b-%y\")\n",
        "        y_values = final_df1.loc[:, \"norm_level\"]\n",
        "        x_values = np.linspace(0,1,len(final_df1.loc[:, \"norm_level\"]))\n",
        "    \n",
        "        #cleaning NaN\n",
        "        idx = np.isfinite(x_values) & np.isfinite(y_values)\n",
        "        slope,b = np.polyfit(x_values[idx], y_values[idx], 1)\n",
        "        \n",
        "        linear_fit = x_values*slope + b\n",
        "        \n",
        "        # PLOTTING\n",
        "        fig, ax = plt.subplots(figsize=(15,5))\n",
        "\n",
        "        gdf_pivot = gdf_melt.pivot_table(index='date',columns='site',values='norm_level')    \n",
        "        plt.plot(gdf_pivot, color='#C0C0C0')     # all the monitoring curves\n",
        "\n",
        "        plt.plot(final_df1.loc[:, \"date\"], final_df1.loc[:, \"norm_level\"], color='#000000') # the composite hydrograph with gaps\n",
        "        \n",
        "        final_df2=final_df1.dropna() # the composite hydrograph without gaps\n",
        "        plt.plot(final_df2.loc[:, \"date\"], final_df2.loc[:, \"norm_level\"],color='black', linestyle=\"dotted\") \n",
        "\n",
        "        plt.plot(final_df1.loc[:, \"date\"], linear_fit, color='#000000', linestyle='dashed') # the linear regression\n",
        "\n",
        "        # axes\n",
        "        ax.set_xlabel('Time', fontsize = 12)\n",
        "        ax.set_ylabel('Normalized groundwater level [masl]', fontsize = 12)\n",
        "        ax.grid()\n",
        "        plt.xlim([datetime.date(2013, 1, 1), datetime.date(2022, 12, 31)])\n",
        "\n",
        "        # title\n",
        "        plt.title(\"Slope: {}. {} observation points available, {} selected\".format(round(slope/10, 3), nb_sites_input, nb_sites))\n",
        "           # the slope is calculated over 0 --> 1, it has to be 2013.01.01 --> 2022.12.31 (10 years)\n",
        "        plt.suptitle(aquifer_name)\n",
        "\n",
        "        # save chart as an image\n",
        "        fig.savefig(path + aquifer_name + \".pdf\", bbox_inches='tight')\n",
        "\n",
        "        # plt.show()\n",
        "        \n",
        "        # CALCULATE PERCENTILE\n",
        "        year_df = final_df1\n",
        "        year_df.set_index('date', inplace=True)\n",
        "        year_df = year_df.resample('Y').mean()\n",
        "        year_df.dropna(inplace = True)\n",
        "        year_df['percent_rank'] = year_df.loc[:,'norm_level'].rank(pct=True)\n",
        "        percentile = year_df.iloc[-1]['percent_rank']\n",
        "        percentile_year = year_df.index[-1]\n",
        "        print(\"percentile: {}\\nyear: {}\".format(percentile, percentile_year))\n",
        "    \n",
        "        # UPDATE THE LIST OF OUTPUTS\n",
        "        list_aquifers.append((aquifer_name, perc_missing, nb_sites_input, nb_sites, round(slope, 3), percentile, percentile_year))\n",
        "    \n",
        "    # IF NO SITE HAS BEEN SELECTED\n",
        "    else:\n",
        "        list_aquifers.append((aquifer_name, perc_missing, nb_sites_input, nb_sites, \"not enough data\")) \n",
        "    \n",
        "# SAVING THE OUTPUT FILE\n",
        "output_df = pd.DataFrame(list_aquifers, columns=['aquifer', 'threshold', 'nb_sites_input', 'nb_sites_selected', 'slope', 'percentile', 'percentile_year'])\n",
        "output_df.to_csv(path + '0_Overview.csv', index=True, encoding='utf-8-sig')\n"
      ],
      "outputs": [],
      "execution_count": null
    }
  ],
  "metadata": {
    "anaconda-cloud": {},
    "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.6.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}