NOCA_runPaper_LandlabLandslide.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Watershed Dynamics Model \n",
    "\n",
    "<img src=\"https://www.washington.edu/brand/files/2014/09/W-Logo_Purple_Hex.png\" style=\"float:right;width:300px;padding:20px\">   \n",
    "\n",
    "Here we explore recharge-driven shallow landsliding in the North Cascades National Park Complex (NOCA) from a HydroShare Observatory. \n",
    "\n",
    "This Jupyter Notebook runs the Landlab LandslideProbability component on a 30-m digital elevation model (DEM) for NOCA using recharge from Variable Infiltration Capacity (VIC) hydrologic model as described in the paper: <br />\n",
    "#### Strauch et al. 2018. A hydro-climatological approach to predicting regional landslide probability using Landlab. Earth Surface Dynamics, 6, 1-26. <br />\n",
    "\n",
    "This notebook performs the following functions to replicate the paper findings:<br >\n",
    "1) Import libraries and set HydroShare variables<br />\n",
    "2) Review data needed as input for the landslide model<br />\n",
    "3) Create a RasterModelGrid based on a 30-m DEM<br />\n",
    "4) Access and assign data fields used to calculate landslide probability<br />\n",
    "5) Set Number of iterations to run Monte Carlo simulation<br />\n",
    "6) Specify recharge option as _data driven spatial_ and access Python dictionaries to generate recharge distributions<br /> \n",
    "7) Run Landlab LandslideProbability component<br /> \n",
    "8) Display and visualize results of stability analysis<br /> \n",
    "<br /> \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## To run this notebook:\n",
    "\n",
    "To run this example, click in each shaded cell below and \"shift + enter\" to run each cell. 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": [
    "## 1.  HydroShare Setup and Preparation\n",
    "\n",
    "To run this notebook, we must import several libraries.\n",
    "The hs_utils library provides functions for interacting with HydroShare, including resource querying, dowloading. and creation. Additional libraries support the functions of Landlab."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import Python utilities for calculating and plotting\n",
    "import six\n",
    "import os\n",
    "import matplotlib as mpl\n",
    "mpl.use('agg')\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import cPickle as pickle\n",
    "\n",
    "#import utilities for importing and exporting to HydroShare\n",
    "from utilities import hydroshare\n",
    "# set variables for interacting with HydroShare from this notebook\n",
    "hs=hydroshare.hydroshare()\n",
    "# Create object to map the home directory\n",
    "homedir = r'/home/jovyan/work/notebooks/data/' + str(os.environ[\"HS_RES_ID\"]) + '/' + str(os.environ[\"HS_RES_ID\"]) + '/data/contents/'\n",
    "print homedir\n",
    "\n",
    "# Import Landlab libraries\n",
    "import landslide_probability\n",
    "from landslide_probability import LandslideProbability\n",
    "from landlab import imshow_grid_at_node\n",
    "# from landlab.plot.imshow import imshow_node_grid\n",
    "from landlab.io import read_esri_ascii\n",
    "from landlab.io import write_esri_ascii\n",
    "from collections import defaultdict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are curious about where the data is being downloaded, click on the Jupyter Notebook dashboard icon in upper rigth corner to see a File System view.  The homedir directory location printed above is where you can find the data and contents you will download to a HydroShare JupyterHub server.  At the end of this work session, you can migrate this data to the HydroShare iRods server as a Generic Resource. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import data\n",
    " This data was preprocessed for the North Cascades National Park Complex case study and is on HydroShare as [Regional landslide hazard using Landlab - NOCA Data](https://www.hydroshare.org/resource/a5b52c0e1493401a815f4e77b09d352b/). Click on the link to see the data repository on HydroShare."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set a variable name for the data resource using the HydroShare resource ID."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Data_ResourceID='a5b52c0e1493401a815f4e77b09d352b'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Execute the next cell to download data from another HydroShare resource - this may take a minute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hs.getResourceFromHydroShare(Data_ResourceID)\n",
    "data_folder = '/home/jovyan/work/notebooks/data/'+ Data_ResourceID +'/'+Data_ResourceID+'/data/contents/Data files/'\n",
    "print('This is the location on the HydroShare JupyterHub server where the data has just been downloaded:')\n",
    "print data_folder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ## 2. Review data needed as input for the landslide model\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the list of data inputs that the component needs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sorted(LandslideProbability.input_var_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Review the details of what each variable represents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LandslideProbability._var_doc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the units of each variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LandslideProbability._var_units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will establish a RasterModelGrid based on a DEM for assigning our variables to.\n",
    "Nodes are the center point of grid cells or pixels that are 30 m by 30 m in this example. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Create a RasterModelGrid based on a 30-m DEM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load  landslide model inputs from ASCII textfile (ArcGIS raster conversion) into Landlab grid\n",
    "\n",
    "Load DEM elevation grid<br />\n",
    "Set_nodata_nodes_to_inactive that have no data (e.g., -9999), which establishes boundary conditions<br />\n",
    "This might take a few minutes as the park is large (2,757 km2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid, z) = read_esri_ascii(data_folder+'/elevation.txt',name='topographic__elevation')\n",
    "grid.at_node.keys()     # loads DEM grid with elevation\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['topographic__elevation'], -9999) # set boundary conditions closed where no data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Confirm the size of the grid, nodes located every 30 m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.number_of_nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will attach data to this grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Access and assign data fields to grid used to calculate landslide probability and set boundary conditions\n",
    "\n",
    "#### For each input below\n",
    "1. Load data from ascii text file\n",
    "2. Assign grid to Landlab node\n",
    "3. Set boundary conditions\n",
    "\n",
    "For the eSurf paper with NOCA extent, this takes ~60 sec to load each file on NCSA ROGER super computer. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load slope\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, slope) = read_esri_ascii(data_folder+'/slope_tang17d.txt')\n",
    "grid.add_field('node', 'topographic__slope', slope)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['topographic__slope'], -9999)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['topographic__slope'], 0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load contributing area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, ca) = read_esri_ascii(data_folder+'/cont_area.txt')\n",
    "grid.add_field('node', 'topographic__specific_contributing_area', ca)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['topographic__specific_contributing_area'], -9999)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load transmissivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, T) = read_esri_ascii(data_folder+'/transmis.txt')\n",
    "grid.add_field('node', 'soil__transmissivity', T)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['soil__transmissivity'], -9999)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load cohesion (mode, min, and max) - this takes ~3 minutes because 3 cohesion fields are provided to create more flexibility in how cohesion is distributed on the landscape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, C) = read_esri_ascii(data_folder+'/cohesion_mode.txt')\n",
    "C[C == 0.0] = 1.0  # ensure minimum is >0 Pa for use in distributions generation\n",
    "grid.add_field('node', 'soil__mode_total_cohesion', C)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['soil__mode_total_cohesion'], -9999)\n",
    "\n",
    "(grid1, C_min) = read_esri_ascii(data_folder+'/cohesion_min.txt')\n",
    "grid.add_field('node', 'soil__minimum_total_cohesion', C_min)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['soil__minimum_total_cohesion'], -9999)\n",
    "\n",
    "(grid1, C_max) = read_esri_ascii(data_folder+'/cohesion_max.txt')\n",
    "grid.add_field('node', 'soil__maximum_total_cohesion', C_max)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['soil__maximum_total_cohesion'], -9999)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load internal angle of friction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, phi) = read_esri_ascii(data_folder+'/frict_angle.txt')\n",
    "grid.add_field('node', 'soil__internal_friction_angle', phi)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['soil__internal_friction_angle'], -9999)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set soil density value and assign to all nodes as a constant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid['node']['soil__density'] = 2000*np.ones(grid.number_of_nodes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load soil thickness or depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, hs) = read_esri_ascii(data_folder+'/soil_depth.txt')\n",
    "grid.add_field('node', 'soil__thickness', hs)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['soil__thickness'], -9999)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load analysis mask and actual landslide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, mask) = read_esri_ascii(data_folder+'/exclud_mask.txt')\n",
    "grid.add_field('node', 'exclusion_mask', mask)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['exclusion_mask'], -9999)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load observed landslide inventory. Class 1-5 are landslides, 8 is no landslide mapped for later plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, slides) = read_esri_ascii(data_folder+'/landslide_type.txt')\n",
    "grid.add_field('node', 'landslides', slides)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Specify recharge option as _data driven spatial_ and access Python dictionaries to generate recharge distributions "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recharge is this model represents the annual maximum recharge in mm/day. This corresponds to the wettest conditions expected annually, which is the severest soil-saturated conditions likely to occur at least once a year. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify distribution to use in simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distribution = 'data_driven_spatial'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load pre-processed routed flows dictionaries containing HSD_id and fractional drainage at each node and recharge dictionaries.  HSD is the Hydrologic Source Domain, which is the VIC data in this case study at ~5x6 km2 grid size.  The 'pickle' utility loads existing dictionaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dict of node id (key) and HSD_ids (values)\n",
    "HSD_id_dict = pickle.load(open(data_folder+'/dict_uniq_ids.p', 'rb'))\n",
    "# dict of node id (key) and fractions (values)\n",
    "fract_dict = pickle.load(open(data_folder+'/dict_coeff.p', 'rb'))\n",
    "# dict of HSD id (key) with arrays of recharge (values)\n",
    "HSD_dict = pickle.load(open(data_folder+'/HSD_dict.p', 'rb'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combine dictionaries into __ordered__ parameters required for _data driven spatial_ distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HSD_inputs = [HSD_dict,HSD_id_dict, fract_dict]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Set Number of iterations to run Monte Carlo simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The landslide component employes the infinite slope model to calculate factor-of-safety index values using a Monte Carlo simulation, which randomly selects input values from parameter distributions. You can specify the number of iterations to run Monte Carlo simulation, but the default is 250. The higher the number of iteration, the longer the program runs, but the more precise the probability of failure results become."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iterations = 3000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Run the Landlab LandslideProbability Component"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To run the landslide model, we first instantiate the LandslideProbability component with the above parameters, as well as the grid and number of iterations we specified before. Instantiate creates an instance of a class. (For example, the iPhone is a class and each phone is an instance.)\n",
    "\n",
    "No outputs are generated by this command as it is setting up the recharge and instantiating the component."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LS_prob = LandslideProbability(grid,\n",
    "    number_of_iterations=iterations,\n",
    "    groudwater__recharge_distribution=distribution,\n",
    "    groudwater__recharge_HSD_inputs=HSD_inputs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the component has been instantiated, we generate outputs from running the component by calling the component's 'calculate_landslide_probability' method using the class instance (e.g., LS_prob). The following cell runs the model; in the following section we will assessing the results. These calculations will take a few minutes given the size of the modeling domain represented by core nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LS_prob.calculate_landslide_probability()\n",
    "print('Landslide probability successfully calculated')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The outputs of landslide model simulation are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sorted(LS_prob.output_var_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the arrays as variables by 'attaching the fields to the grid' and view the outputs. **component already does this**\n",
    "\n",
    "This simulation generates a probability value for each core node. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LS_prob_probability_of_failure = grid.at_node['landslide__probability_of_failure']\n",
    "grid.at_node['landslide__probability_of_failure']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This simulation generates a mean relative wetness value for each core node as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LS_prob_relative_wetness = grid.at_node['soil__mean_relative_wetness']\n",
    "grid.at_node['soil__mean_relative_wetness']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Display and visualize results of stability analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set plotting parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpl.rcParams['xtick.labelsize'] = 15\n",
    "mpl.rcParams['ytick.labelsize'] = 15\n",
    "mpl.rcParams['lines.linewidth'] = 1\n",
    "mpl.rcParams['axes.labelsize'] = 18\n",
    "mpl.rcParams['legend.fontsize'] = 15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot elevation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure('Elevations from the DEM [m]')\n",
    "imshow_grid_at_node(grid, 'topographic__elevation', cmap='terrain',\n",
    "                 grid_units=('coordinates', 'coordinates'),\n",
    "                 shrink=0.75, var_name='Elevation', var_units='m')\n",
    "#plt.savefig('NOCA_elevation.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Excluded areas from the analysis are shown in black, including outside the park and inside the park areas that are water bodies, snow, glaciers, wetlands, exposed bedrock, and slopes <= 17 degrees. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot slope overlaid with mapped landslide types. Takes about a few minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure('Landslides')\n",
    "ls_mask1 = grid.at_node['landslides'] != 1.0\n",
    "ls_mask2 = grid.at_node['landslides'] != 2.0\n",
    "ls_mask3 = grid.at_node['landslides'] != 3.0\n",
    "ls_mask4 = grid.at_node['landslides'] != 4.0\n",
    "overlay_landslide1 = np.ma.array(grid.at_node['landslides'], mask=ls_mask1)\n",
    "overlay_landslide2 = np.ma.array(grid.at_node['landslides'], mask=ls_mask2)\n",
    "overlay_landslide3 = np.ma.array(grid.at_node['landslides'], mask=ls_mask3)\n",
    "overlay_landslide4 = np.ma.array(grid.at_node['landslides'], mask=ls_mask4)\n",
    "imshow_grid_at_node(grid, 'topographic__slope', cmap='pink',\n",
    "                 grid_units=('coordinates', 'coordinates'), vmax=2.,\n",
    "                 shrink=0.75, var_name='Slope', var_units='m/m')\n",
    "imshow_grid_at_node(grid, overlay_landslide1, color_for_closed='None',\n",
    "                 allow_colorbar=False, cmap='cool')\n",
    "imshow_grid_at_node(grid, overlay_landslide2, color_for_closed='None',\n",
    "                 allow_colorbar=False, cmap='autumn')\n",
    "imshow_grid_at_node(grid, overlay_landslide3, color_for_closed='None',\n",
    "                 allow_colorbar=False, cmap='winter')\n",
    "imshow_grid_at_node(grid, overlay_landslide4, color_for_closed='None',\n",
    "                 allow_colorbar=False,cmap='summer')\n",
    "#plt.savefig('NOCA_Landslides_on_Slope.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Legend to mapped landslides: blue - debris avalanches, cyan - falls/topples, red - debris torrents, and green - slumps/creeps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot of soil depth (m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure('Soil Thickness')\n",
    "imshow_grid_at_node(grid, 'soil__thickness', cmap='copper_r',\n",
    "                 grid_units=('coordinates', 'coordinates'), shrink=0.75,\n",
    "                 var_name='Soil Thickness', var_units='m')\n",
    "#plt.savefig('NOCA_SoilDepth.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot probability of saturation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure('Probability of Saturation')\n",
    "imshow_grid_at_node(grid, 'soil__probability_of_saturation', cmap='YlGnBu',\n",
    "                 limits=((0), (1)),\n",
    "                 grid_units=('coordinates', 'coordinates'),\n",
    "                 shrink=0.75, var_name='Probability of Saturation',\n",
    "                 var_units='no units')\n",
    "#plt.savefig('NOCA_ProbabilityofSaturation.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This map shows the probability of saturation as high throughout much of the area because we modeled the annual maximum recharge, which is esssentially the worst case conditions that might lead to instability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot probability of failure; Compare this with the elevation and slope maps."
   ]
  },
  {
   "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=0.75,\n",
    "                 var_name='Probability of Failure', var_units='no units')\n",
    "#plt.savefig('NOCA_ProbabilityofFailure.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The map of probability of failure shows higher probabilities at higher elevations below retreating glaciers where vegation is sparse and shallow unconsolidated sediment is prevalent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To review the fields assigned to the grid, simply execute the following command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.at_node"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Export data from model run: FS probability, mean Reletive wetness, probability of saturation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "core_nodes = grid.core_nodes\n",
    "data_extracted = {'PF_3S_SD': np.array(\n",
    "                 grid.at_node['landslide__probability_of_failure'][grid.core_nodes]),\n",
    "                 'mean_RW': np.array(grid.at_node['soil__mean_relative_wetness']\n",
    "                 [grid.core_nodes]),'prob_sat': np.array(\n",
    "                 grid.at_node['soil__probability_of_saturation'][grid.core_nodes])}\n",
    "headers = ['PF_3S_SD','mean_RW','prob_sat']\n",
    "df = pd.DataFrame(data_extracted, index=core_nodes, columns=(headers))\n",
    "df.to_csv('FS3k_SSD_demo.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make ascii files for raster creation in GIS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "write_esri_ascii('prbF_3kSSD_demo.txt',grid,names='landslide__probability_of_failure')\n",
    "write_esri_ascii('mRW_3kSSD_demo.txt',grid,names='soil__mean_relative_wetness')\n",
    "write_esri_ascii('prbSat_3kSSD_demo.txt',grid,names='soil__probability_of_saturation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2.7",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}