Effective Discharge App.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the map to select a current or historic gaging station. Hovering over station provides station name and USGS ID. Clicking the station runs the effective discharge analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "%matplotlib widget\n",
    "from matplotlib import pyplot as plt\n",
    "plt.ioff()\n",
    "import matplotlib.dates as mdates\n",
    "from ipyleaflet import Map, Marker, MarkerCluster, WMSLayer, Unicode, basemaps\n",
    "from ipysheet import from_dataframe\n",
    "import numpy as np\n",
    "import ipywidgets as widgets\n",
    "wms = WMSLayer(\n",
    "    url='https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryTopo/MapServer/tile/{z}/{y}/{x}',\n",
    "    attribution='US Geological Survey'\n",
    ")\n",
    "\n",
    "m = Map(basemap=basemaps.CartoDB.Positron, center=(30.661, -88.645), zoom=5)\n",
    "m.add_layer(wms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,(daily_flow_plot, inst_flow_plot, ssc_plot) = plt.subplots(3,1, figsize=(12, 8),constrained_layout=True);\n",
    "fig_hist,(flow_hist, flow_hist_log) = plt.subplots(1,2, figsize=(12, 8),constrained_layout=True);\n",
    "fig_rating_curve, rating_curve = plt.subplots(1,1, figsize=(12, 8),constrained_layout=True);\n",
    "fig_eff_discharge, (eff_discharge, eff_discharge_log) = plt.subplots(2,1,figsize=(12,8),constrained_layout=True);\n",
    "cum_discharge = eff_discharge.twinx()  # instantiate a second axes that shares the same x-axis\n",
    "cum_discharge_log = eff_discharge_log.twinx()  # instantiate a second axes that shares the same x-axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "98fdcccb7b844529acea0c0740c228c0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=1.0, bar_style='info', max=1.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8d7409ef6dc141a0822104bd8555e6b0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Accordion(children=(Map(center=[30.661, -88.645], controls=(ZoomControl(options=['position', 'zoom_in_text', '…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from tqdm.notebook import tqdm\n",
    "pbar = tqdm()\n",
    "children = [m,\n",
    "            fig.canvas,\n",
    "            fig_hist.canvas, \n",
    "            fig_rating_curve.canvas, \n",
    "            fig_eff_discharge.canvas,\n",
    "            widgets.VBox(),\n",
    "            widgets.VBox()]\n",
    "tab = widgets.Accordion(children=children)\n",
    "DISCRETE_TAB=5\n",
    "DAILY_TAB=6\n",
    "tab.set_title(0,'Map')\n",
    "tab.set_title(1,'Data Plots')\n",
    "tab.set_title(2,'Flow Frequency')\n",
    "tab.set_title(3,'Rating Curve')\n",
    "tab.set_title(4,'Effective Discharge')\n",
    "tab.set_title(DISCRETE_TAB,'Discrete Data')\n",
    "tab.set_title(DAILY_TAB,'Daily Data')\n",
    "tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests,pandas,io\n",
    "state_data=[]\n",
    "for state in ['la','ms','ar','al','tn','mo']:\n",
    "    resp = requests.post('https://waterservices.usgs.gov/nwis/site/?format=rdb&stateCd={0}&siteType=ST&siteStatus=active&hasDataTypeCd=dv,iv'.format(state))\n",
    "    #resp = requests.post('https://waterservices.usgs.gov/nwis/site/?format=rdb&siteType=ST&siteStatus=active&hasDataTypeCd=dv,iv')\n",
    "    #resp = requests.post('https://waterservices.usgs.gov/nwis/site/?format=rdb&stateCd={0}&parameterCd=80154&siteStatus=active'.format(state))\n",
    "    #resp = requests.post('https://waterservices.usgs.gov/nwis/site/?format=rdb&stateCd={0}&parameterCd=80154&siteStatus=all&hasDataTypeCd=iv,dv'.format(state))\n",
    "    ncomments=0\n",
    "    for line in resp.text.split('\\n'):\n",
    "        if '#' in line:\n",
    "            ncomments+=1\n",
    "        else:\n",
    "            break\n",
    "    state_data.append(pandas.read_table(io.BytesIO(resp._content),comment='#',skiprows=[ncomments+1]))\n",
    "f = pandas.concat(state_data, axis=0)\n",
    "#f.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_station_daily(site):\n",
    "    iterable = range(13)\n",
    "    pbar.reset(total=len(iterable))\n",
    "    site_no=site[0]\n",
    "    fig.suptitle(site[1])\n",
    "    fig_hist.suptitle(site[1])\n",
    "    fig_rating_curve.suptitle(site[1])\n",
    "    fig_eff_discharge.suptitle(site[1])\n",
    "    daily_flow_plot.clear()\n",
    "    inst_flow_plot.clear()\n",
    "    ssc_plot.clear()\n",
    "    flow_hist.clear()\n",
    "    flow_hist_log.clear()\n",
    "    rating_curve.clear()\n",
    "    eff_discharge.clear()\n",
    "    eff_discharge_log.clear()\n",
    "    cum_discharge_log.clear()\n",
    "    cum_discharge.clear()\n",
    "    import re\n",
    "    import requests \n",
    "    import xml.etree.ElementTree as ET\n",
    "    from zipfile import ZipFile\n",
    "    from io import BytesIO\n",
    "    import numpy as np\n",
    "    import math\n",
    "    global station_zipfile\n",
    "    # defining the api-endpoint\n",
    "    URL = \"https://cida.usgs.gov/\"\n",
    "    API_ENDPOINT = \"https://cida.usgs.gov/sediment/data\"\n",
    "    ogcfilter = ET.Element(\"ogc:Filter\",attrib={\"xmlns:ogc\":\"http://www.opengis.net/ogc\"})\n",
    "    property_is = ET.SubElement(ogcfilter, \"ogc:PropertyIsEqualTo\", attrib={\"matchCase\":\"true\"})\n",
    "    ET.SubElement(property_is, \"ogc:PropertyName\").text = \"SITE_NO\"\n",
    "    ET.SubElement(property_is, \"ogc:Literal\").text = '{0:08d}'.format(site_no)\n",
    "    ogcfilter_string = ET.tostring(ogcfilter).decode()  \n",
    "    data = {'email' : 'cekees@gmail.com',\n",
    "            'format': 'tsv',\n",
    "            'directDownload':'true',\n",
    "            'dataTypes': 'sites_data_daily_discrete_',\n",
    "            'dailyFilter' : ogcfilter_string,\n",
    "            'discreteFilter' : ogcfilter_string\n",
    "           }\n",
    "    # sending post request and saving response as response object\n",
    "    pbar.update(1)\n",
    "    search = requests.post(url = API_ENDPOINT, data = data)  \n",
    "    download = requests.get(url=URL+search.text)\n",
    "    pbar.update(2)\n",
    "    station_zipfile = ZipFile(BytesIO(download._content))\n",
    "    def get_comments(filename):\n",
    "        ncomments=0\n",
    "        comments=''\n",
    "        with station_zipfile.open(filename) as f:\n",
    "            commentline = f.readline().decode()\n",
    "            while commentline:\n",
    "                if commentline[0] == '#' or commentline[:2] == r'\"#':\n",
    "                    ncomments +=1\n",
    "                    comments += commentline.strip(r'\"').strip('#')\n",
    "                    commentline=f.readline().decode()\n",
    "                else:\n",
    "                    commentline=None\n",
    "        return (ncomments,comments)   \n",
    "    pbar.update(3)\n",
    "    if station_zipfile.filelist:\n",
    "        has_daily_file=False\n",
    "        if 'daily_data.tsv' in station_zipfile.namelist():\n",
    "            has_daily_file=True\n",
    "            (ncomments, comments) = get_comments('daily_data.tsv')\n",
    "            f = pandas.read_table(station_zipfile.open('daily_data.tsv'),header=ncomments)\n",
    "            f.datetime = pandas.to_datetime(f.datetime)\n",
    "            drop_columns = []\n",
    "            for col in f.columns:\n",
    "                if col not in ['DAILY_FLOW','datetime']:\n",
    "                    drop_columns.append(col)\n",
    "            tab.children[DAILY_TAB].children = [from_dataframe(f.drop(columns=drop_columns))]\n",
    "            #print(f.head())\n",
    "            if 'DAILY_FLOW' in f.columns:\n",
    "                daily_flow_plot.plot(f.datetime,f.DAILY_FLOW)\n",
    "                daily_flow_plot.set_ylabel('Discharge [cfs]')\n",
    "                daily_flow_plot.set_xlabel('Date')\n",
    "                daily_flow_plot.set_title('Daily')\n",
    "                daily_flow_plot.grid()\n",
    "                daily_flow = []\n",
    "                for flow in f.DAILY_FLOW:\n",
    "                    if math.isfinite(flow):\n",
    "                        daily_flow.append(flow)\n",
    "                daily_flow = np.array(daily_flow)\n",
    "            else:\n",
    "                daily_flow = None\n",
    "        pbar.update(4)\n",
    "        if 'discrete_data.tsv' in station_zipfile.namelist():\n",
    "            (ncomments, comments) = get_comments('discrete_data.tsv')\n",
    "            f = pandas.read_table(station_zipfile.open('discrete_data.tsv'),header=ncomments)\n",
    "            drop_columns = []\n",
    "            for col in f.columns:\n",
    "                if col not in ['DAILYFLOW','DATETIME','INSTFLOW','SSC','P63']:\n",
    "                    drop_columns.append(col)\n",
    "            tab.children[DISCRETE_TAB].children = [from_dataframe(f.drop(columns=drop_columns))]\n",
    "            #for col in f.columns: \n",
    "            #    print(col) \n",
    "            f.DATETIME = pandas.to_datetime(f.DATETIME)\n",
    "            #print(f.head())\n",
    "            if 'DAILYFLOW' in f.columns:\n",
    "                if not has_daily_file:\n",
    "                    daily_flow_plot.plot(f.DATETIME, f.DAILYFLOW)\n",
    "                    daily_flow_plot.set_ylabel('Discharge [cfs]')\n",
    "                    daily_flow_plot.set_xlabel('Date')\n",
    "                    daily_flow_plot.set_title('Daily Flow (from discrete file)')\n",
    "                    daily_flow_plot.grid()\n",
    "                    daily_flow = []\n",
    "                    for flow in f.DAILYFLOW:\n",
    "                        if math.isfinite(flow):\n",
    "                            daily_flow.append(flow)\n",
    "                    daily_flow = np.array(daily_flow)\n",
    "            pbar.update(5)\n",
    "            isfinite = np.isfinite(f.DAILYFLOW)\n",
    "            daily_flow_inst_file = f.DAILYFLOW[isfinite]\n",
    "            if 'DATETIME' in f.columns and 'INSTFLOW' in f.columns:\n",
    "                inst_flow_plot.plot(f.DATETIME,f.INSTFLOW)\n",
    "                inst_flow_plot.set_title(\"Discrete\")\n",
    "                inst_flow_plot.set_xlabel('Date')\n",
    "                inst_flow_plot.set_ylabel('Discharge [cfs]')\n",
    "                inst_flow_plot.grid()\n",
    "            pbar.update(6)\n",
    "            if 'SSC' in f.columns and 'P63' in f.columns and 'INSTFLOW' in f.columns and 'DAILYFLOW' in f.columns:\n",
    "                flow_sand=[]\n",
    "                p63_conc=[]\n",
    "                ssc=[]\n",
    "                sand=[]\n",
    "                for i, (insti, dailyi, p63i, ssci) in enumerate(zip(f.INSTFLOW, f.DAILYFLOW, f.P63, f.SSC)):\n",
    "                    hasData=True\n",
    "                    insti = float(insti)\n",
    "                    dailyi = float(dailyi)\n",
    "                    if type(p63i) is str:\n",
    "                        p63i = float(re.sub('[^\\d*\\.?\\d]+','', p63i))\n",
    "                    else:\n",
    "                        p63i = float(p63i)\n",
    "                    if type(ssci) is str:\n",
    "                        ssci = float(re.sub('[^\\d*\\.?\\d]+','', ssci))\n",
    "                    else:\n",
    "                        ssci = float(ssci)\n",
    "                    if math.isfinite(ssci) and math.isfinite(p63i) and (math.isfinite(insti) or math.isfinite(dailyi)):\n",
    "                        ssc.append(ssci)\n",
    "                        p63_conc.append(ssci*(p63i/100.0))\n",
    "                        sandi = ssci*(1.0 - p63i/100.0)\n",
    "                        sand.append(sandi)\n",
    "                        if math.isfinite(dailyi):\n",
    "                            flow_sand.append(dailyi)\n",
    "                        else:\n",
    "                            flow_sand.append(insti)       \n",
    "                pbar.update(7)\n",
    "                flow_sand = np.array(flow_sand)\n",
    "                sand = np.array(sand)\n",
    "                p63_con = np.array(p63_conc)\n",
    "                ssc = np.array(ssc)\n",
    "                ssc_plot.scatter(flow_sand,ssc,label='SSC')\n",
    "                ssc_plot.scatter(flow_sand,p63_con, label='Silt')\n",
    "                ssc_plot.scatter(flow_sand,sand, label='Sand')\n",
    "                ssc_plot.set_xlabel('Discharge [cfs]')\n",
    "                ssc_plot.set_ylabel('Sediment [mg/l]')\n",
    "                ssc_plot.set_title(\"Discrete Sediment Data\")\n",
    "                ssc_plot.grid()\n",
    "                ssc_plot.legend()\n",
    "                pbar.update(8)\n",
    "                if daily_flow.shape[0] < flow_sand.shape[0]:\n",
    "                    daily_flow = flow_sand\n",
    "                (n,bins,patches) = flow_hist.hist(daily_flow,bins=25)\n",
    "                flow_hist.set_ylabel('Incidence')\n",
    "                flow_hist.set_xlabel('Discharge [cfs]');\n",
    "                flow_hist.set_title(\"Flow Frequency\");\n",
    "                annual_frequency = 365.0*n/n.sum()\n",
    "                rep_discharge = 0.5*(bins[0:-1] + bins[1:])\n",
    "                log_daily_flow = np.log10(daily_flow)\n",
    "                (n_log,bins_log,patches_log) = flow_hist_log.hist(log_daily_flow,bins=25)\n",
    "                flow_hist_log.set_ylabel('Incidence')\n",
    "                flow_hist_log.set_xlabel(r'Log Discharge [log(cfs)]');\n",
    "                flow_hist_log.set_title(\"Flow Frequency\");\n",
    "                pbar.update(9)\n",
    "                annual_frequency_log = 365.0*n_log/n_log.sum()\n",
    "                discharge_bins = 10**(bins_log)\n",
    "                rep_discharge_log = 0.5*(discharge_bins[0:-1] + discharge_bins[1:])\n",
    "                from scipy import stats\n",
    "                K=.0027 #conversion factor to get Q_s from discharge and concentration in mg/l\n",
    "                q_sand = flow_sand*sand*K\n",
    "                sand_indeces = np.where(sand > 0.0)\n",
    "                log_flow_sand = np.log10(flow_sand[sand_indeces])\n",
    "                log_sand = np.log10(sand[sand_indeces])\n",
    "                sand_pow = stats.linregress(x=log_flow_sand,y=log_sand)\n",
    "                pow = pandas.DataFrame([sand_pow], \n",
    "                       columns=[\"slopes\", \n",
    "                          \"intercepts\", \n",
    "                          \"r_value\", \n",
    "                          \"p_value\", \n",
    "                          \"std err\"],\n",
    "                       index=['Sand'])\n",
    "                import math\n",
    "                rating_curve.scatter(flow_sand,sand,c='b')\n",
    "                rating_curve.plot(flow_sand,10**(sand_pow[1]+sand_pow[0]*np.log10(flow_sand)),'b')\n",
    "                rating_curve.set_title('Rating Curve')\n",
    "                rating_curve.set_xlabel(r'Discharge [(cfs)]')\n",
    "                rating_curve.set_ylabel(r'Sand Load [(tons/day)]');\n",
    "                rating_curve.set_yscale('log')\n",
    "                rating_curve.set_xscale('log')\n",
    "                rating_curve.grid()\n",
    "                pbar.update(10)\n",
    "                cumulative_load=[0.0]\n",
    "                discharge=[0.0]\n",
    "                for q,fr in zip(rep_discharge,annual_frequency):\n",
    "                    discharge.append(q)\n",
    "                    cumulative_load.append(cumulative_load[-1]+fr*10**(np.log10(q)*sand_pow[0]+sand_pow[1]))\n",
    "                percent_cumulative_load = [100*l/cumulative_load[-1] for l in cumulative_load]\n",
    "                \n",
    "                eff_discharge.set_xlabel('Discharge[cfs]')\n",
    "                eff_discharge.set_ylabel('Sand Load [tons/year]')\n",
    "                eff_discharge.plot(rep_discharge, annual_frequency*(10**(np.log10(rep_discharge)*sand_pow[0]+sand_pow[1])),color='tab:red')\n",
    "                eff_discharge.tick_params(axis='y', labelcolor='tab:red')\n",
    "                eff_discharge.grid()\n",
    "                pbar.update(11)\n",
    "                cum_discharge.set_ylabel('Cummulative Percent Sand Load', color='tab:blue')  # we already handled the x-label with ax1\n",
    "                cum_discharge.plot(discharge,percent_cumulative_load)\n",
    "                points = np.interp([20,80],percent_cumulative_load, discharge)\n",
    "                cum_discharge.plot([points[0]],[20],'ro')\n",
    "                cum_discharge.plot([points[1]],[80],'ro')\n",
    "                cum_discharge.tick_params(axis='y', labelcolor='tab:blue')\n",
    "                eff_discharge.set_title(\"Effective Discharge\")\n",
    "                \n",
    "                cumulative_load_log=[0.0]\n",
    "                discharge_log=[0.0]\n",
    "                for q,fr in zip(rep_discharge_log,annual_frequency_log):\n",
    "                    discharge_log.append(q)\n",
    "                    cumulative_load_log.append(cumulative_load_log[-1]+fr*10**(math.log10(q)*sand_pow[0]+sand_pow[1]))\n",
    "                percent_cumulative_load_log = [100*l/cumulative_load_log[-1] for l in cumulative_load_log]\n",
    "                \n",
    "                eff_discharge_log.set_xlabel('Discharge[cfs]')\n",
    "                eff_discharge_log.set_ylabel('Sand Load [tons/year]')\n",
    "                eff_discharge_log.plot(rep_discharge_log, annual_frequency_log*(10**(np.log10(rep_discharge_log)*sand_pow[0]+sand_pow[1])),color='tab:red')\n",
    "                eff_discharge_log.tick_params(axis='y', labelcolor='tab:red')\n",
    "                eff_discharge_log.grid()\n",
    "                pbar.update(12)\n",
    "                cum_discharge_log.set_ylabel('Cummulative Percent Sand Load', color='tab:blue')  # we already handled the x-label with ax1\n",
    "                cum_discharge_log.plot(discharge_log,percent_cumulative_load_log)\n",
    "                points = np.interp([20,80],percent_cumulative_load_log, discharge_log)\n",
    "                cum_discharge_log.plot([points[0]],[20],'ro')\n",
    "                cum_discharge_log.plot([points[1]],[80],'ro')\n",
    "                cum_discharge_log.tick_params(axis='y', labelcolor='tab:blue')\n",
    "                eff_discharge_log.set_title(\"Effective Discharge\")\n",
    "                pbar.update(13)\n",
    "            else:\n",
    "                pbar.update(13)\n",
    "    else:\n",
    "        #print(station_zipfile.filelist,station_zipfile)\n",
    "        fig.suptitle(\"INSUFFICIENT DATA AT {0}\".format(site[1]))\n",
    "        fig_hist.suptitle(\"INSUFFICIENT DATA AT {0}\".format(site[1]))\n",
    "        fig_rating_curve.suptitle(\"INSUFFICIENT DATA AT {0}\".format(site[1]))\n",
    "        fig_eff_discharge.suptitle(\"INSUFFICIENT DATA AT {0}\".format(site[1]))\n",
    "        tab.children[DISCRETE_TAB].children = []\n",
    "        tab.children[DAILY_TAB].children = []\n",
    "        pbar.reset(0)\n",
    "    fig.canvas.draw()\n",
    "    fig_hist.canvas.draw()\n",
    "    fig_rating_curve.canvas.draw()\n",
    "    fig_eff_discharge.canvas.draw()\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "stations = []\n",
    "site_dict={}\n",
    "\n",
    "def station_on_click_handler(event=None, id=None, properties=None, **args):\n",
    "    plot_station_daily(site_dict[tuple(args['coordinates'])])\n",
    "\n",
    "for site_no, station_nm, lat, long in zip(f.site_no,f.station_nm,f.dec_lat_va, f.dec_long_va):\n",
    "\n",
    "    stations.append(Marker(location=(float(lat),float(long)), \n",
    "                           draggable=False,\n",
    "                           title=\"Site: {0}\\nStation: {1}\".format(site_no,station_nm)))\n",
    "    site_dict[tuple(stations[-1].location)] = (site_no, station_nm)\n",
    "    stations[-1].on_click(station_on_click_handler)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "marker_cluster = MarkerCluster(\n",
    "    markers=tuple(stations)\n",
    ")\n",
    "\n",
    "m.add_layer(marker_cluster);"
   ]
  },
  {
   "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.8.1"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": false,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": false,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}