king_co_fish_passage_v1.1.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernel_info": {
      "name": "white-box"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "nteract": {
      "version": "0.15.0"
    },
    "colab": {
      "name": "Copy of Demo_maptool_v0_2.ipynb",
      "private_outputs": true,
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/cnilsen/kcfp/blob/master/king_co_fish_passage_v1.1.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7NsNYv2iySQW"
      },
      "source": [
        "#King County Fish Passage Hydrology Tool\n",
        "\n",
        "> Copyright (c) Geosyntec Consultants, 2019 <br>\n",
        "<div>The copyright holders grant the freedom to copy, modify, convey, adapt, and/or redistribute this work under the terms of the\n",
        "<a href=\"https://spdx.org/licenses/MPL-2.0.html\">Mozilla Public Licesnse 2.0</a>\n",
        "\n",
        "</div>\n",
        "\n",
        "### **Note** \n",
        "Access to this the cloud tools in this notebook is restricted to members of the \n",
        "Stormwater Heatmap Google Users Group. \n",
        "Visit: https://groups.google.com/forum/#!forum/stormwaterheatmap to request access. \n",
        "\n",
        "### **Overview**\n",
        "This tool supports King County's comprehensive inventory and assessment of assets' on watercourses to\n",
        "determine if the assets are barriers to fish passage. King County is performing the assessment using\n",
        "methodology outlined in the Washington Department of Fish and Wildlife's Fish Passage Inventory,\n",
        "Assessment, and Prioritization Manual (WDFW, 2019).\n",
        "The Level B assessment outlined in the WDFW\n",
        "2019 Manual provides for a hydraulic analysis to determine if a culvert meets the velocity and depth\n",
        "requirements for fish passage between low and high fish passage flows. Per the WDFW 2019 Manual:\n",
        "For an instream feature to be considered a non-barrier, it should not obstruct upstream migration\n",
        "at any time between the low and high fish passage flows at that location:\n",
        "\n",
        "* Low fish passage flow is the 95% exceedance flow; and\n",
        "* High fish passage flow is the 10% exceedance flow during the months of adult upstream migration.\n",
        "\n",
        "### **Data Sources**\n",
        "\n",
        "Data used in this tool has been provided from the Nature Conservancy's Stormwater Heatmap Project. Documentation on the Stormwater Heatmap is in progress. Below is a brief summary of data sources: \n",
        "\n",
        "* Digital Elevation Model - [USGS National Hydrography (NHD Plus High Resolution Raster Data)](https://www.usgs.gov/core-science-systems/ngp/national-hydrography/nhdplus-high-resolution)\n",
        "* Soils - [USDA Gridded SURGO Database](https://www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/home/?cid=nrcs142p2_053628)\n",
        "* Slope - USGS National Elevation Dataset (https://ned.usgs.gov)\n",
        "* Level B Sites - Spatial data provided by King County \n",
        "* Land cover and hydrology- The Nature Conservancy Stormwater Heatmap \n",
        "\n",
        "Documentation for the Stormwater Heatmap is in progress. Hydrology simulations were performed using [PyHSPF](https://github.com/djlampert/PyHSPF). Watershed factors used in this simulation were selected based on stormwater modeling guidance from the Department of Ecology. See the [Western Washington Hydrology Model Users Manual for more information.](https://fortress.wa.gov/ecy/madcap/wq/2014SWMMWWinteractive/2014%20SWMMWW.htm#Topics/VolumeIII2014/VolIII%20AppB%202014.htm%3FTocPath%3D2014%2520SWMMWW%7CVolume%2520III%2520-%2520Hydrologic%2520Analysis%2520and%2520Flow%2520Control%2520BMPs%7C)\n",
        "\n",
        "###**Limitations**\n",
        "\n",
        "\n",
        "**Results not calibrated**\n",
        "> Hydrology simulations were performed on a hydrologic response unit basins using regional parameters. Parameters have not been calibrated to actual data. A model-based gridded precipitation dataset was used (see [Mauger et al 2018](https://cig.uw.edu/news-and-events/publications/new-projections-of-changing-heavy-precipitation-in-king-county/) ). \n",
        "\n",
        "\n",
        "**Model does not include routing or storage**\n",
        ">Modeling was performed on a lumped parameter basis. Attenuation due to routing or storage is not reflected in the model. For large watersheds (greater than ~ 1 square mile) high flows may not be accurate. \n",
        "\n",
        "**Model does not include deep groundwater baseflow**\n",
        ">The modeling results include all hspf flow paths (surface flow, interflow, and groundwater flow), however baseflow from deep groundwater is not included. Consistent with WDFW guidance, results in this notebook have assumed a minimum baseflow of 1 cfs. \n",
        "---\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-_8Sv0oXCsOM"
      },
      "source": [
        "\n",
        "# Project set up and functions\n",
        "\n",
        "Execute the cells below to set up the notebook. \n",
        "If cell throws an error, you can safely ignore. \n",
        "\n",
        "If prompted to restart the runtime, click ```restart runtime```\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LHUDB7ROKkQe"
      },
      "source": [
        "#@title Update Folium Library { vertical-output: true, display-mode: \"form\" }\n",
        "#@markdown Please run this cell. Double click to show code.  \n",
        "\n",
        "# folium update\n",
        "!pip uninstall folium -y\n",
        "!pip install folium\n",
        "import folium\n",
        "!pip install geopandas"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ERVL7FtT_0xN"
      },
      "source": [
        "#@title Global Functions { run: \"auto\", display-mode: \"form\" }\n",
        "#@markdown Please run this cell. Double click to show code.  \n",
        "\n",
        "#@title Import Libraries\n",
        "import geopandas as gpd\n",
        "import pandas as pd\n",
        "import folium\n",
        "import branca\n",
        "import requests\n",
        "import json\n",
        "import math\n",
        "from itertools import product\n",
        "\n",
        "import numpy\n",
        "import altair\n",
        "import pandas_gbq\n",
        "from google.colab import widgets\n",
        "from datetime import datetime\n",
        "import ee\n",
        "from folium.features import GeoJson, GeoJsonTooltip\n",
        "\n",
        "DEFAULT_PROJECT = \"stormwaterheatmap-hydrology\"\n",
        "SQM_TO_ACRES = 4046.8564224\n",
        "_VALID_SOIL_TYPES = {0: 'A/B', 1: 'C', 2: 'D'}\n",
        "_VALID_LANDUSES = {0: 'forest', 1: 'pasture', 2: 'lawn', 5: 'impervious'}\n",
        "_VALID_SLOPE_CLASSES = {0: 'flat', 1: 'moderate', 2: 'steep'}\n",
        "ALL_HRUs = [\n",
        "    f'hru{soil}{landuse}{slope}'\n",
        "    for soil, landuse, slope in product(_VALID_SOIL_TYPES, _VALID_LANDUSES,\n",
        "                                        _VALID_SLOPE_CLASSES)\n",
        "    if (landuse != 5) or (landuse == 5 and soil == 2)\n",
        "]\n",
        "\n",
        "ALL_MONTHS = list(range(1, 13))\n",
        "\n",
        "def _get_soil(hruname):\n",
        "    stype = int(hruname[0])\n",
        "    return _VALID_SOIL_TYPES[stype]\n",
        "\n",
        "def _get_land_cover(hruname):\n",
        "    lutype = int(hruname[1])\n",
        "    return _VALID_LANDUSES[lutype]\n",
        "\n",
        "def _get_slope(hruname):\n",
        "    slopetype = int(hruname[2])\n",
        "    return _VALID_SLOPE_CLASSES[slopetype]\n",
        "\n",
        "# 'hru description': lambda df: df.apply(lambda row: (row['soil'])+\", \"+(row['land cover'])+\", \"+(row['slope']), axis=1),\n",
        "def _hru_descr(hruname):\n",
        "    return ', '.join([fxn(hruname) for fxn in (_get_soil, _get_land_cover, _get_slope)])\n",
        "\n",
        "def get_HRUs_in_watershed(ee, watershed_feature, tile_scale=10, image_file=None):\n",
        "    \"\"\"\n",
        "    Retrieve a list of HRUs within a watershed from\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    ee : Initialized Earth Engine object\n",
        "    watershed_geom : shapely Polygon\n",
        "        Shapely object representing the watershed\n",
        "    tile_ecale : int (optional, default = 10)\n",
        "        Scale to which the detail of the image file should be reduced\n",
        "    image_file : str (optional)\n",
        "        EE image name that contains the HRUs. When not provided, falls back to:\n",
        "        \"users/stormwaterheatmap/hruOut_fixed\".\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    HRUs : dict\n",
        "        Dictionary where the keys are the HRU and the values are the areas in\n",
        "        square meters\n",
        "\n",
        "    \"\"\"\n",
        "    if not image_file:\n",
        "        image_file = \"users/stormwaterheatmap/public/hrusJan2020Mode\"\n",
        "    \n",
        "    hrus = ee.Image(image_file)\n",
        "    area = ee.Image.pixelArea()\n",
        "    img = ee.Image.cat(area,hrus)\n",
        "    regionReduce = img.reduceRegion(\n",
        "        ee.Reducer.sum().group(1, 'hru'),\n",
        "        watershed_feature, 2\n",
        "    ).get('groups').getInfo()\n",
        "    df = pd.DataFrame(regionReduce).rename(columns={\"sum\": \"sq.m\"})\n",
        "    return df\n",
        "\n",
        "\n",
        "ALL_HRUs = [\n",
        "    f'hru{soil}{landuse}{slope}'\n",
        "    for soil, landuse, slope in product(_VALID_SOIL_TYPES, _VALID_LANDUSES,\n",
        "                                        _VALID_SLOPE_CLASSES)\n",
        "    if (landuse != 5) or (landuse == 5 and soil == 2)\n",
        "]\n",
        "ALL_MONTHS = list(range(1, 13))\n",
        "\n",
        "def _get_soil(hruname):\n",
        "    stype = int(hruname[0])\n",
        "    return _VALID_SOIL_TYPES[stype]\n",
        "\n",
        "\n",
        "def _get_land_cover(hruname):\n",
        "    lutype = int(hruname[1])\n",
        "    return _VALID_LANDUSES[lutype]\n",
        "\n",
        "\n",
        "def _get_slope(hruname):\n",
        "    slopetype = int(hruname[2])\n",
        "    return _VALID_SLOPE_CLASSES[slopetype]\n",
        "\n",
        "# 'hru description': lambda df: df.apply(lambda row: (row['soil'])+\", \"+(row['land cover'])+\", \"+(row['slope']), axis=1),\n",
        "def _hru_descr(hruname):\n",
        "    return ', '.join([fxn(hruname) for fxn in (_get_soil, _get_land_cover, _get_slope)])\n",
        "\n",
        "def get_HRUs_in_watershed(ee, watershed_id, scale = 10, tile_scale=10, image_file=None):\n",
        "    \"\"\"\n",
        "    Retrieve a list of HRUs within a watershed from\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    ee : Initialized Earth Engine object\n",
        "    watershed_geom : shapely Polygon\n",
        "        Shapely object representing the watershed\n",
        "    tile_ecale : int (optional, default = 10)\n",
        "        Scale to which the detail of the image file should be reduced\n",
        "    image_file : str (optional)\n",
        "        EE image name that contains the HRUs. When not provided, falls back to:\n",
        "        \"users/stormwaterheatmap/hruOut_fixed\".\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    HRUs : dict\n",
        "        Dictionary where the keys are the HRU and the values are the areas in\n",
        "        square meters\n",
        "\n",
        "    \"\"\"\n",
        "    if not image_file:\n",
        "        image_file = \"users/stormwaterheatmap/public/hrusJan2020Mode\"\n",
        "    \n",
        "    hrus = ee.Image(image_file)\n",
        "    area = ee.Image.pixelArea()\n",
        "    img = ee.Image.cat(area,hrus)\n",
        "    regionReduce = img.reduceRegion(\n",
        "        ee.Reducer.sum().group(1, 'hru'),\n",
        "        watershed, scale\n",
        "    ).get('groups').getInfo()\n",
        "    df = pd.DataFrame(regionReduce).rename(columns={\"sum\": \"sq.m\"})\n",
        "    return df\n",
        "\n",
        "\n",
        "def process_HRUs(df):\n",
        "    hru_mapper = dict(zip(\n",
        "        [40, 41, 42, 140, 141, 142, 240, 241, 242, 50, 51, 52, 150, 151, 152],\n",
        "        [250, 251, 252, 250, 251, 252, 250, 251, 252, 250, 251, 252, 250, 251, 252]\n",
        "    ))\n",
        "\n",
        "    hru_df = (\n",
        "        df.assign(hru=lambda df: df['hru'].replace(hru_mapper))\n",
        "        .assign(hruname=lambda df: df['hru'].map(lambda x: f'{x:03d}'))\n",
        "        .loc[lambda df: ~df['hruname'].str.contains('3')]\n",
        "        .groupby(by=['hru', 'hruname']).sum()\n",
        "        .reset_index()\n",
        "        .assign(**{\n",
        "            'acres': lambda df: df['sq.m'] / SQM_TO_ACRES,\n",
        "            'soil': lambda df: df['hruname'].apply(_get_soil),\n",
        "            'land cover': lambda df: df['hruname'].apply(_get_land_cover),\n",
        "            'slope': lambda df: df['hruname'].apply(_get_slope),\n",
        "            'hru description': lambda df: df['hruname'].apply(_hru_descr),\n",
        "        })\n",
        "    )\n",
        "    return hru_df\n",
        "\n",
        "\n",
        "\n",
        "def run_query(datadir=None, x=None, y=None):\n",
        "    #Authenticate to bigquery\n",
        "    #key_path = datadir / \"stormwaterheatmap-hydrology-8308154a1232.json\"\n",
        "    pandas_gbq.context.project = \"stormwaterheatmap-hydrology\"\n",
        "    pandas_gbq.context.credentials = service_account.Credentials.from_service_account_file(\n",
        "        str(key_path))\n",
        "\n",
        "    #Generate query statement\n",
        "    gridquery = \"\"\"\n",
        "        with location as (\n",
        "            SELECT\n",
        "            st_contains(st_geogfromtext(wkt_geom),\n",
        "                st_geogpoint(\n",
        "                    {lon},\n",
        "                    {lat}\n",
        "                    ) ) as grid,\n",
        "                    geohash\n",
        "            FROM\n",
        "            `stormwaterheatmap-hydrology.geometry.geohash` )\n",
        "            select geohash from location\n",
        "            where grid\n",
        "        \"\"\".format(lon=x, lat=y)\n",
        "        #run query to find precip grid geohash\n",
        "    project_id = \"stormwaterheatmap-hydrology\"\n",
        "    bqOutput = pd.read_gbq(gridquery, project_id=project_id)\n",
        "    gridID = bqOutput.values[0][0]\n",
        "    namesAreas = []\n",
        "    for hru in hrustoQuery.keys():\n",
        "        namesAreas.append(\n",
        "        hru + \"*\" + str(hrustoQuery.get(hru))\n",
        "        )\n",
        "    quantstmt = \"\"\"\n",
        "    approx_quantiles( ' + '\n",
        "    + '.join({namesAreas}) + ', '\n",
        "    +str(numofQuantiles)+ ')\n",
        "    AS Quantiles'\n",
        "    \"\"\".format(namesAreas=namesAreas)\n",
        "    #generate quantile query statement\n",
        "    qry = \"\"\"\n",
        "        Select\n",
        "        {selectstmt}\n",
        "        FROM   `stormwaterheatmap-hydrology.gfdl.{grid_id}`\n",
        "        WHERE  year BETWEEN {year0} AND    {yearN}\n",
        "        AND    month IN {months}\n",
        "        \"\"\".format(\n",
        "            grid_id=gridID, selectstmt=quantstmt, year0=year0, yearN=yearN, months=months)\n",
        "    #Run query\n",
        "    bqOutput = pd.read_gbq(qry, project_id=project_id)\n",
        "    return bqOutput\n",
        "\n",
        "\n",
        "def process_quantiles(flow_quantiles, name):\n",
        "    flows = numpy.array(flow_quantiles.Quantiles[0]) * 9.80962964e-6  # mm/hr*m2 to cfs\n",
        "    #flows = numpy.delete(flows, -1)  # remove the 0 value\n",
        "    n = (flows.shape[0]-1)\n",
        "    probs = numpy.arange(0, 1+(1 / n), 1 / n)\n",
        "    excd = (1 - probs).round(2)\n",
        "    return pd.DataFrame({'Exceedance': excd, name: flows})\n",
        "\n",
        "\n",
        "\n",
        "def get_geohash_of_point(centroid, project_id=None):\n",
        "    if not project_id:\n",
        "        project_id = DEFAULT_PROJECT\n",
        "\n",
        "    sql = '''with test as (\n",
        "    SELECT\n",
        "    st_contains(st_geogfromtext(wkt_geom),\n",
        "        st_geogpoint({lon},\n",
        "            {lat}) ) as grid,\n",
        "            geohash\n",
        "    FROM\n",
        "    `stormwaterheatmap-hydrology.geometry.geohash` )\n",
        "    select geohash from test\n",
        "    where grid\n",
        "    '''.format(lon=centroid[0], lat=centroid[1])\n",
        "    return pd.read_gbq(sql, project_id=project_id).geohash[0]\n",
        "\n",
        "\n",
        "def get_flow_quantiles(hru_df, geohash, year0, yearN, nquantiles,\n",
        "                       *months, project_id=None):\n",
        "    if not project_id:\n",
        "        project_id = DEFAULT_PROJECT\n",
        "\n",
        "    names_and_areas = hru_df[['hruname', 'sq.m']].apply(\n",
        "        lambda r: f\"(hru{r['hruname']} * {r['sq.m']})\", axis=1\n",
        "    ).tolist()\n",
        "\n",
        "    quantstmt = 'approx_quantiles(({names_areas}), {nquantiles}) AS Quantiles'.format(\n",
        "        names_areas=' + '.join(names_and_areas),\n",
        "        nquantiles=nquantiles\n",
        "    )\n",
        "\n",
        "    qry = \"\"\"\n",
        "    SELECT\n",
        "    {selectstmt}\n",
        "    FROM   `stormwaterheatmap-hydrology.gfdl.{grid_id}`\n",
        "    WHERE  year BETWEEN {year0} AND {yearN}\n",
        "    AND    month IN ({months})\n",
        "    \"\"\".format(\n",
        "        grid_id=geohash, selectstmt=quantstmt,\n",
        "        year0=year0, yearN=yearN,\n",
        "        months=(', '.join(str(x) for x in months))\n",
        "    )\n",
        "    return pd.read_gbq(qry, project_id=project_id)\n",
        "\n",
        "\n",
        "\n",
        "def prob_plot(quants):\n",
        "    ymax = math.ceil(quants['Qhigh'].max(0))\n",
        "    qplot = altair.Chart(quants).transform_fold(\n",
        "        ['Qlow', 'Qhigh']\n",
        "    ).mark_line().encode(\n",
        "        x=altair.X('Exceedance:Q', axis=altair.Axis(format='%')),\n",
        "        y=altair.Y(\n",
        "            'value:Q',\n",
        "            scale=altair.Scale(\n",
        "                type='log',\n",
        "                domain=[1,ymax],\n",
        "                clamp=True),\n",
        "                axis=altair.Axis(\n",
        "                    title = 'Discharge (cfs)'\n",
        "                )),\n",
        "        color='key:N',\n",
        "        tooltip=['Exceedance:Q', 'value:Q'], \n",
        "        )\n",
        "    return qplot\n",
        "\n",
        "\n",
        "def hilo_flow_widgets(quants):\n",
        "    low = (quants.loc[quants['Exceedance'] == 0.95]['Qlow']).values[0]\n",
        "    high =round((quants.loc[quants['Exceedance'] == 0.10]['Qhigh']).values[0],2)\n",
        "\n",
        "    low_label = widgets.Label(\n",
        "        value=f'Low Flow: {low:.3g} cfs',\n",
        "        description='Low Flow (cfs)',\n",
        "    )\n",
        "\n",
        "    high_label = widgets.Label(\n",
        "        value=f'High Flow: {round(high,3)} cfs',\n",
        "        description='Low Flow (cfs)',\n",
        "    )\n",
        "\n",
        "    return widgets.VBox([\n",
        "        widgets.Label('Results'),\n",
        "        low_label,\n",
        "        high_label\n",
        "    ])\n",
        "\n",
        "\n",
        "\n",
        "    #add ee method to folium\n",
        "    folium.Map.add_ee_layer = add_ee_layer\n",
        "\n",
        "    # Set visualization parameters.\n",
        "    LCvis_params = {'min': 0, 'max': 5,\n",
        "                    \"palette\":[\"55775e\",\"dacd7f\",\"7e9e87\",\"b3caff\",\"844c8b\",\"ead1ff\"]};\n",
        "    soils_vis = {'min': 0, 'max': 2,\n",
        "                    \"palette\":['#564138', '#69995D', '#F06543']}\n",
        "    slope_vis = {'min': 0, 'max': 2, 'dimensions': 400,\n",
        "                    \"palette\": ['#009B9E', '#F1F1F1', '#C75DAB']}\n",
        "\n",
        "    #Read watershed information from the watershed dictionary\n",
        "    #center = watershed_dict[\"pour_point\"]\n",
        "    #watershed_gdf = GeoData(geo_dataframe=(watershed_dict[\"watershed_geometry\"]), name='watershed')\n",
        "    watershed = table3.filter(ee.Filter.eq(\"OBJECTID\",2878))\n",
        "\n",
        "    # Create a folium map object.\n",
        "    eeMap = folium.Map(location=[center.y, center.x], zoom_start=13)\n",
        "\n",
        "    #download stormwaterheatmap data from earth engine\n",
        "\n",
        "    #landcover is a stormwaterheatmap asset:\n",
        "    landcover = ee.Image(\"users/stormwaterheatmap/landcover_2m\").clip(watershed)\n",
        "\n",
        "    #soils data is a stormwaterheatmap asset:\n",
        "    soils = ee.Image(\"users/stormwaterheatmap/soils2m\").clip(watershed)\n",
        "\n",
        "    #we generate slope information from the USGS National Elevation dataslet:\n",
        "    NED = ee.Image(\"USGS/NED\").clip(watershed)\n",
        "\n",
        "    #classify by WWHM thresholds:\n",
        "    thresholds = ee.Image([5.0, 15, 100]);\n",
        "    slopeZone = ee.Terrain.slope(NED).gt(thresholds).reduce('sum');\n",
        "    slope = slopeZone.clip(watershed)\n",
        "\n",
        "    #add to the folium Map\n",
        "    eeMap.add_ee_layer(ee, soils, soils_vis, 'Soils')\n",
        "    eeMap.add_ee_layer(ee, slope, slope_vis, 'Slope')\n",
        "    eeMap.add_ee_layer(ee, landcover, LCvis_params, 'Land cover')\n",
        "\n",
        "    # Add a layer control panel to the map.\n",
        "    eeMap.add_child(folium.LayerControl())\n",
        "    return eeMap\n",
        "\n",
        "def add_ee_layer(self, ee, ee_image_object, vis_params, name):\n",
        "    map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)\n",
        "    folium.raster_layers.TileLayer(\n",
        "        tiles=map_id_dict['tile_fetcher'].url_format,\n",
        "        attr=\n",
        "        'Map Data &copy; <a href=\"https://earthengine.google.com/\">Google Earth Engine</a>',\n",
        "        name=name,\n",
        "        overlay=True,\n",
        "        control=True).add_to(self)\n",
        "\n",
        "%matplotlib inline"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-VMwvXDyySQm"
      },
      "source": [
        "## Earth Engine Authentication\n",
        "\n",
        "This notebook requires a Google Earth Engine Account to query and analyze spatial data. Visit https://earthengine.google.com/signup to sign up for an account.  \n",
        "\n",
        "After you have registered your account. Run the `ee.Authenticate` function to authenticate your access to Earth Engine servers and `ee.Initialize` to initialize it. Upon running the following cell you'll be asked to grant Earth Engine access to your Google account. Follow the instructions printed to the cell. (Sometimes this takes awhile). "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ErYMTAMBySQo"
      },
      "source": [
        "ee.Authenticate() \n",
        "ee.Initialize()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fI8iVGXYKk-W"
      },
      "source": [
        "## Big Query Authentication \n",
        "Likewise, you will need to authenticate to Big Query by running the cell below. \n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sGWYXpMCKo0j"
      },
      "source": [
        "#@title  { vertical-output: true, display-mode: \"form\" }\n",
        "#@markdown Please run this cell. Double click to show code.  \n",
        "project_id = \"stormwaterheatmap-hydrology\"\n",
        "pandas_gbq.context.project = project_id\n",
        "#test authentication by running a sample query\n",
        "sql = \"\"\"\n",
        "SELECT geohash FROM `stormwaterheatmap-hydrology.geometry.geohash` LIMIT 1\n",
        "\"\"\"\n",
        "geohash_test_auth = pandas_gbq.read_gbq(\n",
        "    sql,\n",
        "    project_id=project_id\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sFzLqvOgL-_K"
      },
      "source": [
        "<style>\n",
        "div.blue { background-color:#e6f0ff; border-radius: 5px; padding: 16px 20px 8px 20px;}\n",
        "</style>\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "svF9ZQoCK2pe"
      },
      "source": [
        "# Enter simulation parameters\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nnSJiRujySQf"
      },
      "source": [
        "#@title Input for Calculating Flow Quantiles { run: \"auto\", vertical-output: true, form-width: \"500px\", display-mode: \"form\" }\n",
        "#@markdown Edit the variables below for your particular case:   \n",
        "from google.colab import widgets\n",
        "\n",
        "# Now we can create a grid, optional header_row and header_column\n",
        "# control whether we want header elements in the grid\n",
        "\n",
        "#@markdown **First Year of Simulation:** \n",
        "year0 = 1970 #@param {type:\"slider\", min:1970, max:2000, step:1}\n",
        "#@markdown **Last Year of Simulation:**  \n",
        "yearN = 2000  #@param {type:\"slider\", min:1970, max:2000, step:1}\n",
        "\n",
        "#low flow is set for all months \n",
        "low_flow_months = ALL_MONTHS \n",
        "#@markdown **Months to use for High Flow Migration Period:** *```hard coded for now```*\n",
        "#@markdown \n",
        "#@markdown ```January```\n",
        "\n",
        "#@markdown **Months to use for Low Flow Migration Period:** *```hard coded for now```*  \n",
        "#@markdown \n",
        "#@markdown ```January-December```\n",
        "\n",
        "#high flow is set for upstream migration  \n",
        "high_flow_months = list(range(1,2)) #high flow is calculated for January only \n",
        "#@markdown **Number of Quantile Values to Calculate:** \n",
        "\n",
        "numofQuantiles = 100 #@param {type:\"integer\"}"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MEnhnc5gySQ6"
      },
      "source": [
        "# Site Selection\n",
        "Next, genearate a map of the watersheds. Since some libraries are limited in colab, the way to select a site from the map is to find the watershed ID, ```AU_ID``` from the map and enter it as text. This method is being used to demonstrate the quantile calcuation functions. \n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0PHOqvBb1FtO"
      },
      "source": [
        "#@title Display and choose a watershed { run: \"auto\", form-width: \"400px\", display-mode: \"form\" }\n",
        "#@markdown This notebook can use any watershed feature class in  ```.geojson``` format. For now, use the geojson uploaded to github. \n",
        "\n",
        "#@markdown The default watersheds are from \n",
        "#@markdown Ecology's [Puget Sound Watershed Characterization Project](https://fortress.wa.gov/ecy/coastalatlas/wc/landingpage.html)\n",
        "#https://raw.githubusercontent.com/cnilsen/kcfp/master/01-data/king_co_au_single.geojson\n",
        "#@markdown **Enter URL below** (or leave as-is to accept the default): \n",
        "watershed_fileName = \"https://raw.githubusercontent.com/cnilsen/kcfp/master/01-data/king_co_au_single.geojson\" #@param {type:\"string\"}\n",
        "#@markdown Enter the name of the field to use as the watershed id\n",
        "idField = \"AU_ID\" #@param {type:\"string\"}\n",
        "#@markdown Enter the projection information (default is WGS84)\n",
        "crs = \"EPSG:4326\"#@param {type:\"string\"}\n",
        "response = requests.get(watershed_fileName)\n",
        "data = response.json()\n",
        "sheds = gpd.GeoDataFrame.from_features(\n",
        "    data, crs=crs)\n",
        "sheds[\"Colab_ID\"]= sheds[idField].astype(str)\n",
        "sheds = sheds.set_index(\"Colab_ID\",drop=False)\n",
        "\n",
        "# Create a folium map object.\n",
        "folium.Map(prefer_canvas=True)\n",
        "m = folium.Map(location=[47.3903, -122.0454], zoom_start=10, height=500)\n",
        "\n",
        "#add tooltip\n",
        "tooltip = GeoJsonTooltip(\n",
        "    fields=[\"Colab_ID\"],\n",
        "    aliases=[\"WatershedID\"],\n",
        "    sticky=False,\n",
        "    labels=True,\n",
        "    style=\"\"\"\n",
        "        background-color: #F0EFEF;\n",
        "        border: 2px solid black;\n",
        "        border-radius: 3px;\n",
        "        box-shadow: 3px;\n",
        "    \"\"\"\n",
        ")\n",
        "g = folium.GeoJson(\n",
        "    sheds,\n",
        "    tooltip = tooltip\n",
        ").add_to(m)\n",
        "\n",
        "m"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Qv-ZUIuG2HYF"
      },
      "source": [
        "#type in AU ID\n",
        "#@title Enter Selected Watershed AU ID { run: \"auto\", display-mode: \"form\" }\n",
        "#@markdown From the above map, enter the watershed ID (AU ID) in the field below: \n",
        "WatershedID = \t\"6505\"  #@param {type: \"string\"}\n",
        " \n",
        "aoi = sheds.loc[WatershedID,\"geometry\"]\n",
        "pd.DataFrame(sheds.loc[WatershedID,:]).drop(\"geometry\", axis=0)\n",
        "\n",
        "eeWS = gpd.GeoSeries([aoi]).__geo_interface__\n",
        "#list of coordinates to pass to earth engine\n",
        "eegeometry = list(eeWS['features'][0]['geometry']['coordinates'][0]) \n",
        "watershed = ee.Algorithms.GeometryConstructors.Polygon(eegeometry).getInfo()\n",
        "HRUs = get_HRUs_in_watershed(ee, watershed)\n",
        "\n",
        "def detail_map(watershed, ee=None):\n",
        "    #add ee method to folium\n",
        "    folium.Map.add_ee_layer = add_ee_layer\n",
        "\n",
        "    # Set visualization parameters.\n",
        "    LCvis_params = {'min': 0, 'max': 5,\n",
        "                    \"palette\":[\"55775e\",\"dacd7f\",\"7e9e87\",\"b3caff\",\"844c8b\",\"ead1ff\"]};\n",
        "    soils_vis = {'min': 0, 'max': 2,\n",
        "                    \"palette\":['#564138', '#69995D', '#F06543']}\n",
        "    slope_vis = {'min': 0, 'max': 2, 'dimensions': 400,\n",
        "                    \"palette\": ['#009B9E', '#F1F1F1', '#C75DAB']}\n",
        "\n",
        "    #Read watershed information from the watershed dictionary\n",
        "    center =  ee.Geometry(watershed).centroid()\n",
        "    \n",
        "\n",
        "    # Create a folium map object.\n",
        "    eeMap = folium.Map(location=[center.y, center.x], zoom_start=13)\n",
        "\n",
        "    #download stormwaterheatmap data from earth engine\n",
        "\n",
        "    #landcover is a stormwaterheatmap asset:\n",
        "    landcover = ee.Image(\"users/stormwaterheatmap/landcover_2m\").clip(watershed)\n",
        "\n",
        "    #soils data is a stormwaterheatmap asset:\n",
        "    soils = ee.Image(\"users/stormwaterheatmap/soils2m\").clip(watershed)\n",
        "\n",
        "    #we generate slope information from the USGS National Elevation dataslet:\n",
        "    NED = ee.Image(\"USGS/NED\").clip(watershed)\n",
        "\n",
        "    #classify by WWHM thresholds:\n",
        "    thresholds = ee.Image([5.0, 15, 100]);\n",
        "    slopeZone = ee.Terrain.slope(NED).gt(thresholds).reduce('sum');\n",
        "    slope = slopeZone.clip(watershed)\n",
        "\n",
        "    #add to the folium Map\n",
        "    eeMap.add_ee_layer(ee, soils, soils_vis, 'Soils')\n",
        "    eeMap.add_ee_layer(ee, slope, slope_vis, 'Slope')\n",
        "    eeMap.add_ee_layer(ee, landcover, LCvis_params, 'Land cover')\n",
        "\n",
        "    # Add a layer control panel to the map.\n",
        "    eeMap.add_child(folium.LayerControl())\n",
        "    return eeMap\n",
        "\n",
        "folium.Map.add_ee_layer = add_ee_layer\n",
        "\n",
        "# Set visualization parameters.\n",
        "LCvis_params = {'min': 0, 'max': 5,\n",
        "                    \"palette\":[\"55775e\",\"dacd7f\",\"7e9e87\",\"b3caff\",\"844c8b\",\"ead1ff\"]};\n",
        "soils_vis = {'min': 0, 'max': 2,\n",
        "                    \"palette\":['#564138', '#69995D', '#F06543']}\n",
        "slope_vis = {'min': 0, 'max': 2, 'dimensions': 400,\n",
        "                    \"palette\": ['#009B9E', '#F1F1F1', '#C75DAB']}\n",
        "\n",
        "    #Read watershed information from the watershed dictionary\n",
        "center_obj =  ee.Geometry(watershed).centroid().getInfo()\n",
        "center = center_obj.get('coordinates')\n",
        "# Create a folium map object.\n",
        "eeMap = folium.Map(location=[center[1], center[0]], zoom_start=13)\n",
        "\n",
        "    #download stormwaterheatmap data from earth engine\n",
        "\n",
        "    #landcover is a stormwaterheatmap asset:\n",
        "landcover = ee.Image(\"users/stormwaterheatmap/landcover_2m\").clip(watershed)\n",
        "\n",
        "    #soils data is a stormwaterheatmap asset:\n",
        "soils = ee.Image(\"users/stormwaterheatmap/soils2m\").clip(watershed)\n",
        "\n",
        "    #we generate slope information from the USGS National Elevation dataslet:\n",
        "NED = ee.Image(\"USGS/NED\").clip(watershed)\n",
        "\n",
        "#classify by WWHM thresholds:\n",
        "thresholds = ee.Image([5.0, 15, 100]);\n",
        "slopeZone = ee.Terrain.slope(NED).gt(thresholds).reduce('sum');\n",
        "slope = slopeZone.clip(watershed)\n",
        "\n",
        "    #add to the folium Map\n",
        "eeMap.add_ee_layer(ee, soils, soils_vis, 'Soils')\n",
        "eeMap.add_ee_layer(ee, slope, slope_vis, 'Slope')\n",
        "eeMap.add_ee_layer(ee, landcover, LCvis_params, 'Land cover')\n",
        "\n",
        "    # Add a layer control panel to the map.\n",
        "eeMap.add_child(folium.LayerControl())\n",
        "\n",
        "\n",
        "  #@title\n",
        "def hru_barchart(hru_df):\n",
        "    chart = altair.Chart(hru_df)\n",
        "    c1 = chart.mark_text().encode(\n",
        "        y=altair.Y('hru description', axis=None)\n",
        "\n",
        "    )\n",
        "    c4 = chart.mark_bar().encode(\n",
        "        x =altair.X('acres:Q'),\n",
        "        y=altair.Y('hru description:N', sort=altair.EncodingSortField(field='acres', order='descending')),\n",
        "        color = 'hru description:N'\n",
        "        #color='land use'\n",
        "    ).properties(width=300)\n",
        "    #desc = c1.encode(text='hru description:N').properties(title='descriptions')\n",
        "    #ac = c1.encode(text='acres:Q').properties(title='acres')\n",
        "    #tab = altair.hconcat(desc, ac)\n",
        "    return c4\n",
        "\n",
        "def process_HRUs(df):\n",
        "    hru_df = (\n",
        "        #df.assign(hru=lambda df: df['hru'].replace(hru_mapper))\n",
        "        df.assign(hruname=lambda df: df['hru'].map(lambda x: f'{x:03d}'))\n",
        "        #.loc[lambda df: ~df['hruname'].str.contains('3')]\n",
        "        .groupby(by=['hru', 'hruname']).sum()\n",
        "        .reset_index()\n",
        "        .assign(**{\n",
        "            'acres': lambda df: df['sq.m'] / SQM_TO_ACRES,\n",
        "            'soil': lambda df: df['hruname'].apply(_get_soil),\n",
        "            'land cover': lambda df: df['hruname'].apply(_get_land_cover),\n",
        "            'slope': lambda df: df['hruname'].apply(_get_slope),\n",
        "            'hru description': lambda df: df['hruname'].apply(_hru_descr),\n",
        "        })\n",
        "    )\n",
        "    return hru_df\n",
        "\n",
        "HRUs = HRUs.assign(hruname=lambda df: df['hru'].map(lambda x: f'{x:03d}'))\n",
        "#HRUs.loc[lambda df: ~df['hruname'].str.contains('3')]\n",
        "HRUs.drop(HRUs[HRUs.hruname.str[1] == \"9\"].index,inplace=True)\n",
        "HRUs.drop(HRUs[HRUs.hruname == \"255\"].index,inplace=True)\n",
        "HRUs= process_HRUs(HRUs)\n",
        "\n",
        "centroid = ee.Geometry(watershed).centroid().getInfo().get('coordinates')\n",
        "geohash = get_geohash_of_point(centroid, project_id=project_id)\n",
        "\n",
        "\n",
        "low_flow_quants = get_flow_quantiles(HRUs, geohash, year0, yearN, numofQuantiles, *low_flow_months)\n",
        "high_flow_quants = get_flow_quantiles(HRUs, geohash, year0, yearN, numofQuantiles, *high_flow_months)\n",
        "\n",
        "quants = process_quantiles(low_flow_quants, 'Qlow').merge(\n",
        "    process_quantiles(high_flow_quants, 'Qhigh'), on=['Exceedance']\n",
        ")\n",
        "\n",
        "qplot = prob_plot(quants)\n",
        "\n",
        "summary_table = quants.style.apply(\n",
        "    lambda x: [\"background: #1e88e5\" if x['Exceedance']==0.10  else(\"background: #F58518\" if x['Exceedance']==0.95 else \"\") for v in x], axis = 1).hide_index()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Smy4vowONDpV"
      },
      "source": [
        ""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1TaWYQauZLjs"
      },
      "source": [
        "# Display Results\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "glsdnBhQjNZf"
      },
      "source": [
        "#@title  { run: \"auto\",  form-width: \"400px\", display-mode: \"form\" }\n",
        "#rounding \n",
        "decimals = pd.Series([2, 2, 2], index = ['Exceedance','Qlow','Qhigh'])\n",
        "high = quants.round(decimals)[(quants.Exceedance == 0.1)][\"Qhigh\"] \n",
        "high = high.values[0]\n",
        "low = quants.round(decimals)[(quants.Exceedance == 0.95)][\"Qlow\"] \n",
        "low = low.values[0]\n",
        "lowflag = \"\"\n",
        "if low < 1.0:\n",
        "  low = 1.0 \n",
        "  lowflag = \"   Modeled low-flow below 1.0 cfs threshold.\"\n",
        "%load_ext google.colab.data_table\n",
        "from IPython.core.display import display, HTML\n",
        "\n",
        "html_report = f'''\n",
        "  <style type=\"text/css\">\n",
        "  div.test\n",
        "  {{\n",
        "  padding: 4px;\n",
        "  border-top-style: none;\n",
        "  border-right-style: none;\n",
        "  border-bottom-style: none;\n",
        "  border-left-style: solid;\n",
        "  border-color: #333333; \n",
        "  margin: 0px 0px 0px 50px;\n",
        "  }}\n",
        "  </style>\n",
        "  <div class=\"test\">\n",
        "  <b> High Passage Flow:</b> {high} cfs\n",
        "  <br> \n",
        "  <b>Low Passage Flow:</b> {low} cfs.\n",
        "  <br> \n",
        "  <br> <em>{lowflag}</em>\n",
        "     '''\n",
        "t = widgets.TabBar([\"Land use summary\",\"Hydrology Summary\",\"Map\"])\n",
        "\n",
        "## Land use summary - [0]\n",
        "with t.output_to(0):\n",
        "  grid = widgets.Grid(1, 2) \n",
        "\n",
        "  with grid.output_to(0,0):\n",
        "       display(hru_barchart(HRUs).properties(\n",
        "    title='Hydrologic Response Units'))\n",
        "       \n",
        "  with grid.output_to(0,1):\n",
        "       display(HRUs)\n",
        "\n",
        "## Hydrology Summary - [1]\n",
        "with t.output_to(1):\n",
        "  grid2 = widgets.Grid(1, 2) \n",
        "  with grid2.output_to(0,0):\n",
        "    display(qplot.mark_line() + qplot.mark_point().properties(\n",
        "    title='Flow Duration Curve'\n",
        "    ))\n",
        "    display(HTML(html_report))\n",
        "  with grid2.output_to(0,1):\n",
        "    display((quants.round(decimals)))\n",
        "    \n",
        "\n",
        "## Map - [2]\n",
        "with t.output_to(2):\n",
        "  display(eeMap)\n",
        "  \n",
        "\n",
        "#GRID --- \n",
        "\n",
        "  \n",
        "  "
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}