2. WNTR geospatial tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# WNTR Geospatial Tutorial\n",
    "The following tutorial illustrates the use of the `wntr.gis` module to use geospatial data in resilience analysis.  The objective of this tutorial is to 1) quantify water service disruptions that could occur from pipes damaged in landslides and 2) identify the social vulnerability of populations impacted by the service disruptions.\n",
    "\n",
    "To simplify the tutorials, it is assumed that pipes within a 1000 ft buffer of each landslide susceptible region are damaged in that landslide.\n",
    "This assumption could be replaced with detailed landslide analysis that includes slope, soil type, weather conditions, and pipe material.\n",
    "Social vulnerability data could also be replaced by datasets that describe other attributes of the population and critical services."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Imports\n",
    "Install and import WNTR and additional Python packages that are needed for the tutorial\n",
    "- Geopandas is used to load geospatial data\n",
    "- Shapely is used to define a region of interest to crop data\n",
    "- Matplotlib is used to create subplots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install geopandas==0.14.4 shapely matplotlib wntr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "from shapely.geometry import box\n",
    "import matplotlib.pylab as plt\n",
    "import wntr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create directories to hold results if they do not already exist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "if not os.path.isdir('geojsons'):\n",
    "    os.makedirs('geojsons')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Units\n",
    "WNTR uses SI (International System) units (length in meters, time in seconds, mass in kilograms), **with the exception of the landslide buffer which is in feet to match the coordinate reference system of the geospatial data**.  See https://usepa.github.io/WNTR/units.html for more details on WNTR units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following line defines coordinates used to zoom in on network graphics\n",
    "zoom_coords = [(5.75e6, 5.79e6), (3.82e6, 3.85e6)] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Water Network Model\n",
    "The following section creates a `WaterNetworkModel` object from an EPANET INP file and converts the model to GeoDataFrames for use in geospatial analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a WaterNetworkModel from an EPANET INP file\n",
    "The water distribution network model used in this tutorial was downloaded from the [UKnowledge Water Distribution Systems Research Database](https://uknowledge.uky.edu/wdsrd/). KY10 was selected for the analysis. The following section creates a `WaterNetworkModel` from an EPANET INP file and computes some general attributes of the model.\n",
    "\n",
    "*Citation: Hoagland, Steven, \"10 KY 10\" (2016). Kentucky Dataset. 12. https://uknowledge.uky.edu/wdst/12. Accessed on 4/4/2024.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a WaterNetworkModel from an EPANET INP file\n",
    "inp_file = 'networks/ky10.inp'\n",
    "wn = wntr.network.WaterNetworkModel(inp_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print a basic description of the model. \n",
    "# The level can be 0, 1, or 2 and defines the level of detail included in the description.\n",
    "wn.describe(level=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute total pipe length\n",
    "length = wn.query_link_attribute('length')\n",
    "total_length = length.sum() # m\n",
    "print('Total pipe length =', total_length, 'm, =', total_length*3.28084, 'ft')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute average expected demand per day \n",
    "average_expected_demand = wntr.metrics.average_expected_demand(wn) # m^3/s\n",
    "average_volume_per_day = average_expected_demand*(24*3600) # m^3\n",
    "total_water_use = average_volume_per_day.sum() # m^3\n",
    "print('Total water use =', total_water_use, 'm^3, =', total_water_use*264.172/1e6, 'million gallons')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate population using the default average volume of water consumed per capita per day of 200 gallons/day\n",
    "population = wntr.metrics.population(wn) \n",
    "total_population = population.sum()\n",
    "print('Total population =', total_population)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create a basic network graphic, showing junction elevation\n",
    "# Note, the remaining graphics in this tutorial are created from the geospatial data directly, rather than the `plot_network` function.\n",
    "# The `plot_network` function currently does not include vertices.\n",
    "ax = wntr.graphics.plot_network(wn, node_attribute='elevation', node_range=(175, 300), title='ky10 elevation')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert the WaterNetworkModel to GIS data\n",
    "The `WaterNetworkModel` is converted to a collection of GIS compatible GeoDataFrames and the coordinate reference system (CRS) is set to **EPSG:3089 (NAD83 / Kentucky Single Zone (ftUS)**, see https://epsg.io/3089 for more details).  Data for junctions, tanks, reservoirs, pipes, pumps, and valves are stored in separate GeoDataFrames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Convert the WaterNetworkModel to GIS data and set the CRS\n",
    "wn_gis = wn.to_gis()\n",
    "wn_gis.junctions.head()\n",
    "#wn_gis.tanks.head()\n",
    "#wn_gis.reservoirs.head()\n",
    "#wn_gis.pipes.head()\n",
    "#wn_gis.pumps.head()\n",
    "#wn_gis.tanks.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Set the CRS to EPSG:3089 (NAD83 / Kentucky Single Zone (ftUS))\n",
    "crs = 'EPSG:3089'\n",
    "wn_gis.set_crs(crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the GIS data to create a figure of the network\n",
    "fig, ax = plt.subplots(figsize=(5,5))\n",
    "ax = wn_gis.pipes.plot(column='diameter', linewidth=1, label='pipes', alpha=0.8, ax=ax, zorder=1)\n",
    "ax = wn_gis.reservoirs.plot(color='k', marker='s', markersize=60, label='reservoirs', ax=ax)\n",
    "ax = wn_gis.tanks.plot(color='r', markersize=20, label='tanks', ax=ax)\n",
    "ax = wn_gis.pumps.centroid.plot(color='b', markersize=20, label='pumps', ax=ax)\n",
    "ax = wn_gis.valves.centroid.plot(color='c', markersize=20, label='valves', ax=ax)\n",
    "tmp = ax.axis('off')\n",
    "# Comment/uncomment the following 2 lines to change the zoom on the network graphic\n",
    "#tmp = ax.set_xlim(zoom_coords[0])\n",
    "#tmp = ax.set_ylim(zoom_coords[1])\n",
    "tmp = plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save the GIS data to GeoJSON or Shapefile files\n",
    "The GIS data can be written to GeoJSON files or Shapefile files.  One file is created for junctions, tanks, reservoirs, pipes, pumps, and valves (ky10_junctions.geojson, ky10_tanks.geojson, etc.).  The GeoJSON or Shapefile files can be loaded into GIS software platforms for further analysis. **Note that controls, patterns, curves, and options are not included in the GIS formatted data files.** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Store the WaterNetworkModel as a collection of GeoJSON files\n",
    "wn_gis.write_geojson('geojsons/ky10')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# External GIS Data\n",
    "The external data used in this tutorial includes landslide inventory and social vulnerability data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load landslide GIS data\n",
    "The landslide data used in this tutorial was downloaded from the [UKnowledge Kentucky Geological Survey Research Data](https://uknowledge.uky.edu/kgs_data/).  The Kentucky Geological Survey Landslide Inventory from March 2023 was selected for the analysis.  The data contains locations of known landslides and areas susceptible to debris flows, derived from aerial photography. \n",
    "\n",
    "*Citation: Crawford, M.M., 2023. Kentucky Geological Survey landslide inventory [2023-03]: Kentucky Geological Survey Research Data, https://uknowledge.uky.edu/kgs_data/7/, Accessed on 4/4/2024.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# To reduce the file size checked into the WNTR repository, the following code was run on the raw data file.\n",
    "# The region of interest (ROI) includes a 5000 ft buffer surrounding all pipes. The function `box` was imported from shapely.\n",
    "\"\"\"\n",
    "bounds = wn_gis.pipes.total_bounds # total_bounds returns the upper and lower bounds on x and y\n",
    "geom = box(*bounds)\n",
    "ROI = geom.buffer(5000) # feet\n",
    "\n",
    "landslide_file = '../data/KGS_Landslide_Inventory_exp.gdb'\n",
    "landslide_data = gpd.read_file(landslide_file, driver=\"FileGDB\", layer='Areas_derived_from_aerial_photography')\n",
    "print(landslide_data.crs)\n",
    "landslide_data = landslide_data.clip(ROI)\n",
    "landslide_data.to_file(\"../data/ky10_landslide_data.geojson\", index=True, driver='GeoJSON')\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Load the landslide data from file and print the CRS to ensure it is in EPSG:3089.  \n",
    "# The methods `to_crs` and `set_crs` can be used to change coordinate reference systems if needed.\n",
    "landslide_file = 'data/ky10_landslide_data.geojson'\n",
    "landslide_data = gpd.read_file(landslide_file).set_index('index') \n",
    "print(landslide_data.crs)\n",
    "\n",
    "landslide_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the landslide data along with pipes\n",
    "ax = landslide_data.plot(color='red', label='Landslide data')\n",
    "ax = wn_gis.pipes.plot(color='black', linewidth=1, ax=ax)\n",
    "ax.set_title('Landslide and pipe data')\n",
    "tmp = ax.axis('off')\n",
    "# Comment/uncomment the following 2 lines to change the zoom on the network graphic\n",
    "tmp = ax.set_xlim(zoom_coords[0])\n",
    "tmp = ax.set_ylim(zoom_coords[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Social Vulnerability Index (SVI) GIS data\n",
    "The social vulnerability data used in this tutorial was downloaded from the [Centers for Disease Control and Prevention/Agency for Toxic Substances and Disease Registry](https://www.atsdr.cdc.gov/placeandhealth/svi/index.html). The data contains census and social vulnerability metrics for each census tract. \n",
    "\n",
    "The quantity of interest used in this analysis is `RPL_THEMES` which ranks vulnerability across socioeconomic status, household characteristics, racial and ethnic minority status, and housing type and transportation.  The value ranges between 0 and 1, where higher values are associated with higher vulnerability.\n",
    "\n",
    "*Citation: Centers for Disease Control and Prevention/Agency for Toxic Substances and Disease Registry/Geospatial Research, Analysis, and Services Program. CDC/ATSDR Social Vulnerability Index 2020 Database Kentucky. https://www.atsdr.cdc.gov/placeandhealth/svi/data_documentation_download.html. Accessed on 4/4/2024.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# To reduce the file size checked into the WNTR repository, the following code was run on the raw data file. \n",
    "# The region of interest (ROI) was defined above.\n",
    "\"\"\"\n",
    "svi_file = '../data/SVI2020_KENTUCKY_tract.gdb'\n",
    "svi_data = gpd.read_file(svi_file, driver=\"FileGDB\", layer='SVI2020_KENTUCKY_tract')\n",
    "print(svi_data.crs)\n",
    "svi_data.to_crs(crs, inplace=True)\n",
    "svi_data = svi_data.clip(ROI)\n",
    "svi_data.to_file(\"../data/ky10_svi_data.geojson\", index=True, driver='GeoJSON')\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Load the SVI data from file and print the CRS to ensure it is in EPSG:3089.  \n",
    "# The methods `to_crs` and `set_crs` can be used to change coordinate reference systems if needed.\n",
    "svi_file = 'data/ky10_svi_data.geojson'\n",
    "svi_data = gpd.read_file(svi_file).set_index('index') \n",
    "print(svi_data.crs)\n",
    "\n",
    "svi_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot SVI data and pipes (higher values of SVI are associated with higher vulnerability)\n",
    "ax = svi_data.plot(column='RPL_THEMES', label='SVI data', cmap='RdYlGn_r', vmin=0, vmax=1, legend=True)\n",
    "ax = wn_gis.pipes.plot(color='black', linewidth=1, ax=ax)\n",
    "ax.set_title('SVI and pipe data')\n",
    "tmp = ax.axis('off')\n",
    "# Comment/uncomment the following 2 lines to change the zoom on the network graphic\n",
    "tmp = ax.set_xlim(zoom_coords[0])\n",
    "tmp = ax.set_ylim(zoom_coords[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expand the size of each landslide using a buffer\n",
    "Each landslide is extended to include the surrounding 1000 ft, to create a region that might be impacted by an individual landslide.  The distance unit for buffering matches the distance unit of the CRS (ft).\n",
    "This assumption could be replaced with detailed landslide analysis that includes slope, soil type, weather conditions, and pipe material. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create a GeoDataFrame to hold information used in landslide scenarios (initially copied from landslide_data)\n",
    "# Buffer each landslide polygon by 1000 ft\n",
    "landslide_scenarios = landslide_data.copy()\n",
    "landslide_scenarios['geometry'] = landslide_data.buffer(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add a prefix to the landslide scenario index to indicate the scenario name\n",
    "landslide_scenarios.index = 'LS-' + landslide_scenarios.index.astype(str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot the landslide data, region included in each landslide scenario, and pipes\n",
    "ax = landslide_scenarios.plot(color='gray', alpha=0.5)\n",
    "ax = landslide_data.plot(color='red', label='Landslide data', ax=ax)\n",
    "ax = wn_gis.pipes.plot(color='black', linewidth=1, ax=ax)\n",
    "ax.set_title('Landslide scenario and pipe data')\n",
    "tmp = ax.axis('off')\n",
    "# Comment/uncomment the following 2 lines to change the zoom on the network graphic\n",
    "tmp = ax.set_xlim(zoom_coords[0])\n",
    "tmp = ax.set_ylim(zoom_coords[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geospatial Intersects\n",
    "In this section, landslide scenario and SVI data are interested with pipes and junctions in the `WaterNetworkModel`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Identify pipes that intersect each landslide\n",
    "Landslide polygons are intersected with pipes to obtain a list of pipes that intersect each landslide.  <font color='red'>This information is used to to define the pipes that are closed in each landslide scenario.</font> The pipe attribute `length` is also included in the intersection to gather statistics on the pipe length that intersects each landslide.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Use the intersect function to determine pipes and pipe length that intersects each landslide\n",
    "A = landslide_scenarios\n",
    "B = wn_gis.pipes\n",
    "B_value = 'length'\n",
    "landslide_intersect = wntr.gis.intersect(A, B, B_value)\n",
    "\n",
    "# Print results in order of descending total pipe length\n",
    "landslide_intersect.sort_values('sum', ascending=False).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the intersection results to the landslide scenario data\n",
    "landslide_scenarios[['intersections', 'n', 'total pipe length']] = landslide_intersect[['intersections', 'n', 'sum']]\n",
    "\n",
    "# Print results in order of descending total pipe length\n",
    "landslide_scenarios.sort_values('total pipe length', ascending=False).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot intersection results\n",
    "fig, axes = plt.subplots(1,2, figsize=(15,5))\n",
    "\n",
    "wn_gis.pipes.plot(color='gray', linewidth=1, ax=axes[0])\n",
    "landslide_scenarios.plot(column='n', vmax=10, legend=True, ax=axes[0])\n",
    "tmp = axes[0].set_title('Number of pipes that intersect each landslide')\n",
    "tmp = axes[0].axis('off')\n",
    "# Comment/uncomment the following 2 lines to change the zoom on the network graphic\n",
    "tmp = axes[0].set_xlim(zoom_coords[0])\n",
    "tmp = axes[0].set_ylim(zoom_coords[1])\n",
    "\n",
    "wn_gis.pipes.plot(color='gray', linewidth=1, ax=axes[1])\n",
    "landslide_scenarios.plot(column='total pipe length', vmax=10000, legend=True, ax=axes[1])\n",
    "tmp = axes[1].set_title('Length of pipe that intersect each landslide')\n",
    "tmp = axes[1].axis('off')\n",
    "# Comment/uncomment the following 2 lines to change the zoom on the network graphic\n",
    "tmp = axes[1].set_xlim(zoom_coords[0])\n",
    "tmp = axes[1].set_ylim(zoom_coords[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Identify landslides that intersect each pipe\n",
    "Pipes are intersected with landslides to obtain a list of landslides that intersect each pipe. The landslide attribute `Confidence_Ranking` is also included in the intersection to gather statistics on landslide confidence for each pipe.  <font color='red'>While this information is not used in the analysis below, this type of information could be used to inform uncertainty or probability of damage.</font>\n",
    "\n",
    "**Note that `Confidence_Ranking` has a value of 3 (\"Landslide likely at or near the specified location\") for each landslide in region of interest. Since the values are uniform in this dataset, the intersected sum, min, max, and mean are all the same value.** More information on Confidence ranking can be found at https://kgs.uky.edu/kgsmap/helpfiles/landslide_help.shtm.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Use the intersect function to determine landslides and landslide confidence ranking that intersects each pipe\n",
    "A = wn_gis.pipes\n",
    "B = landslide_scenarios\n",
    "B_value = 'Confidence_Ranking'\n",
    "pipe_intersect = wntr.gis.intersect(A, B, B_value)\n",
    "\n",
    "# Print results in order of descending number of intersections.\n",
    "pipe_intersect.sort_values('n', ascending=False).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the intersection results to the GIS pipe data\n",
    "wn_gis.pipes[['intersections', 'n', 'Confidence_Ranking']] = pipe_intersect[['intersections', 'n', 'mean']]\n",
    "\n",
    "# Print results in order of descending number of intersections\n",
    "wn_gis.pipes.sort_values('n', ascending=False).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot intersection results\n",
    "ax = wn_gis.pipes.plot(color='gray', linewidth=1, zorder=1)\n",
    "wn_gis.pipes[wn_gis.pipes['n'] > 0].plot(column='n', vmax=20, legend=True, ax=ax)\n",
    "tmp = ax.set_title('Number of landslide scenarios that intersect each pipe')\n",
    "tmp = ax.axis('off')\n",
    "# Comment/uncomment the following 2 lines to change the zoom on the network graphic\n",
    "tmp = ax.set_xlim(zoom_coords[0])\n",
    "tmp = ax.set_ylim(zoom_coords[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assign social vulnerability to each junction\n",
    "Junctions are intersected with SVI to determine the social vulnerability of the population at each junction.  The SVI data attribute `RPL_THEMES` is included in the intersection. <font color='red'>This information is used to determine the social vulnerability of individuals that experience water service disruptions.</font>\n",
    "\n",
    "The SVI data column `RPL_THEMES` ranks vulnerability across socioeconomic status, household characteristics, racial and ethnic minority status, and housing type and transportation. The value ranges between 0 and 1, where higher values are associated with higher vulnerability.\n",
    "\n",
    "**Note that since each junction only intersects one census tract, the SVI sum, min, max, and mean are all the same value.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Use the intersect function to determine SVI of each junction.  \n",
    "A = wn_gis.junctions\n",
    "B = svi_data\n",
    "B_value = 'RPL_THEMES'\n",
    "junction_svi = wntr.gis.intersect(A, B, B_value)\n",
    "\n",
    "junction_svi.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the intersection results (SVI value) to the GIS junction data\n",
    "wn_gis.junctions['RPL_THEMES'] = junction_svi['mean']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot SVI for each census tract and SVI assigned to each junction\n",
    "fig, axes = plt.subplots(1,2, figsize=(15,5))\n",
    "\n",
    "svi_data.plot(column='RPL_THEMES', label='SVI data', vmin=0, vmax=1, legend=True, ax=axes[0])\n",
    "wn_gis.pipes.plot(color='black', linewidth=1, ax=axes[0])\n",
    "tmp = axes[0].set_title('Census tract SVI and pipe data')\n",
    "tmp = axes[0].axis('off')\n",
    "# Comment/uncomment the following 2 lines to change the zoom on the network graphic\n",
    "tmp = axes[0].set_xlim(zoom_coords[0])\n",
    "tmp = axes[0].set_ylim(zoom_coords[1])\n",
    "\n",
    "wn_gis.pipes.plot(color='gray', linewidth=1, ax=axes[1])\n",
    "wn_gis.junctions.plot(column='RPL_THEMES', vmin=0, vmax=1, legend=True, ax=axes[1])\n",
    "tmp = axes[1].set_title('SVI value at each junction')\n",
    "tmp = axes[1].axis('off')\n",
    "# Comment/uncomment the following 2 lines to change the zoom on the network graphic\n",
    "tmp = axes[1].set_xlim(zoom_coords[0])\n",
    "tmp = axes[1].set_ylim(zoom_coords[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Hydraulic Simulations\n",
    "The following section runs hydraulic simulations for the baseline (no landslide) and landslide scenarios. A subset of landslide scenarios is run to simply the tutorial.  Simulation results are stored for later analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create a function to setup the WaterNetworkModel for hydraulic simulations\n",
    "def model_setup(inp_file):\n",
    "    wn = wntr.network.WaterNetworkModel(inp_file)\n",
    "    wn.options.hydraulic.demand_model = 'PDD'\n",
    "    wn.options.hydraulic.required_pressure = 20 # m\n",
    "    wn.options.hydraulic.minimum_pressure  = 0 # m\n",
    "    wn.options.time.duration = 48*3600 # s (48 hour simulation)\n",
    "    return wn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run baseline simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Run a baseline simulation, with no landslide or damage.  \n",
    "wn = model_setup(inp_file)\n",
    "sim = wntr.sim.EpanetSimulator(wn)\n",
    "baseline_results = sim.run_sim()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View a subset of the simulation results\n",
    "baseline_results.node['pressure'].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Run landslide scenarios\n",
    "Landslide scenarios are downselected by identifying the set of landslides that impact a unique set of pipes.  Scenarios are further downselected to 6 scenarios to simplify the tutorial.  A hydraulic simulation is run for each landslide scenario, where pipes that intersect the landslide are closed for 48 hours.  Results from each scenario are stored for later analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Remove scenarios with no intersecting pipes\n",
    "landslide_scenarios = landslide_scenarios[landslide_scenarios['n'] > 0]\n",
    "landslide_scenarios = landslide_scenarios[~landslide_scenarios['n'].isna()]\n",
    "\n",
    "# Downselect landslide scenarios that impact a unique set of pipes\n",
    "duplicated_intersections = landslide_scenarios['intersections'].astype(str).duplicated()\n",
    "landslide_scenarios = landslide_scenarios.loc[~duplicated_intersections, :]\n",
    "\n",
    "print('Number of unique landslide scenarios', landslide_scenarios.shape[0])\n",
    "landslide_scenarios.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Further downselect the landslide scenarios for demonstration purposes. Choose one of the following 4 options.\n",
    "# Option 1. 6 scenarios that illustrate a wide range of impact\n",
    "landslide_scenarios_downselect = landslide_scenarios.loc[['LS-4495', 'LS-7003', 'LS-7111', 'LS-5086', 'LS-6966', 'LS-7058'],:] \n",
    "\n",
    "# Option 2. 6 scenarios with the highest intersecting pipe length\n",
    "#landslide_scenarios_downselect = landslide_scenarios.sort_values('total pipe length', ascending=False).iloc[0:6,:]\n",
    "\n",
    "# Option 3. 6 scenarios with the highest number of intersecting pipes\n",
    "#landslide_scenarios_downselect = landslide_scenarios.sort_values('n', ascending=False).iloc[0:6,:]\n",
    "\n",
    "# Option 4. Random selection of 6 scenarios\n",
    "#landslide_scenarios_downselect = landslide_scenarios.sample(n=6, random_state=1234)\n",
    "\n",
    "landslide_scenarios_downselect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the location of landslides used in the analysis\n",
    "ax = landslide_scenarios_downselect.plot(color='blue')\n",
    "wn_gis.pipes.plot(color='gray', linewidth=1, ax=ax)\n",
    "tmp = ax.set_title('Landslide scenarios')\n",
    "tmp = ax.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run a hydraulic simulation for each landslide scenario, store results in a dictionary\n",
    "# Each scenario closes all pipes that intersect the landslide for the 48 hour simulation\n",
    "results = {}\n",
    "for i, scenario in landslide_scenarios_downselect.iterrows():\n",
    "    wn = model_setup(inp_file)\n",
    "    for pipe_i in scenario['intersections']:\n",
    "        pipe_object = wn.get_link(pipe_i)\n",
    "        pipe_object.initial_status = 'CLOSED'\n",
    "    sim = wntr.sim.EpanetSimulator(wn)\n",
    "    results[i] = sim.run_sim()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Analysis Results\n",
    "The following section computes and plots analysis results, including water service availability (WSA) and the social vulnerability index (SVI) of impacted junctions for each scenario."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Water Service Availability (WSA)\n",
    "Water service availability (WSA) is the ratio of delivered demand to the expected demand.  WSA is computed for each junction (alternatively, WSA can be computed for each timestep, or for each junction and timestep).  A value below 1 indicates that expected demand it me, while a value of 0 indicates that the expected demand is not met. \n",
    "\n",
    "**Note that WSA can be > 1 and < 0 due to numerical differences in expected and actual demand. For certain types of analysis, WSA should be truncated to values between 0 and 1.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute expected demand for each junction and timestep\n",
    "expected_demand = wntr.metrics.expected_demand(wn)\n",
    "\n",
    "expected_demand.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute total expected demand at each junction (axis 0 is the time index)\n",
    "expected_demand_j = expected_demand.sum(axis=0)\n",
    "\n",
    "expected_demand_j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute baseline WSA for each junction\n",
    "demand_baseline = baseline_results.node['demand'].loc[:,wn.junction_name_list]\n",
    "demand_baseline_j = demand_baseline.sum(axis=0) # total demand at each junction\n",
    "wsa_baseline_j = wntr.metrics.water_service_availability(expected_demand_j, demand_baseline_j)\n",
    "\n",
    "wsa_baseline_j.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Add WSA from the base simulation to the junction GIS data\n",
    "wn_gis.junctions['baseline'] = wsa_baseline_j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot WSA from the base simulation\n",
    "ax = wn_gis.pipes.plot(color='black', linewidth=1)\n",
    "ax = wn_gis.junctions.plot(column='baseline', cmap='RdYlGn', vmin=0, vmax=1, legend=True, ax=ax)\n",
    "tmp = ax.set_title('Baseline WSA')\n",
    "tmp = ax.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute WSA associated with each landslide scenarios\n",
    "for i, scenario in landslide_scenarios_downselect.iterrows():\n",
    "    demand = results[i].node['demand'].loc[:,wn.junction_name_list]\n",
    "    demand_j = demand.sum(axis=0) # total demand at each junction\n",
    "    wsa_j = wntr.metrics.water_service_availability(expected_demand_j, demand_j)\n",
    "    \n",
    "    # Add WSA to the junction GIS data\n",
    "    wn_gis.junctions[i] = wsa_j\n",
    "    print(i, len(scenario['intersections']), wsa_j.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Extract WSA for each scenario \n",
    "wsa_results = wn_gis.junctions[landslide_scenarios_downselect.index]\n",
    "\n",
    "wsa_results.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot distribution of WSA for each scenario. Note that WSA can be > 1 and < 0 due to numerical differences in expected and actual demand. \n",
    "# For certain types of analysis, the WSA should be truncated to values between 0 and 1.\n",
    "ax = wsa_results.boxplot()\n",
    "tmp = ax.set_ylim(-0.25, 1.25)\n",
    "tmp = ax.set_ylabel('WSA')\n",
    "tmp = ax.set_title('Distribution of WSA for each scenario')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot WSA for each scenario\n",
    "fig, axes = plt.subplots(2,3, figsize=(15,10))\n",
    "axes = axes.flatten()\n",
    "\n",
    "for i, scenario in enumerate(wsa_results.columns):\n",
    "    wn_gis.pipes.plot(color='gray', linewidth=1, ax=axes[i]) # pipes\n",
    "    wn_gis.junctions.plot(column=scenario, cmap='RdYlGn', vmin=0, vmax=1, legend=True, ax=axes[i]) # junction wsa\n",
    "    tmp = axes[i].set_title('WSA '+scenario)\n",
    "    tmp = axes[i].axis('off')\n",
    "    if i >= 6: # axes is defined to have 6 subplots\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SVI of impacted junctions\n",
    "In this analysis, impacted junctions are defined as junctions where WSA falls below 0.5 (50% of the expected water was received) at any time during the simulation. Other criteria could also be used to defined impact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract junctions that are impacted by WSA < 0.5 for each scenario\n",
    "impacted_junctions = {}\n",
    "for scenario in wsa_results.columns:\n",
    "    filter = wsa_results[scenario] < 0.5\n",
    "    impacted_junctions[scenario] = wsa_results.index[filter]\n",
    "\n",
    "impacted_junctions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot SVI of impacted junctions for each scenario\n",
    "fig, axes = plt.subplots(2,3, figsize=(15,10))\n",
    "axes = axes.flatten()\n",
    "\n",
    "for i, scenario in enumerate(wsa_results.columns):\n",
    "    j = impacted_junctions[scenario]\n",
    "    wn_gis.pipes.plot(color='gray', linewidth=1, alpha=0.5, ax=axes[i]) # pipes\n",
    "    if len(j) > 0:\n",
    "        wn_gis.junctions.loc[j,:].plot(column='RPL_THEMES', cmap='RdYlGn_r', vmin=0, vmax=1, legend=True, ax=axes[i]) # junction wsa\n",
    "    tmp = axes[i].set_title('SVI '+scenario)\n",
    "    tmp = axes[i].axis('off')\n",
    "    if i >= 6: # axes is defined to have 6 subplots\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save analysis results to GIS files\n",
    "The analysis above stored WSA results for each scenario to the `wn_gis` object, which can be saved to GIS formatted files and loaded into GIS software platforms for further analysis.  **Note that lists (such as the information stored in 'intersections') is not JSON serializable and must first be removed.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del wn_gis.pipes['intersections']\n",
    "wn_gis.write_geojson('geojsons/ky10_analysis_results')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  },
  "toc-autonumbering": true,
  "vscode": {
   "interpreter": {
    "hash": "45e5fc2add644b17ca9d771f06d4d0791b1baf8fc0ad9b13d75afaf563ce7ce4"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}