Landlab_landslide_probability_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example of Landlab 'landslide' component"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial provides an example driver code that runs Landlab's 'landslide' component.\n",
    "What a driver does is gathers input data, run the Landlab landslide component, and visualizes data and outputs.\n",
    "Landlab is a Python-based landscape modeling environment and the landslide component is one of many components available for users to access and link together to build their own landscape model.\n",
    "For more information about Landlab, see http://landlab.github.io/#/.\n",
    "\n",
    "Input data is provide by the user and consists of elevation from a DEM to provide topographic traits such as slope, contributing area, and flow direction. The user also supplies soil characteristics derived from a soil survey,\n",
    "land cover, or other sources, including transmissivity, cohesion, internal angle of friction, density, and thickness. \n",
    "Data for this example can be acquired from https://www.HydroShare.org under the resource \"Thunder Creek Landlab Landslide Example\".\n",
    "\n",
    "Method calculates factor-of-safety stability index by using node specific parameters, creating distributions of these parameters, and calculating the index by sampling these distributions 'n' times. The index is calculated from the 'infinite slope stabilty factor-of-safety equation' in the format of Pack RT, Tarboton DG, and Goodwin CN (1998) 'The SINMAP approach to terrain stability mapping.'\n",
    "\n",
    "Output includes figures of relative wetness, mean factor-of-safety, and probability of failure based on \n",
    "factor-of-safety calculations within a Monte Carlo simulation.\n",
    "\n",
    "This version allows users to provide just a minimum and maximum recharge that us used in a uniform distribution rather than a spatially distributed recharge field. Thus, it is for testing and teaching purposes. It is design to run from the\n",
    "same directory where the data files are located.\n",
    "\n",
    "@author: R.Strauch and E.Istanbulluoglu - Univerity of Washington Created on Thu Aug 20 16:47:11 2015 Last edit July 20, 2016\n",
    "\n",
    "Strauch, R., Istanbulluoglu, E., Nudurupati, S. S., Bandaragoda, C., Gasparini, N. M., and Tucker, G. E.: A hydroclimatological approach to predicting regional landslide probability using Landlab, Earth Surf. Dynam., 6, 49–75, https://doi.org/10.5194/esurf-6-49-2018, 2018. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To run example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To run this example, click in each shaded cell below and \"shift + enter\" to run each cell.\n",
    "Alternatively, you can run groups of cells by clicking \"Cell\" on the menu above and selecting you run option.  This is also where you can clear outputs from previous runs.\n",
    "\n",
    "If an error occurs, try \"Restart\" the kernel by clicking \"Kernel\" on the menu above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import libraries and components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "%matplotlib inline\n",
    "import pylab\n",
    "\n",
    "from landlab.plot.imshow import imshow_grid\n",
    "from landlab.plot.imshow import imshow_grid_at_node\n",
    "\n",
    "# read ESRI\n",
    "from landlab.io import read_esri_ascii\n",
    "\n",
    "#from landlab.components.flow_accum import flow_accumulation\n",
    "from landlab.components.flow_accum import find_drainage_area_and_discharge\n",
    "\n",
    "#from landlab.components.flow_accum import FlowAccumulator, SinkFiller\n",
    "from landlab.components import FlowAccumulator, SinkFiller\n",
    "\n",
    "\n",
    "from landlab.components.landslides import LandslideProbability\n",
    "\n",
    "import time\n",
    "st = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE - Warnings about matplotlib and gFlex are OKAY and won't affect running the component."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load grid, fields, and constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reading a esri_ascii file from ArcGIS sets up the RasterModelGrid and assign elevation field to nodes. Nodes are the center point of grid cells or pixels that are 30 m by 30 m in this example.\n",
    "The ascii (txt) files are located in the same directory where this notebook is run.\n",
    "This might take 30 seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid,Z)=read_esri_ascii('data/PNW_watershed_dem.asc', name='topographic__elevation')  # read the DEM again  #NFN_21m.asc\n",
    "\n",
    "grid.set_watershed_boundary_condition('topographic__elevation') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_width=grid.dx\n",
    "cell_width"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Zmin=np.min(grid.at_node['topographic__elevation'][grid.core_nodes])\n",
    "Zmax=np.max(grid.at_node['topographic__elevation'][grid.core_nodes])\n",
    "\n",
    "print(\"Min elevation of core nodes=\", Zmin) # print minimum elevation of core nodes\n",
    "print(\"Max elevation of core nodes=\", Zmax) # print maximum elevation of core nodes\n",
    "\n",
    "imshow_grid(grid,'topographic__elevation', plot_name = 'Topographic Elevation', \n",
    "            var_name = 'Elevation', var_units = 'm', grid_units = ('m','m'), \n",
    "            cmap = 'cividis', limits = (259, 853))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Flow accumulations and drainage are calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now fa is your object to run the FlowAccumulator method with the D8 flow directions and fill the depressions\n",
    "fa = FlowAccumulator(grid,\n",
    "                     surface='topographic__elevation',\n",
    "                     flow_director='FlowDirectorD8',\n",
    "                     runoff_rate=None,\n",
    "                     depression_finder='DepressionFinderAndRouter')\n",
    "# depression_finder='DepressionFinderAndRouter')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(da, q) = fa.accumulate_flow()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('max drainge area (km^2)=', np.max(grid.at_node['drainage_area'][grid.nodes])/1000000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Map contributing area of each cell**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow_grid(grid,'drainage_area', plot_name = 'Contributing Area', \n",
    "            var_name = 'Contributing Area', var_units = 'm^2', grid_units = ('m','m'), \n",
    "            cmap = 'viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow_grid(grid, np.log10(da), plot_name = 'Log(Catchment Area)', \n",
    "            var_name = 'Catchment Area', var_units = 'm^2', grid_units = ('m','m'), \n",
    "            cmap = 'viridis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Slope - Area plot of the watershed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "core_mask = np.where(grid.core_nodes)\n",
    "boundary_neighbor_mask = grid.node_has_boundary_neighbor(grid.nodes.flatten())\n",
    "core_interior = np.all([core_mask, np.invert(boundary_neighbor_mask)], 0)\n",
    "\n",
    "\n",
    "plt.loglog(\n",
    "    grid.at_node['drainage_area'][core_interior],\n",
    "    grid.at_node['topographic__steepest_slope'][core_interior],\n",
    "    'k.'\n",
    ")\n",
    "\n",
    "plt.xlabel('drainage area (m^2)')\n",
    "plt.ylabel('slope (m/m)')\n",
    "_ = plt.title('Slope-Area Plot')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see landscape concavity increases where draige area is greater than 1.5.10^4, we can identify channels using this area threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Channels=np.zeros(grid.number_of_nodes)\n",
    "_=grid.add_field('node','Channels',Channels,clobber=True) #field to how saturated grids\n",
    "\n",
    "grid.at_node['Channels'][np.where(da>15000)] = 1   # cont\n",
    "\n",
    "imshow_grid(grid, 'Channels', plot_name = 'Channels=1')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Prepare  Input for Landslide probability component**\n",
    "\n",
    "Check what other fields you need for the landslide probability component. Fields are data values that are assigned to each node, providing spatial variability to landscape. You check the input parameter names and input fields by commenting out the lines below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sorted(LandslideProbability.input_var_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#LandslideProbability?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check for Variable Units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('soil mode total cohesion :', LandslideProbability.var_units('soil__mode_total_cohesion'))\n",
    "print('saturated hydraulic conductivity :', LandslideProbability.var_units('soil__saturated_hydraulic_conductivity'))\n",
    "print('soil transmissivity :', LandslideProbability.var_units('soil__transmissivity'))\n",
    "print('soil internal friction angle :', LandslideProbability.var_units('soil__internal_friction_angle'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Specific catchment area* is calculate here and assigned as a field using the name searched by the component"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_width=grid.dx\n",
    "\n",
    "ca=grid.at_node['drainage_area']/cell_width  # now calculate SCA\n",
    "\n",
    "_=grid.add_field('node', 'topographic__specific_contributing_area', ca, clobber=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Topographic slope:* we will use the steepest slope calculated as part of flow routing as input to landslide probability using the name \"topographic__slope\" because landslide probability requires this name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow_grid(grid,'topographic__steepest_slope', plot_name = 'Local Slope', \n",
    "            var_name = 'Local Slope', var_units = 'fraction', grid_units = ('m','m'), \n",
    "            cmap = 'plasma', limits = (0, 1.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TS=grid.at_node['topographic__steepest_slope']\n",
    "\n",
    "_=grid.add_field('node', 'topographic__slope', TS, clobber=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Assign model parameters as fields**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# these three parameters are used to paramaterize a triangler distribution used to quantify uncertainty in root cohesion over space and time\n",
    "\n",
    "C_min=3     # min cohesion kPa-- this includes both soil and vegetation cohesion\n",
    "C_max=14    # max cohesion kPa-- this includes both soil and vegetation cohesion\n",
    "Cmode=7    # mode of cohesion kPA\n",
    "\n",
    "soil_density=2000 # kg/m3\n",
    "\n",
    "Phi=30      # internal friction angle in degrees\n",
    "Ks=5        # Ksat (m/day)\n",
    "\n",
    "hs_min=1.5  #soil thickness min\n",
    "hs_max=1.5  #soil thickness max\n",
    "\n",
    "grid.add_field('node', 'soil__minimum_total_cohesion', C_min*1000*np.ones(grid.number_of_nodes), clobber=True)\n",
    "grid.add_field('node', 'soil__maximum_total_cohesion', C_max*1000*np.ones(grid.number_of_nodes), clobber=True)\n",
    "grid.add_field('node', 'soil__mode_total_cohesion', Cmode*1000*np.ones(grid.number_of_nodes), clobber=True)\n",
    "grid.add_field('node', 'soil__density', soil_density*np.ones(grid.number_of_nodes), clobber=True)\n",
    "grid.add_field('node', 'soil__internal_friction_angle', Phi*np.ones(grid.number_of_nodes), clobber=True)\n",
    "grid.add_field('node', 'soil__saturated_hydraulic_conductivity', Ks*np.ones(grid.number_of_nodes), clobber=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Soil Thickness**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soil thickness is low on hillslopes and deeper in the valleys. Here we use a simple model to relate soil thickness to elevation to create a soil thicnkness map "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hs=hs_max-((Z-Zmin)/(Zmax-Zmin))*(hs_max-hs_min)\n",
    "_=grid.add_field('node', 'soil__thickness', hs*np.ones(grid.number_of_nodes), clobber=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Transmissivity:** this is depth integrated saturated hydraulic conductivity, Ks.  This will come from the SSURGU data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "To=Ks*hs\n",
    "_=grid.add_field('node', 'soil__transmissivity', To*np.ones(grid.number_of_nodes), clobber=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recharge input: This is a conceptual range of input for the statical distribution of recharge input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "groundwater__recharge_mean=20; #\n",
    "groundwater__recharge_standard_deviation=5; #"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Run landslide component to calculate probability of failure**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_iterations = 2500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LS_prob=LandslideProbability(grid,\n",
    "                         number_of_iterations=number_of_iterations,\n",
    "                         groundwater__recharge_distribution='lognormal',\n",
    "                         groundwater__recharge_mean=groundwater__recharge_mean,\n",
    "                         groundwater__recharge_standard_deviation=groundwater__recharge_standard_deviation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LS_prob.calculate_landslide_probability()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make plots and figures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot relative wetness with countours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.figure('Mean Relative Wetness')\n",
    "imshow_grid_at_node(grid, 'soil__mean_relative_wetness', cmap='YlGnBu',\n",
    "                 grid_units=('coordinates', 'coordinates'),\n",
    "                 shrink=0.75, var_name='Relative Wetness',\n",
    "                 var_units='no units')\n",
    "# this plots the contour lines\n",
    "#cs = pylab.contour(elev, extent=[0,42500, 0,32500], hold='on', colors='gray')\n",
    "#pylab.clabel(cs, inline=True, fmt='%1i', fontsize=10, manual=manual_locals)\n",
    "#imshow_grid_at_node(grid, overlay_watershed, color_for_closed='None',\n",
    "#                 allow_colorbar=False, cmap='spring')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot probability of failure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure('Probability of Failure')\n",
    "imshow_grid_at_node(grid, 'landslide__probability_of_failure', cmap='OrRd',\n",
    "                 grid_units=('coordinates', 'coordinates'), shrink=1,\n",
    "                 var_name='Probability of Failure',limits = (0, 1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This map shows the probability of failure based on then number of simulations where the factor-of-safety was <= 1.0 out all the simulations. The probability tends to be higer in the upper portions of the watershed where cohesion is less and soils are shallower, quickly saturating."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prob = grid.at_node['landslide__probability_of_failure']\n",
    "plt.style.use('ggplot')\n",
    "plt.hist(prob, bins=[0, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])\n",
    "plt.xlabel('Probability of Failure')\n",
    "plt.ylabel ('Number of Cells')\n",
    "plt.title('Landslide Probability of Failure')\n",
    "plt.ylim = 10000\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Landslide probability plotted with respect to local slope"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "core_mask = np.where(grid.core_nodes)\n",
    "boundary_neighbor_mask = grid.node_has_boundary_neighbor(grid.nodes.flatten())\n",
    "core_interior = np.all([core_mask, np.invert(boundary_neighbor_mask)], 0)\n",
    "\n",
    "\n",
    "plt.loglog(\n",
    "    grid.at_node['topographic__steepest_slope'][core_interior],\n",
    "    grid.at_node['landslide__probability_of_failure'][core_interior],\n",
    "    'k.'\n",
    ")\n",
    "\n",
    "plt.xlabel('local slope')\n",
    "plt.ylabel('Probability of landsliding')\n",
    "_ = plt.title('Slope-Area Plot')\n",
    "\n",
    "plt.xlim(0.2, 1.2)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}