IFT_Step02_determine_IFT_mpof.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Original filename: \"02_mpof_diversion_amounts.py\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the following commands in python to install these packages\n",
    "#pip install statsmodels\n",
    "#pip install loess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import datetime as dt\n",
    "import calendar\n",
    "from statsmodels.nonparametric.smoothers_lowess import lowess as lw\n",
    "from loess.loess_1d import loess_1d as ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paths, input, and output file names\n",
    "\n",
    "outer_dir = os.path.abspath(os.path.join(os.getcwd() ,\"../..\"))\n",
    "refpth = os.path.join(outer_dir, 'IFT_files', 'Reference Files')\n",
    "sfe_char_csv = os.path.join(refpth, 'TNC SFE LOI Characteristics with MAF.csv')\n",
    "sfe_char_xls= os.path.join(refpth, 'TNC SFE LOI Characteristics.xlsx')\n",
    "subset_xls = os.path.join(refpth, 'TNC Subset SFE LOIs.xlsx')\n",
    "#sfe_char_csv = os.path.join(refpth, 'All SFE LOI Characteristics with MAF.csv')\n",
    "#sfe_char_xls= os.path.join(refpth, 'All SFE LOI Characteristics.xlsx')\n",
    "#subset_xls = os.path.join(refpth, 'Subset SFE LOIs.xlsx')\n",
    "\n",
    "unimpath = os.path.join(outer_dir, \"IFT_files\", 'Unimpaired Flow')\n",
    "comid_csv = os.path.join(refpth, 'TNC-POI-COMID-20210804.csv')\n",
    "#comid_csv = os.path.join(refpth, 'SFER-POI-COMID-16Jun2020.csv')\n",
    "wytdir = os.path.join(unimpath, \"Water Year Types\")\n",
    "wmtdir = os.path.join(unimpath, \"Water Month Types\")\n",
    "\n",
    "startdir = os.path.join(outer_dir, 'IFT_files', 'IFT Results')\n",
    "wmtfile = os.path.join(wmtdir, \"LOI 9999 WMT.csv\") #file containing WMTs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "avgwin = 30 #window size of moving average\n",
    "switchdays = 30 # number of days to look from beginning and end of list of dates to determine when to switch between LOESS and moving average\n",
    "\n",
    "startdate = dt.datetime(1950,10,1)\n",
    "stopdate = dt.datetime(2021,9,30)\n",
    "\n",
    "#startdate = dt.datetime(1995,10,1)\n",
    "#stopdate = dt.datetime(2017,9,30)\n",
    "alldates = pd.date_range(startdate,stopdate)\n",
    "\n",
    "# The following variables can have a single value or a list of values (i.e., [#1, #2, #3,...]). The script loops through\n",
    "# however many values are put in\n",
    "exd_perc_flows = [0.1] #[0.1, 0.2, 0.3] #the exceedance percentile flow to create the streamflow baseline (default = 10% or 0.1)\n",
    "divert_ratios = [0.1, 0.2] #[0.1, 0.2, 0.3] #proportion of streamflow baseline for setting diversion allocation (default = 10% or 0.1)\n",
    "prd_pct_reqts = [0.1] #percentile of the length of time the requirement applies, as specified in reqt_time. For example,\n",
    "# if this is 0.1, the 10th percentile of all flows\n",
    "reqt_times = [0] #[0,1,2,3] #length of time for requirements to apply. 0 for daily, 1 for weekly, 2 for semi-monthly (1st and 15th), 3 for monthly\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Original filename: \"00_add_ts_col.py\"\n",
    "def add_ts_col(tab):\n",
    "#adds column for TS as first column of DataFrame as formatted date\n",
    "    tab['TS'] = tab.index.strftime('%m/%d/%Y')\n",
    "    cols = tab.columns.tolist()\n",
    "    cols = cols[-1:] + cols[:-1]\n",
    "    tab = tab[cols]\n",
    "    return tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Original filename: \"00_get_all_sfe_lois.py\"\n",
    "# Read in SFE LOI characteristics table and calculate bankfull flow using cont. Area relation\n",
    "\n",
    "def get_all_sfe_lois():\n",
    "    sfelois = pd.read_excel(sfe_char_xls,index_col=0)\n",
    "    loi = [str(i) for i in sfelois['LOI']]\n",
    "    sfelois['Outlet LOI'] = sfelois['LOI']\n",
    "    sfelois['Contributing Area (mi^2)'] = sfelois['Contributing Area']\n",
    "    #sfelois['MAF'] = sfelois['Mean Annual Flow (cfs)']\n",
    "    sfelois['Qbf'] = 71.5 * sfelois['Contributing Area (mi^2)'] #71.5 cfs/mi^2 according to Darren\n",
    "\n",
    "    subset = pd.read_excel(subset_xls)\n",
    "    sublois = sfelois.loc[subset['SWSID'],:]\n",
    "    return sublois, loi, sfelois"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Original filename: \"00_read_loi_paradigm_flow_v3.py\"\n",
    "def read_loi_paradigm_flow(p):\n",
    "    unimpflowfile = os.path.join(unimpath, str(p) + '.csv')\n",
    "    unimp = pd.read_csv(unimpflowfile,index_col=0)\n",
    "    unimp.index = pd.DatetimeIndex(unimp.index)\n",
    "    return unimp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Oct-06-2021 04:37:23 PM\n",
      "MPOF - Default Complete\n",
      "MPOF - 20% of 10th Percentile Hydrograph, Daily Complete\n",
      "Oct-06-2021 04:37:28 PM\n"
     ]
    }
   ],
   "source": [
    "loitab, loi, fultab = get_all_sfe_lois()\n",
    "print(dt.datetime.now().strftime('%b-%d-%Y %I:%M:%S %p'))\n",
    "#loop through all selections\n",
    "for e in exd_perc_flows:\n",
    "    for d in divert_ratios:\n",
    "        for r in reqt_times:\n",
    "            for p in prd_pct_reqts:\n",
    "                mpofout = pd.DataFrame(columns=loi, index=alldates)\n",
    "                if (not r == 0) or ((r == 0) and (p == prd_pct_reqts[len(prd_pct_reqts)-1])): #don't do daily for every different p\n",
    "                    #determine MPOF variation labels\n",
    "                    exdstr = str(int(e * 100)) + 'th Percentile Hydrograph, '\n",
    "                    pctstr = str(int(d * 100)) + '%'\n",
    "                    prdstr = str(int(p * 100)) + 'th Percentile of '\n",
    "                    if r == 0:\n",
    "                        lenstr = 'Daily'\n",
    "                        prdstr = ''\n",
    "                    elif r == 1:\n",
    "                        lenstr = 'Weekly'\n",
    "                    elif r == 2:\n",
    "                        lenstr = 'Biweekly' #Note: not perfectly biweekly, but new requirements set on the 1st and 15th of each month\n",
    "                    elif r == 3:\n",
    "                        lenstr = 'Monthly'\n",
    "                    methstr = pctstr+ ' of ' +exdstr+ prdstr+lenstr\n",
    "                    if (e == 0.1) & (d == 0.1) & (r == 0):\n",
    "                        methstr = 'Default'\n",
    "                    methstr = 'MPOF - ' + methstr\n",
    "                    for l in loi:\n",
    "                        unimp = read_loi_paradigm_flow(l) #read in unimpaired flow\n",
    "\n",
    "                        # determine exceedance hydrograph from unimpaired flow\n",
    "                        flow90exd = unimp['flow'].groupby(by=[unimp.index.month, unimp.index.day]).quantile(e)\n",
    "                        #fix for water year order of months\n",
    "                        prevyr = flow90exd.loc[np.arange(10,13)]\n",
    "                        flow90exd = prevyr.append(flow90exd.loc[np.arange(1, 10)])\n",
    "\n",
    "                        tempind = [dt.datetime(1995,10,1)+dt.timedelta(x) for x in range(366)]\n",
    "                        flow90exd.index = tempind\n",
    "                        loessind = flow90exd.index\n",
    "                        loessfilt = lw(flow90exd.get_values(), loessind, frac=70 / 365, it=0) #70-day loess filter\n",
    "\n",
    "                        # moving average\n",
    "                        flowmovavg= flow90exd.append(flow90exd)\n",
    "                        mastr = str(avgwin) + '-Day MA'\n",
    "                        movavg_dups = flowmovavg.rolling(avgwin,center=True).mean().dropna().sort_index() #calculate rolling/moving avg\n",
    "                        movavg = movavg_dups.loc[~movavg_dups.index.duplicated(keep='first')] #remove duplicates\n",
    "                        #format into table\n",
    "                        flowbase = pd.DataFrame({'Daily 90% Exceedance': flow90exd, mastr: movavg, 'LOESS': loessfilt[:,1]},index=flow90exd.index)\n",
    "                        methdiff = np.abs(flowbase[mastr]-flowbase['LOESS']) #difference between moving avg and loess\n",
    "                        #determine day within beginning and ending number of days to swtich from loess to moving avg (default = 30)\n",
    "                        headswitch = methdiff.index[1:switchdays+1][methdiff.iloc[1:switchdays+1] == min(methdiff.iloc[1:switchdays+1])]\n",
    "                        tailswitch = methdiff.index[-(switchdays+1):-1][methdiff.iloc[-(switchdays+1):-1] == min(methdiff.iloc[-(switchdays+1):-1])]\n",
    "                        #determine streamflow baseline from combination of methods\n",
    "                        flowbase['Daily Streamflow Baseline (cfs)'] = flowbase.loc[min(flowbase.index):pd.DatetimeIndex(headswitch-pd.DateOffset(1)).to_pydatetime()[0],mastr].append(\n",
    "                            flowbase.loc[headswitch.to_pydatetime()[0]:(tailswitch-pd.DateOffset(1)).to_pydatetime()[0],'LOESS']).append(\n",
    "                            flowbase.loc[tailswitch.to_pydatetime()[0]:max(flowbase.index),mastr])\n",
    "\n",
    "\n",
    "                        if r == 0: #daily, no resampling needed\n",
    "                            reindfb = flowbase.copy()\n",
    "                        else:\n",
    "                            if r == 1: # weekly\n",
    "                                freqdates = [min(flowbase.index) + (dt.timedelta(days=7) * i) for i in range(0, 53)]\n",
    "                                gbinds = np.zeros_like(flowbase.index,dtype='int')\n",
    "                                for i in range(0, len(freqdates)):\n",
    "                                    gbinds[(flowbase.index < freqdates[i]) & (gbinds == 0)] = i\n",
    "                                gbinds[gbinds == 0] = 53\n",
    "                                indfq = 'W'\n",
    "                            elif r == 2: #biweekly\n",
    "                                days = np.array([1,15])\n",
    "                                months = np.append(np.arange(10,13),np.arange(1, 10))\n",
    "                                years = np.append(1995 * np.ones(6), 1996 * np.ones(18)).astype(int) ## NOT sure why 1996 and 1995\n",
    "                                dates = np.array(np.meshgrid(months, days)).T.reshape(-1, 2)\n",
    "                                freqdates = [dt.datetime(years[i],dates[i,0], dates[i,1]) for i in range(len(dates))]\n",
    "                                gbinds = np.zeros_like(flowbase.index,dtype='int')\n",
    "                                for i in range(0, len(freqdates)):\n",
    "                                    gbinds[(flowbase.index < freqdates[i]) & (gbinds == 0)] = i\n",
    "                                gbinds[gbinds == 0] = 24\n",
    "                                indfq = 'SMS'\n",
    "                            elif r == 3: #monthly\n",
    "                                gbinds = flowbase.index.month\n",
    "                                indfq = 'MS'\n",
    "                            reindfb = flowbase.groupby(gbinds, sort=False).quantile(p) #group by time period and take quantile\n",
    "                            reindfb.index = pd.date_range(min(tempind), max(tempind), freq=indfq)\n",
    "                            reindfb = reindfb.reindex(flowbase.index, method='ffill')\n",
    "                        #calculate diverion allocation, both daily and resampled (may be the same if both are daily)\n",
    "                        flowbase['Daily Diversion Allocation (cfs)'] = d * flowbase['Daily Streamflow Baseline (cfs)']\n",
    "                        flowbase['Resampled Streamflow Baseline (cfs)'] = reindfb['Daily Streamflow Baseline (cfs)']\n",
    "                        flowbase['Resampled Diversion Allocation (cfs)'] = d * reindfb['Daily Streamflow Baseline (cfs)']\n",
    "\n",
    "                        #save off calculations performed so results can be analyzed manually\n",
    "                        if not(os.path.exists(os.path.join(startdir, methstr))):\n",
    "                            os.mkdir(os.path.join(startdir, methstr))\n",
    "                        flowbase.to_csv(os.path.join(startdir, methstr, 'LOI ' + l + ' MPOF ' + methstr + ' Data.csv'))\n",
    "\n",
    "                        #now need to take calculated diversion allocation and subtract it  from unimpaired flow to get IFT\n",
    "                        fbst = flowbase.copy()\n",
    "                        fbst.index = flowbase.index.strftime('%m-%d')\n",
    "                        years = np.arange(min(alldates.year), max(alldates.year)+1)\n",
    "                        #need to copy diversion allocation to all years\n",
    "                        fballyr = pd.concat([fbst]*(years[-1]-years[0]))\n",
    "                        alyind = (np.concatenate([np.repeat(years[0],92), \\\n",
    "                                                  np.repeat(years[1:-1],366), \\\n",
    "                                                  np.repeat(years[-1],274)]) \\\n",
    "                                  .astype(str)+fballyr.index)\n",
    "                        fballyr.index = alyind\n",
    "                        for y in years[1:]:\n",
    "                            if ~calendar.isleap(y):\n",
    "                                fballyr.drop(index=str(y)+'02-29',inplace=True) #remove feb 29 where it doesn't exist\n",
    "                        fballyr.index = pd.to_datetime(fballyr.index,format='%Y%m-%d')\n",
    "                        mpofout.loc[mpofout.index,l] = unimp.loc[mpofout.index,'flow'] - fballyr.loc[mpofout.index,'Resampled Diversion Allocation (cfs)']\n",
    "\n",
    "                        mpofout.loc[mpofout[l] < 0, l] = 0 #can't have negative IFT\n",
    "                    #format for WEAP and save\n",
    "                    mpofout = add_ts_col(mpofout)\n",
    "                    savdir = startdir\n",
    "                    if not ((e == 0.1) & (d == 0.1) & (r == 0)):\n",
    "                        savdir = os.path.join(startdir, 'MPOF Variants')\n",
    "                        if not(os.path.exists(savdir)):\n",
    "                            os.mkdir(savdir)\n",
    "                    mpofout.to_csv(os.path.join(savdir, 'All LOI ' +methstr+ ' IFTs.csv'))\n",
    "                    print(methstr + ' Complete')\n",
    "\n",
    "print(dt.datetime.now().strftime('%b-%d-%Y %I:%M:%S %p'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (freshwater)",
   "language": "python",
   "name": "freshwater"
  },
  "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}