SCL_load_spatial_parameters_DataDriven-testing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Seattle City Light Landslide Model \n",
    "\n",
    "Here we explore groundwater-driven shallow landsliding in the Skagit River Watershed above the confluence with Sauk River. \n",
    "\n",
    "This Jupyter Notebook runs the Landlab LandslideProbability component on a 30-m digital elevation model (DEM) using depth to groundwater from the Distributed Hydrology Soil and Vegetation (DHSVM) hydrologic model.\n",
    "\n",
    "This notebook performs the following functions:<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 goundwater option as _data driven spatial_ and access NetCDF to generate groundwater depth distributions<br /> \n",
    "7) Run Landlab LandslideProbability component<br /> \n",
    "8) Adjust root cohession to simulate fire and run Landlab LandslideProbability Component<br />\n",
    "9) 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": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ignore deprecation warnings\n",
      "Successfully established a connection with HydroShare\n",
      "Data will be loaded from and saved to:/home/jovyan/work/notebooks/data/4cac25933f6448409cab97b293129b4f/4cac25933f6448409cab97b293129b4f/data/contents\n",
      "/home/jovyan/work/notebooks/data/4cac25933f6448409cab97b293129b4f/4cac25933f6448409cab97b293129b4f/data/contents\n",
      "You are not going crazy...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/envs/python2/lib/python2.7/site-packages/landlab/__init__.py:27: UserWarning: \n",
      "This call to matplotlib.use() has no effect because the backend has already\n",
      "been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot,\n",
      "or matplotlib.backends is imported for the first time.\n",
      "\n",
      "The backend was *originally* set to 'module://ipykernel.pylab.backend_inline' by the following code:\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/runpy.py\", line 174, in _run_module_as_main\n",
      "    \"__main__\", fname, loader, pkg_name)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/runpy.py\", line 72, in _run_code\n",
      "    exec code in run_globals\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/ipykernel_launcher.py\", line 16, in <module>\n",
      "    app.launch_new_instance()\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/traitlets/config/application.py\", line 658, in launch_instance\n",
      "    app.start()\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/ipykernel/kernelapp.py\", line 486, in start\n",
      "    self.io_loop.start()\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/tornado/ioloop.py\", line 1064, in start\n",
      "    handler_func(fd_obj, events)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/tornado/stack_context.py\", line 276, in null_wrapper\n",
      "    return fn(*args, **kwargs)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/zmq/eventloop/zmqstream.py\", line 450, in _handle_events\n",
      "    self._handle_recv()\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/zmq/eventloop/zmqstream.py\", line 480, in _handle_recv\n",
      "    self._run_callback(callback, msg)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/zmq/eventloop/zmqstream.py\", line 432, in _run_callback\n",
      "    callback(*args, **kwargs)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/tornado/stack_context.py\", line 276, in null_wrapper\n",
      "    return fn(*args, **kwargs)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/ipykernel/kernelbase.py\", line 283, in dispatcher\n",
      "    return self.dispatch_shell(stream, msg)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/ipykernel/kernelbase.py\", line 233, in dispatch_shell\n",
      "    handler(stream, idents, msg)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/ipykernel/kernelbase.py\", line 399, in execute_request\n",
      "    user_expressions, allow_stdin)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/ipykernel/ipkernel.py\", line 208, in do_execute\n",
      "    res = shell.run_cell(code, store_history=store_history, silent=silent)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/ipykernel/zmqshell.py\", line 537, in run_cell\n",
      "    return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/IPython/core/interactiveshell.py\", line 2714, in run_cell\n",
      "    interactivity=interactivity, compiler=compiler, result=result)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/IPython/core/interactiveshell.py\", line 2818, in run_ast_nodes\n",
      "    if self.run_code(code, result):\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/IPython/core/interactiveshell.py\", line 2878, in run_code\n",
      "    exec(code_obj, self.user_global_ns, self.user_ns)\n",
      "  File \"<ipython-input-3-68bb3205dfc9>\", line 7, in <module>\n",
      "    get_ipython().magic(u'matplotlib inline')\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/IPython/core/interactiveshell.py\", line 2160, in magic\n",
      "    return self.run_line_magic(magic_name, magic_arg_s)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/IPython/core/interactiveshell.py\", line 2081, in run_line_magic\n",
      "    result = fn(*args,**kwargs)\n",
      "  File \"<decorator-gen-105>\", line 2, in matplotlib\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/IPython/core/magic.py\", line 188, in <lambda>\n",
      "    call = lambda f, *a, **k: f(*a, **k)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/IPython/core/magics/pylab.py\", line 100, in matplotlib\n",
      "    gui, backend = self.shell.enable_matplotlib(args.gui)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/IPython/core/interactiveshell.py\", line 2966, in enable_matplotlib\n",
      "    pt.activate_matplotlib(backend)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/IPython/core/pylabtools.py\", line 315, in activate_matplotlib\n",
      "    matplotlib.pyplot.switch_backend(backend)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/matplotlib/pyplot.py\", line 233, in switch_backend\n",
      "    _backend_mod, new_figure_manager, draw_if_interactive, _show = pylab_setup()\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/matplotlib/backends/__init__.py\", line 62, in pylab_setup\n",
      "    [backend_name], 0)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/ipykernel/pylab/backend_inline.py\", line 165, in <module>\n",
      "    _enable_matplotlib_integration()\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/ipykernel/pylab/backend_inline.py\", line 155, in _enable_matplotlib_integration\n",
      "    activate_matplotlib(backend)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/IPython/core/pylabtools.py\", line 315, in activate_matplotlib\n",
      "    matplotlib.pyplot.switch_backend(backend)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/matplotlib/pyplot.py\", line 231, in switch_backend\n",
      "    matplotlib.use(newbackend, warn=False, force=True)\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/matplotlib/__init__.py\", line 1410, in use\n",
      "    reload(sys.modules['matplotlib.backends'])\n",
      "  File \"/opt/conda/envs/python2/lib/python2.7/site-packages/matplotlib/backends/__init__.py\", line 16, in <module>\n",
      "    line for line in traceback.format_stack()\n",
      "\n",
      "\n",
      "  matplotlib.use('Agg')\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "You are not going crazy...\n"
     ]
    }
   ],
   "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 pickle\n",
    "\n",
    "import warnings\n",
    "print('Ignore deprecation warnings')\n",
    "warnings.filterwarnings(action='ignore', category=DeprecationWarning)\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 = os.getcwd()\n",
    "print('Data will be loaded from and saved to:'+homedir)\n",
    "os.chdir('/home/jovyan/work/notebooks/data/4cac25933f6448409cab97b293129b4f/4cac25933f6448409cab97b293129b4f/data/contents')\n",
    "print(homedir)\n",
    "\n",
    "# Import Landlab libraries\n",
    "import landslide_probability\n",
    "\n",
    "from landslide_probability_20191202 import LandslideProbability\n",
    "from landlab import RasterModelGrid\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\n",
    "from landlab.io.netcdf import read_netcdf\n",
    "from landlab.io.netcdf import write_netcdf\n",
    "\n",
    "import time\n",
    "st = time.time()"
   ]
  },
  {
   "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 Seattle City Light case study and is on HydroShare as [Slippery Future Data: Predicting future regional landslide probability using soil saturation](https://www.hydroshare.org/resource/01b486f301864828ba2cd9ab7ac77c4e/). 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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Data_ResourceID='01b486f301864828ba2cd9ab7ac77c4e'"
   ]
  },
  {
   "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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This resource already exists in your userspace.\n",
      "01b486f301864828ba2cd9ab7ac77c4e/\n",
      "|-- 01b486f301864828ba2cd9ab7ac77c4e/\n",
      "|   |-- bagit.txt\n",
      "|   |-- manifest-md5.txt\n",
      "|   |-- readme.txt\n",
      "|   |-- tagmanifest-md5.txt\n",
      "|   |-- data/\n",
      "|   |   |-- resourcemap.xml\n",
      "|   |   |-- resourcemetadata.xml\n",
      "|   |   |-- contents/\n",
      "|   |   |   |-- ASCII_Files/\n",
      "|   |   |   |   |-- glacier_scl_exclusion.txt\n",
      "|   |   |   |   |-- goodell_fire.txt\n",
      "|   |   |   |   |-- noca_landslidetype.txt\n",
      "|   |   |   |   |-- phi28-40_30m.txt\n",
      "|   |   |   |   |-- scl_coh_max.txt\n",
      "|   |   |   |   |-- scl_coh_min.txt\n",
      "|   |   |   |   |-- scl_coh_mod.txt\n",
      "|   |   |   |   |-- scl_dem_30m.txt\n",
      "|   |   |   |   |-- scldomain_header_30m.txt\n",
      "|   |   |   |   |-- scl_domain_mask_30m.txt\n",
      "|   |   |   |   |-- scl_lulc_8c.txt\n",
      "|   |   |   |   |-- scl_slope_tan.txt\n",
      "|   |   |   |   |-- scl_soildepth_1-3m.txt\n",
      "\n",
      "Do you want to overwrite these data [Y/n]? n\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<b>Found the following file(s) associated with this HydroShare resource.</b>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "ASCII_Files"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "These files are stored in a dictionary called <b>hs.content</b> for your convenience.  To access a file, simply issue the following command where MY_FILE is one of the files listed above: <pre>hs.content[\"MY_FILE\"] </pre> "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This is the location on the HydroShare JupyterHub server where the data has just been downloaded:\n",
      "/home/jovyan/work/notebooks/data/01b486f301864828ba2cd9ab7ac77c4e/01b486f301864828ba2cd9ab7ac77c4e/data/contents/ASCII_Files/\n"
     ]
    }
   ],
   "source": [
    "hs.getResourceFromHydroShare(Data_ResourceID)\n",
    "data_folder = '/home/jovyan/work/notebooks/data/'+ Data_ResourceID +'/'+Data_ResourceID+'/data/contents/ASCII_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": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['soil__density',\n",
       " 'soil__internal_friction_angle',\n",
       " 'soil__maximum_total_cohesion',\n",
       " 'soil__minimum_total_cohesion',\n",
       " 'soil__mode_total_cohesion',\n",
       " 'soil__saturated_hydraulic_conductivity',\n",
       " 'soil__thickness',\n",
       " 'soil__transmissivity',\n",
       " 'topographic__slope',\n",
       " 'topographic__specific_contributing_area']"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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 are is large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid, z) = read_esri_ascii(data_folder+'/scl_dem_30m.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 SCL extent, this takes ~60 sec to load each file on NCSA ROGER super computer. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load slope (this is in tan theta)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, slope) = read_esri_ascii(data_folder+'/scl_slope_tan.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 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+'/scl_coh_mod.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+'/scl_coh_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+'/scl_coh_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+'/phi28-40_30m.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+'/scl_soildepth_1-3m.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 to exclude areas outside the Skagit River Watershed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, mask) = read_esri_ascii(data_folder+'/scl_domain_mask_30m.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 within North Cascades National Park for later plotting. There are 6 classes of landslides: rock fall or topple (1), debris torrent (2), debris avalanche (3), slump or creep (4), snow avalanche-impacted (5), and sackung (6)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, slides) = read_esri_ascii(data_folder+'/noca_landslidetype.txt')\n",
    "grid.add_field('node', 'landslides', slides)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load landcover classes to use for examining patterns in results.  There are 8 classes: forest (41), shrubs (52), herbaceous (71), wetland (90), developed (21), barren (31), ice & snow (12), and water (11).  Exclude the last 2 classes from analysis as landslides based on this approach are not modeled in water bodies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, lulc) = read_esri_ascii(data_folder+'/scl_lulc_8c.txt')\n",
    "grid.add_field('node', 'landcover', lulc)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['landcover'], 12)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['landcover'], 11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load area of Goodell Creek Fire (2015) located within Skagit River Watershed. Use this area to adjust root cohesion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, fire) = read_esri_ascii(data_folder+'/goodell_fire.txt')\n",
    "grid.add_field('node', 'fire_mask', mask)\n",
    "#grid.set_nodata_nodes_to_closed(grid.at_node['fire_mask'], -9999)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load area of mapped glaciers and exclude from analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(grid1, glaciers) = read_esri_ascii(data_folder+'/glacier_scl_exclusion.txt')\n",
    "grid.add_field('node', 'glacier_mask', glaciers)\n",
    "grid.set_nodata_nodes_to_closed(grid.at_node['glacier_mask'], 1)"
   ]
  },
  {
   "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": 1,
   "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": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "key=np.arange(grid.number_of_nodes)\n",
    "yrs=44\n",
    "HSD_dict_mean=np.zeros(grid.number_of_nodes)\n",
    "HSD_dict_stddev=[]\n",
    "HSD_dict_annualmaxDWT=zip(np.random.rand(yrs),np.random.rand(yrs)*1.5,np.random.rand(yrs)*0.5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "def statHSD(HSDout,dsi,x1,x2,y1,y2):\n",
    "        #Create python dictionary of depth to water tabls\n",
    "        number_of_nodes=7412000\n",
    "        keys=np.arange(number_of_nodes)\n",
    "\n",
    "        HSD_dict_stat_annualmaxDWT={}\n",
    "        counter=0\n",
    "        for j in np.arange(x1,x2):\n",
    "            for k in np.arange(y1,y2):\n",
    "                one_location=dsi.isel(x=[j],y=[k]).to_array()\n",
    "                loc1list=np.array(one_location.variable)\n",
    "                b = list(itertools.chain(*loc1list))\n",
    "                c = list(itertools.chain(*b))\n",
    "                d = list(itertools.chain(*c))\n",
    "                mean_HSD_DWT=d.mean()\n",
    "                std_HSD_DWT=np.std(keys)\n",
    "                HSD_dict_stat_annualmaxDWT[counter] = {keys[counter]:zip(mean_HSD_DWT,std_HSD_DWT)}\n",
    "                \n",
    "                counter=counter+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.98713525  0.21162708  0.86321152 ...,  0.56937298  0.80983611\n",
      "  0.68153639]\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "'numpy.ndarray' object has no attribute 'to_list'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0mTraceback (most recent call last)",
      "\u001b[0;32m<ipython-input-61-ece55c8632c2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0mcounter\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcounter\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;32mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mHSD_dict_mean\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mHSD_dict_annualmaxDWT\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mHSD_dict_mean\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_list\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m: 'numpy.ndarray' object has no attribute 'to_list'"
     ]
    }
   ],
   "source": [
    "counter=0\n",
    "for i in range(grid.number_of_nodes):\n",
    "    HSD_dict_mean[counter]=np.random.rand(1)\n",
    "    \n",
    "    counter=counter+1\n",
    "print(HSD_dict_mean)\n",
    "HSD_dict_annualmaxDWT = {key:HSD_dict_mean.to_list()} "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combine dictionaries into __ordered__ parameters required for _data driven spatial_ distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0mTraceback (most recent call last)",
      "\u001b[0;32m<ipython-input-13-e4690993e40f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;31m#i=1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgrid\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnumber_of_nodes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m     \u001b[0mHSD_dict_annualmaxDWT\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mHSD_dict_annualmaxDWT\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;32mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mHSD_dict_annualmaxDWT\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: list index out of range"
     ]
    }
   ],
   "source": [
    "\n",
    "LS_prob1=LandslideProbability(grid, number_of_iterations=n,\n",
    "                     groundwater__depth_distribution='data_driven_spatial',\n",
    "                     groundwater__depth_HSD_inputs=[HSD_dict_annualmaxDWT])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input read for groundwater__recharge_distribution\n",
      "None\n",
      "Input read for  groundwater__depth_distribution\n",
      "lognormal_spatial\n"
     ]
    }
   ],
   "source": [
    "distribution1 = 'lognormal_spatial'\n",
    "grid_size = grid.number_of_nodes\n",
    "\n",
    "LS_prob1=LandslideProbability(grid, number_of_iterations=250,\n",
    "                     groundwater__depth_distribution='lognormal_spatial',\n",
    "                     groundwater__depth_mean=np.random.randint(2, 5, grid_size),\n",
    "                     groundwater__depth_standard_deviation=np.random.rand(grid_size))\n"
   ]
  },
  {
   "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 = 50"
   ]
  },
  {
   "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__depth_distribution=distribution,\n",
    "    groudwater__depth_HSD_inputs=[HSD_dict_annualmaxDWT])"
   ]
  },
  {
   "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. Stability analysis following fire"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a coy of original grid and then reduce root cohesion 30% where fire burned vegetation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "grid_fire=copy.deepcopy(grid)\n",
    "fire_reduction = 0.7\n",
    "grid_fire.at_node['soil__mode_total_cohesion']=grid.at_node['soil__mode_total_cohesion']*fire_reduction\n",
    "grid_fire.at_node['soil__minimum_total_cohesion']=grid.at_node['soil__minimum_total_cohesion']*fire_reduction\n",
    "grid_fire.at_node['soil__maximum_total_cohesion']=grid.at_node['soil__maximum_total_cohesion']*fire_reduction\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the cohesion values of the 2 grids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.max(grid.at_node['soil__minimum_total_cohesion'][grid.core_nodes]))\n",
    "print(np.min(grid.at_node['soil__minimum_total_cohesion'][grid.core_nodes]))\n",
    "print(np.max(grid_fire.at_node['soil__minimum_total_cohesion'][grid.core_nodes]))\n",
    "print(np.min(grid_fire.at_node['soil__minimum_total_cohesion'][grid.core_nodes]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll run the landslide component with the adjusted cohesion, everything else kept constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LS_probFire = LandslideProbability(grid_fire,number_of_iterations=iterations,\n",
    "    groudwater__depth_distribution=distribution,\n",
    "    groudwater__depth_HSD_inputs=[HSD_dict_annualmaxDWT])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LS_probFire.calculate_landslide_probability()\n",
    "print('Landslide probability successfully calculated')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9. 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('SCL_elevation.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Excluded areas from the analysis are shown in black, including outside the Skagit River Watershed and inside the watershed that are water, snow & ice, glaciers, and wetlands, and slopes=0 degrees. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot slope overlaid with mapped 4 primary landslide types within national park. 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": [
    "Let's compare the Probability of Failure between the before and after fire"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure('Probability of Failures')\n",
    "xticks = np.arange(-0.1, 1.6, 0.8)\n",
    "\n",
    "\n",
    "ax1 = fig.add_subplot(221)\n",
    "ax1.xaxis.set_visible(False)\n",
    "imshow_grid_at_node(grid, 'landslide__probability_of_failure', \n",
    "                    plot_name='Probability of Failure - Pre-fire',\n",
    "                    allow_colorbar=True, cmap='OrRd',\n",
    "                    grid_units=('coordinates', 'coordinates'), shrink=0.75,\n",
    "                    var_name='Probability of Failure', var_units='no units')\n",
    "\n",
    "ax2 = fig.add_subplot(212)\n",
    "ax2.set_xticks(xticks)\n",
    "imshow_grid_at_node(grid_fire, 'landslide__probability_of_failure', \n",
    "                    plot_name='Probability of Failure - Post-fire',\n",
    "                    allow_colorbar=True, cmap='OrRd',\n",
    "                    grid_units=('coordinates', 'coordinates'), shrink=0.75,\n",
    "                    var_name='Probability of Failure', var_units='no units')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How long did the code above take to run?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print 'Elapsed time is %3.2f seconds' % (time.time() - st)"
   ]
  },
  {
   "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')"
   ]
  }
 ],
 "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}