volume_eval_processing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Read NWM output and USGS data, process monthly and daily volume\n",
    "time series, store for subsequent evaluation and plotting\n",
    "\n",
    "Note: This script assumes csv files have been generated that contain\n",
    "NWM AnA hourly channel-rt output times series and USGS (native time step)\n",
    "per specified location. Code would require altering to read time series\n",
    "from alternative source or format.\n",
    "'''\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "### base working directory\n",
    "base_dir = 'C:\\\\repos\\\\git\\\\nwm_eval_projects\\\\annual_runoff'\n",
    "\n",
    "### define list of sites to include in analysis by USGS ID\n",
    "site_list_fname = 'Volumetric_Eval_Sites_3.csv'\n",
    "site_id_header = 'USGS_ID'\n",
    "site_num_header = 'Site_Num' # numbered roughly from East to West\n",
    "\n",
    "### NWM directory (location of hourly AnA time series csv files)\n",
    "nwm_dir = 'ts-data\\\\nwm_ana'\n",
    "nwm_nheadlines = 19\n",
    "\n",
    "### NWM gap summary and aggregated time series output files\n",
    "nwm_gapsum_fname = \"NWM_gaps.csv\"\n",
    "nwm_gaplist_fname = \"NWM_gaps_list.csv\"\n",
    "nwm_accum19_fname = \"NWM_DailyAccum_WY19.csv\"\n",
    "nwm_accum20_fname = \"NWM_DailyAccum_WY20.csv\"\n",
    "nwm_monthly_fname = \"NWM_MonthlyVolume_WYs19-20.csv\"\n",
    "\n",
    "### USGS data directory (raw USGS data csv files)\n",
    "usgs_dir = 'ts-data\\\\usgs'\n",
    "usgs_nheadlines = 21\n",
    "\n",
    "### USGS gap summary and aggregated time series output files\n",
    "usgs_gapsum_fname = \"USGS_gaps.csv\"\n",
    "usgs_gaplist_fname = \"USGS_gaps_list.csv\"\n",
    "usgs_accum19_fname = \"USGS_DailyAccum_WY19.csv\"\n",
    "usgs_accum20_fname = \"USGS_DailyAccum_WY20.csv\"\n",
    "usgs_monthly_fname = \"USGS_MonthlyVolume_WYs19-20.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "### read analysis site information\n",
    "### sites selected by external review, stored in a csv table\n",
    "\n",
    "read_path = os.path.join(base_dir, site_list_fname)\n",
    "sites_df = pd.read_csv(read_path)\n",
    "site_ids =sites_df[site_id_header].astype(\"string\")\n",
    "for i in range(len(site_ids)):\n",
    "    site_ids[i] = site_ids[i].zfill(8) #pad string with zeros to 8 digits\n",
    "sites_df[site_id_header] = site_ids\n",
    "sites_df = sites_df.set_index('Site_Num')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_csv_to_df_dict(read_path, nheadlines, sitelist, q_unit):\n",
    "    '''\n",
    "    read all output/data time series in csv files in read_path \n",
    "    store in dictionary of dataframes if location is in sitelist \n",
    "    return a dictionary containing all data to be processed and plotted\n",
    "    '''\n",
    "    tsdict = {}\n",
    "    q_unit = 'q_' + q_unit\n",
    "    for fname in os.listdir(read_path):\n",
    "        if fname.lower().endswith('.csv'):\n",
    "            # get the location from the filename\n",
    "            extract_start = fname.find(\"_\")+1\n",
    "            location = fname[extract_start:fname.find(\"_\", extract_start+1)]\n",
    "            if not location.isdigit():\n",
    "                continue\n",
    "            if not location in sitelist.values:\n",
    "                continue\n",
    "            fpath = os.path.join(read_path, fname)\n",
    "            df = pd.read_csv(fpath, names=['datetime', q_unit], index_col = 'datetime', skiprows=nheadlines, parse_dates=True)\n",
    "            tsdict[location] = df\n",
    "    return tsdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_data(sites_df, q_unit, tsdict, write_path, filenames):\n",
    "    '''\n",
    "    loop through all sites \n",
    "    check for missing data (>60 min gap), write to file max gap size per location \n",
    "    calculate average daily flow, fill NAN by interpolation, calculate cumulative daily volume by year\n",
    "    '''\n",
    "    # itialize gaplist dataframe\n",
    "    gapstats = sites_df.iloc[:,:1].copy()\n",
    "    gapstats['n_gaps'] = 0\n",
    "    gapstats['tot_days'] = 0\n",
    "    gapstats['max_gap_days'] = 0\n",
    "    gapstats = gapstats.set_index('USGS_ID')\n",
    "    site_ids = sites_df['USGS_ID']\n",
    "\n",
    "    if q_unit == 'cfs':\n",
    "        q_to_AF = 43560.000443512\n",
    "    else:\n",
    "        q_to_AF = 1233.48185532\n",
    "    \n",
    "    q_unit = 'q_' + q_unit\n",
    "    \n",
    "    for location, ts in tsdict.items():\n",
    "\n",
    "        # remove rows with missing data (-999999 or other value < 0)\n",
    "        ts = ts.loc[ts[q_unit] >= 0]\n",
    "        \n",
    "        # calculate deltas between each time step\n",
    "        diff = np.diff(ts.index) \n",
    "        secs = diff / np.timedelta64(1, 's')\n",
    "        minutes = secs / 60\n",
    "        hrs = secs / 3600\n",
    "        days = hrs / 24\n",
    "        ts_with_timedelta = ts[1:].copy() # copy original dataframe beginning at 2nd index as deltas have one less entry\n",
    "        ts_with_timedelta['gap_minutes'] = minutes\n",
    "        ts_with_timedelta['gap_hours'] = hrs\n",
    "        ts_with_timedelta['gap_days'] = days\n",
    "\n",
    "        # find dates with data gaps (delta > 60 min)\n",
    "        dates_w_gaps = ts_with_timedelta[ts_with_timedelta['gap_minutes'].gt(60)].index\n",
    "        gaplist = ts_with_timedelta.loc[dates_w_gaps].copy()\n",
    "        gaplist.insert(0,'prior_q',gaplist[q_unit]) \n",
    "        gaplist.insert(0,'ID', location)\n",
    "\n",
    "        # store max gap size per location in gapstats \n",
    "        gapstats.loc[location, :]  = [len(dates_w_gaps), gaplist['gap_days'].sum(axis = 0), gaplist['gap_days'].max(axis = 0)]\n",
    "        # list each gap per location, with surrounding values and gapsize in ts_w_gaps \n",
    "        for i in range(len(dates_w_gaps)):\n",
    "            ind = np.where(ts_with_timedelta.index==dates_w_gaps[i])\n",
    "            to_col = gaplist.columns.get_loc('prior_q')\n",
    "            from_col = ts_with_timedelta.columns.get_loc(q_unit) \n",
    "            gaplist.iat[i,to_col] = ts_with_timedelta.iat[int(ind[0]) - 1, from_col]\n",
    "\n",
    "        # calculate average daily flow, fill NAN by interpolation\n",
    "        # calculate monthly volume and cumulative daily volume time series\n",
    "        ts_daily_mean = ts_with_timedelta[q_unit].resample('D').mean()  # get mean daily flow\n",
    "        ts_daily_df = ts_daily_mean.to_frame() # write to dataframe\n",
    "        ts_daily_df['interp'] = ts_daily_mean.interpolate() # fill NAN values by interpolation\n",
    "        ts_daily_df['vol_AF'] = ts_daily_df['interp'] * 3600 * 24 / q_to_AF # add column of daily volume in AF  \n",
    "        ts_daily_df['cumvol_AF'] = ts_daily_df['vol_AF'].cumsum() # add column of cumulative sum of daily volume\n",
    "        ts_monthly = ts_daily_df['vol_AF'].resample('M').sum() # get total monthly volumes\n",
    "\n",
    "        wy19 = ts_daily_df.loc[:'2019-09-30'].copy()\n",
    "        wy20 = ts_daily_df.loc['2019-10-01':].copy()\n",
    "        wy20['cumvol_AF'] = wy20['vol_AF'].cumsum()\n",
    "\n",
    "        # if first site, initialize dataframes of processed results\n",
    "        # otherwise append subsequent sites's results as new columns\n",
    "        \n",
    "        if location == site_ids.iloc[0]:\n",
    "            gaplist_allsites = gaplist.copy()\n",
    "            wy19_dailyaccum_allsites = pd.DataFrame(wy19['cumvol_AF'].values, columns = [location], index = wy19.index)\n",
    "            wy20_dailyaccum_allsites = pd.DataFrame(wy20['cumvol_AF'].values, columns = [location], index = wy20.index)\n",
    "            monthly_allsites = pd.DataFrame(ts_monthly.values, columns = [location], index = ts_monthly.index)\n",
    "            monthly_allsites.index = monthly_allsites.index.strftime(\"%b %Y\")\n",
    "        else:\n",
    "            gaplist_allsites = gaplist_allsites.append(gaplist)    \n",
    "            wy19_dailyaccum_allsites[location] = wy19['cumvol_AF'].values\n",
    "            wy20_dailyaccum_allsites[location] = wy20['cumvol_AF'].values\n",
    "            monthly_allsites[location] = ts_monthly.values\n",
    "\n",
    "    gaplist_allsites.to_csv(os.path.join(write_path, filenames[0]))\n",
    "    gapstats.to_csv(os.path.join(write_path, filenames[1]))\n",
    "\n",
    "    wy19_dailyaccum_allsites.to_csv(os.path.join(write_path, filenames[2]))\n",
    "    wy20_dailyaccum_allsites.to_csv(os.path.join(write_path, filenames[3]))\n",
    "    monthly_allsites.to_csv(os.path.join(write_path, filenames[4]))  \n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Read NWM\n",
    "nwm_path = os.path.join(base_dir, nwm_dir)\n",
    "nwm_tsdict = read_csv_to_df_dict(nwm_path, nwm_nheadlines, site_ids, 'cms')\n",
    "\n",
    "## Process NWM\n",
    "nwm_filenames = [nwm_gaplist_fname, nwm_gapsum_fname, nwm_accum19_fname, nwm_accum20_fname, nwm_monthly_fname]\n",
    "process_data(sites_df, 'cms', nwm_tsdict, nwm_path, nwm_filenames)\n",
    "\n",
    "### Read USGS\n",
    "usgs_path = os.path.join(base_dir, usgs_dir)\n",
    "usgs_tsdict = read_csv_to_df_dict(usgs_path, usgs_nheadlines, site_ids, 'cfs')\n",
    "\n",
    "## Process USGS\n",
    "usgs_filenames = [usgs_gaplist_fname, usgs_gapsum_fname, usgs_accum19_fname, usgs_accum20_fname, usgs_monthly_fname]\n",
    "process_data(sites_df, 'cfs', usgs_tsdict, usgs_path, usgs_filenames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}