¶

- Load data from a regional Landlab landslide model (Strauch et al., 2018) developed for the North Cascades National Park, WA USA, published on HydroShare.
- Define a geographic subset (Thunder Creek watershed) within the study region.
- Explore landslide probability sensitivity to fire by adjusting the cohesion parameter for Thunder Creek.
- Save results to a new HydroShare resource.

The shallow landslide model you will is based on a spatially distributed Monte Carlo solution of the infinite slope stability model. Detailes of the model and the study site are described in Strauch et al. (2018). Please see the end of this Notebook for Acknowledgements and Citation information.

Click in each shaded cell below and "shift + enter" to run each code block. Alternatively, you can run groups of cells by clicking "Cell" on the menu above and selecting your run options from the pull-down menu. This is also where you can clear outputs from previous runs.

If an error occurs, click on *Kernal* and *Restart and Clear Outputs* in the menu above.

** 1.1 Infinite Slope Factor of Safety Equation **

This equation predicts the ratio of stabilizing to destabilizing forces on a hillslope plane. The properties are assumed to represent an infinte plane, neglecting the boundary conditions around the landslide location. When FS<1, the slope is instable.

C* can be calculated by the ratio of the sum of root cohesion, Cr, and soil cohesion, Cs, to the product of soil depth, density, and gravity.

Relative wetness is the ratio of depth of water subsurface flow above an impervious layer to the depth of soil.

R: recharge rate to water table (m/d) T: soil transmissivity (m^2/d) a: specific catchment area (m).

** 1.2 Monte Carlo solution of the FS equation **

Below the Monte Carlo approach used by Strauch et al. (2018) is illustrated in a figure. At each node of the model grid variables for soil and vegetation are generated from triangular distributions. Recharge to water table is obtained from hydrologic model simulations by selecting the largest daily value for each year. Probability of shallow landslide initiation is defined as the ratio of number of times FS<1 to the total sample size of FS calculations.

This Jupyter Notebook runs the Landlab LandslideProbability component on a 30-m digital elevation model (DEM) for Thunder Creek, a subwaterhsed in the larger study domain developed by Strauch et al. (2018), using *data driven spatial* recharge distribution as described in the paper (https://www.earth-surf-dynam.net/6/1/2018/).

To run a landslide demonstration using the paper data and approach we will:

1) Import data from North Cascades National Park (NOCA: study area)

2) Review data needed as input for the landslide model

3) Create a RasterModelGrid based on a 30-m DEM - subset to Thunder Creek watershed

4) Assign data fields used to calculate landslide probability

5) Specify recharge option as *data driven spatial* and access Python dictionaries to generate recharge distributions

6) Set Number of iterations to run Monte Carlo simulation

7) Run Landlab LandslideProbability component

8) Run the model again to simulate post-fire conditions,

9) Display and visualize results of stability analysis

10) Save Notebook and Results back to HydroShare

The estimated time to run this Notebook is 30-60 minutes with 20 minutes of computation time.

Update or add libraries that are not installed on the CUAHSI JupyterHub server. At the time of publication, the Landlab library was dependent on a different numpy library than installed in the software environment. The code block below was added as an example of how to update a library to a new version and restart the notebook kernel. If this update is not needed, comment with a <#>. If other updates are needed, add them below.

In [1]:

```
# comment this out if this update is not needed
#!conda update -y numpy
#from IPython.display import display_html
#display_html("<script>Jupyter.notebook.kernel.restart()</script>",raw=True)
```

To run this notebook, we must import several libraries. The hs_utils library provides functions for interacting with HydroShare, including resource querying, dowloading. and creation. Additional libraries support the functions of Landlab. The CUAHSI JupyterHub server provides many Python packages and libraries, but to add additional libraries to your personal user space, use the cell below. To request an Installation to the server, visit https://github.com/hydroshare/hydroshare-jupyterhub, create a New Issue, and add the label 'Installation Request'.

In [2]:

```
#Import standard Python utilities for calculating and plotting
import six
import os
import matplotlib as mpl
mpl.use('agg')
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import pandas as pd
import pickle as pickle
from datetime import datetime, timedelta
import geopandas as gpd
# Import Landlab libraries
from landlab.components import LandslideProbability
from landlab import imshow_grid_at_node
from landlab.io import read_esri_ascii
from landlab.io import write_esri_ascii
from landlab.plot import imshow_grid
from landlab import CORE_NODE, CLOSED_BOUNDARY
#Import utilities for importing and exporting to HydroShare
from hstools import hydroshare
# Import general tools
import time
from collections import defaultdict
st = time.time()
%matplotlib inline
```

After importing libraires, we now establish a secure connection with HydroShare by instantiating the hydroshare class that is defined within hs_utils. In addition to connecting with HydroShare, this command also sets and prints environment variables for several parameters that will be useful for saving work back to HydroShare.

In [3]:

```
hs=hydroshare.hydroshare()
homedir=os.getcwd()
#homedir = os.path.join('/home/jovyan/work/notebooks/data', str(homeresid),str(homeresid),'data/contents/')
#print('This is a basic Unix folder structure.')
#print('Data will be loaded from and saved to:'+homedir)
```

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.

Strauch et al (2018) pre-processed the data for the North Cascades National Park Complex case study and is on HydroShare as Regional landslide hazard using Landlab - NOCA Data. Here the first task is find this resource on Hydrohsare. We will click on the link to see the published data repository on HydroShare and collect the resource ID of the data. The resource ID can be found in the "How to cite" box, and it will be the series of numbers and letters following "hs.". Here's the copied resource: http://dx.doi.org/10.4211/hs.a5b52c0e1493401a815f4e77b09d352b citation. Now we copy this ID and introduce it as "Data_ResourceID=" in the code below.

To learn more about this data visit Regional landslide hazard using Landlab - NOCA Data

In [4]:

```
Data_ResourceID='a5b52c0e1493401a815f4e77b09d352b'
```

We will execute the next cell to download data from HydroShare iRods database to your personal user space - this may take a few minutes.

In [5]:

```
hs.getResourceFromHydroShare(Data_ResourceID)
data_folder = os.path.join(os.getcwd(), hs.getContentPath(Data_ResourceID), 'Data files')
print('This HydroShare resource has been downloaded to: %s' % data_folder)
```

This section loads metadata associated with the Landlab component.

To view the code source for this component, visit Landlab on Github or Download the landslide_probability.py python file

Check the list of data inputs that the component needs.

In [6]:

```
sorted(LandslideProbability.input_var_names)
```

Out[6]:

Review the details of what each variable represents.

In [7]:

```
LandslideProbability._var_doc
```

Out[7]:

Check the units of each variable.

In [8]:

```
LandslideProbability._var_units
```

Out[8]:

Now we will establish a RasterModelGrid based on a DEM for assigning our variables to. Nodes are the center point of grid cells or pixels that are 30 m by 30 m in this example.

The shapefile that you will download is a GIS point file with a table containing a grid_code for each 30m DEM cell in the Thunder Creek watershed. We import this table to generate a list of locations we want to keep "Open" as active nodes in this Landlab application for Thunder Creek.

Here we will establish the watershed domain for our modeling study. In step 2.2.2, we have downloaded the DEM of the entire NOCA study region of Strauch et al. (2018) as well as other input data. This includes NoData resource that was created using a mask of Thunder Creek will be downladed and used to set no data nodes as inactive nodes (e.g., -9999). This step will establish boundary conditions. In your final run when you output modeled probability of landslide initiation, you will notice gaps in the model results at some ridge tops and peaks. This results from excluding glaciated areas and bedrock.

This might take a few minutes as the park is large (2,757 km2).

This shapefile resource was uploaded to HydroShare to generate an interoperable Notebook. To learn more about this data visit Thunder Creek Landlab Node ID Point Shapefile

The shapefile table contains Landlab Gridcode value and Albers Conical X and Y values, spatial projection is WGS84.

In [9]:

```
Data_ResourceID='8bf8de77227c4dcba4816dbe15e55687'
hs.getResourceFromHydroShare(Data_ResourceID)
```

In [10]:

```
Node_path = os.path.join(hs.getContentPath(Data_ResourceID), 'Thunder_node_id_WGS.shp')
NodeID_shpfile=gpd.GeoDataFrame.from_file(Node_path)
```

Print a sample portion (bottom) of the table

- Col 1 = FID or POINTID
- Col 2 = landlab grid code for NOCA region
- Col 3 Point_X - the Alber Conical Latitude (original dataset)
- Col 4 Point_Y - the Alber Conical Longitude (original dataset)
- Geometry of lat/lon in WGS coordinate system of this GIS point file for launching in HydroShareGIS App

In [11]:

```
NodeID_shpfile.tail()
```

Out[11]:

Make an array using Grid_Code column to develop a node based mask. No output from this command.

In [12]:

```
filtercriteria = np.array(NodeID_shpfile.GRID_CODE)
```

Create Landlab RasterModelGrid using DEM grid with elevation - this takes approximately 60 sec for the North Cascades National Park (NOCA).

In [13]:

```
(grid, z) = read_esri_ascii(data_folder+'/elevation.txt',name='topographic__elevation')
grid.at_node.keys()
```

Out[13]:

The original NOCA dem dataset has nodes closed outside of NOCA. All nodes insides of NOCA are open. To close all all nodes except those inside Thunder Creek:

- Close all nodes in the landlab RasterModelGrid.
- Open the Landlab RasterModelGrid only within the Thunder Creek watershed.
- Assign Core Nodes (these are the only nodes where computation will be executed).

Plot the elevation grid of NOCA

In [14]:

```
imshow_grid(grid, 'topographic__elevation', limits=(0, 3000))
```

Plot where the watershed is witin NOCA and zoom in.

In [15]:

```
grid.status_at_node[grid.nodes.flatten()] = CLOSED_BOUNDARY
grid.status_at_node[filtercriteria] = CORE_NODE
```

In [16]:

```
# In case the DEM has no data values inside the subset area, set boundary conditions closed.
grid.set_nodata_nodes_to_closed(grid.at_node['topographic__elevation'], -9999)
```

In [17]:

```
fig = plt.figure('Limit NOCA computational nodes to Thunder Creek extent')
xticks = np.arange(-0.1, 0.8, 0.4)
ax1 = fig.add_subplot(221)
ax1.xaxis.set_visible(False)
imshow_grid(grid, 'topographic__elevation', limits=(0, 3000),plot_name='NOCA extent',
allow_colorbar=False,color_for_closed='white')
ax2 = fig.add_subplot(222)
ax2.xaxis.set_visible(False)
imshow_grid(grid, 'topographic__elevation', limits=(0, 3000),plot_name='Thunder Creek extent',color_for_closed='white')
plt.xlim(25000, 55000)
plt.ylim(30000, 55000)
```

Out[17]:

Confirm the size of the grid, nodes located every 30 m.

In [18]:

```
grid.number_of_nodes
```

Out[18]:

Confirm the size of the core nodes where we'll run our model, a subset of the nodes for the watershed.

In [19]:

```
grid.number_of_core_nodes
```

Out[19]:

This will be used to calculate shallow landslide probability and set boundary conditions. THe data we attach in this step was downloaded in Step 2.2.2.

- Load data from ascii text file
- Add this data as node variable to the Thunder Creek grid
- Set boundary conditions

For the entire NOCA extent, this takes ~60 sec to load each file using the CUAHSI JupyterHub server.

In [20]:

```
(grid1, slope) = read_esri_ascii(data_folder+'/slope_tang17d.txt')
grid.add_field('node', 'topographic__slope', slope)
grid.set_nodata_nodes_to_closed(grid.at_node['topographic__slope'], -9999)
grid.set_nodata_nodes_to_closed(grid.at_node['topographic__slope'], 0.0)
```

In [21]:

```
print(np.max(grid.at_node['topographic__slope'][grid.core_nodes]))
print(np.min(grid.at_node['topographic__slope'][grid.core_nodes]))
```

In [22]:

```
(grid1, ca) = read_esri_ascii(data_folder+'/cont_area.txt')
grid.add_field('node', 'topographic__specific_contributing_area', ca)
grid.set_nodata_nodes_to_closed(grid.at_node['topographic__specific_contributing_area'], -9999)
```

In [23]:

```
(grid1, T) = read_esri_ascii(data_folder+'/transmis.txt')
grid.add_field('node', 'soil__transmissivity', T)
grid.set_nodata_nodes_to_closed(grid.at_node['soil__transmissivity'], -9999)
```

This takes ~3 minutes because 3 cohesion fields are provided to create more flexibility in how cohesion is distributed on the landscape with different vegetation.

In [24]:

```
(grid1, C) = read_esri_ascii(data_folder+'/cohesion_mode.txt')
C[C == 0.0] = 1.0 # ensure minimum is >0 Pa for use in distributions generation
grid.add_field('node', 'soil__mode_total_cohesion', C)
grid.set_nodata_nodes_to_closed(grid.at_node['soil__mode_total_cohesion'], -9999)
(grid1, C_min) = read_esri_ascii(data_folder+'/cohesion_min.txt')
grid.add_field('node', 'soil__minimum_total_cohesion', C_min)
grid.set_nodata_nodes_to_closed(grid.at_node['soil__minimum_total_cohesion'], -9999)
(grid1, C_max) = read_esri_ascii(data_folder+'/cohesion_max.txt')
grid.add_field('node', 'soil__maximum_total_cohesion', C_max)
grid.set_nodata_nodes_to_closed(grid.at_node['soil__maximum_total_cohesion'], -9999)
```

In [25]:

```
(grid1, phi) = read_esri_ascii(data_folder+'/frict_angle.txt')
grid.add_field('node', 'soil__internal_friction_angle', phi)
grid.set_nodata_nodes_to_closed(grid.at_node['soil__internal_friction_angle'], -9999)
```

In this example, we assign all nodes a constant value.

In [26]:

```
grid['node']['soil__density'] = 2000*np.ones(grid.number_of_nodes)
```

In [27]:

```
(grid1, hs2) = read_esri_ascii(data_folder+'/soil_depth.txt')
grid.add_field('node', 'soil__thickness', hs2)
grid.set_nodata_nodes_to_closed(grid.at_node['soil__thickness'], -9999)
```

In [28]:

```
(grid1, slides) = read_esri_ascii(data_folder+'/landslide_type.txt')
grid.add_field('node', 'landslides', slides)
```

Out[28]:

Recharge in this model represents the annual maximum recharge in mm/day generated within the upslope contributing area of each model element. This corresponds to the wettest conditions expected annually, which would provide the highest pore-water pressure in a year. Details of this approach can be found in Strauch et al. (2018).

In [29]:

```
distribution = 'data_driven_spatial'
```

These contain HSD_id and fractional drainage at each node and recharge dictionaries. HSD is the Hydrologic Source Domain, which is the VIC hydrologic model data in this case study at ~5x6 km2 grid size. The 'pickle' utility loads existing dictionaries.

In [30]:

```
# dict of node id (key) and HSD_ids (values)
HSD_id_dict = pickle.load(open(data_folder+'/dict_uniq_ids.p', 'rb'),encoding='latin1')
# dict of node id (key) and fractions (values)
fract_dict = pickle.load(open(data_folder+'/dict_coeff.p', 'rb'),encoding='latin1')
# dict of HSD id (key) with arrays of recharge (values)
HSD_dict = pickle.load(open(data_folder+'/HSD_dict.p', 'rb'),encoding='latin1')
```

This sequence of parameters is required for *data driven spatial* distribution in the component. Recharge is this model is unique to each node represented by an array.

In [31]:

```
HSD_inputs = [HSD_dict, HSD_id_dict, fract_dict]
```

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 Monte Carlo samples, but the default is 250. The larger the Monte Carlos sample size, the longer the program runs, but the more precise the probability of failure results become. Strauch et al. (2018) sampled 3,000 times for each parameter in each model grid.

In [32]:

```
iterations = 250
```

To run the landslide model, we first instantiate the LandslideProbability component with the above parameters, as well as the grid and number of samples we specified before. Instantiation creates an instance of a class called LS_prob.

No outputs are generated by this command as it is setting up the recharge and instantiating the component.

In [33]:

```
LS_prob = LandslideProbability(grid,
number_of_iterations=iterations,
groudwater__recharge_distribution=distribution,
groudwater__recharge_HSD_inputs=HSD_inputs)
```

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 cell below 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.

In [34]:

```
LS_prob.calculate_landslide_probability()
print('Landslide probability successfully calculated')
```

What is the maximum probabilty of failure we found?

In [35]:

```
np.max(grid.at_node['landslide__probability_of_failure'][grid.core_nodes])
```

Out[35]:

The outputs of landslide model simulation are:

In [36]:

```
sorted(LS_prob.output_var_names)
```

Out[36]:

This simulation generates a probability value for each core node.

In [37]:

```
grid.at_node['landslide__probability_of_failure']
```

Out[37]:

This simulation generates a probability of saturation value for each core node as well.

In [38]:

```
grid.at_node['soil__probability_of_saturation']
```

Out[38]:

-9999 means there is no data for that cell, it is a closed node.

Make a 'grid_fire' version of 'grid' for which we will give post-fire cohesion parameter values 30% of the original cohesion. This is a crude estimation of the lowest root cohesion after tree removal based on a combined decay and regrowth model (e.g., Sidle, 1992).

Sidle, R. C. (1992), A theoretical model of the effects of timber harvesting on slope stability, Water Resour. Res., 28(7), 1897–1910.

In [39]:

```
import copy
grid_fire=copy.deepcopy(grid)
grid_fire.at_node['soil__mode_total_cohesion']=grid.at_node['soil__mode_total_cohesion']*0.3
grid_fire.at_node['soil__minimum_total_cohesion']=grid.at_node['soil__minimum_total_cohesion']*0.3
grid_fire.at_node['soil__maximum_total_cohesion']=grid.at_node['soil__maximum_total_cohesion']*0.3
```

In [40]:

```
print("this is the highest mode value for coehsion before fire across the domain")
print(np.max(grid.at_node['soil__mode_total_cohesion'][grid.core_nodes]))
print("this is the lowest mode value for coehsion before fire across the domain")
print(np.min(grid.at_node['soil__mode_total_cohesion'][grid.core_nodes]))
print("this is the highest mode value for coehsion after fire across the domain")
print(np.max(grid_fire.at_node['soil__mode_total_cohesion'][grid.core_nodes]))
print("this is the lowest mode value for coehsion after fire across the domain")
print(np.min(grid_fire.at_node['soil__mode_total_cohesion'][grid.core_nodes]))
```

Now we'll run the landslide component with the adjusted cohesion, everything else kept constant.

In [41]:

```
LS_probFire = LandslideProbability(grid_fire,number_of_iterations=iterations,
groudwater__recharge_distribution=distribution,
groudwater__recharge_HSD_inputs=HSD_inputs)
print('Post-fire cohesion successfully instantiated')
```

In [42]:

```
LS_probFire.calculate_landslide_probability()
print('Landslide probability successfully calculated')
```

In [43]:

```
np.max(grid_fire.at_node['landslide__probability_of_failure'][grid.core_nodes])
```

Out[43]:

The outputs of landslide model simulation are:

In [44]:

```
sorted(LS_probFire.output_var_names)
```

Out[44]:

This simulation generates a probability value for each core node.

In [45]:

```
grid_fire.at_node['landslide__probability_of_failure']
```

Out[45]:

Set plotting parameters

In [46]:

```
mpl.rcParams['xtick.labelsize'] = 15
mpl.rcParams['ytick.labelsize'] = 15
mpl.rcParams['lines.linewidth'] = 1
mpl.rcParams['axes.labelsize'] = 18
mpl.rcParams['legend.fontsize'] = 15
```

Plot elevation

In [47]:

```
plt.figure('Elevations from the DEM [m]')
imshow_grid_at_node(grid, 'topographic__elevation', cmap='terrain',
grid_units=('coordinates', 'coordinates'),
shrink=0.75, var_name='Elevation', var_units='m',color_for_closed='white')
plt.xlim(25000, 55000)
plt.ylim(30000, 55000)
plt.savefig('NOCA_elevation.png')
```

Excluded areas from the analysis are shown in white, including outside the park and inside the park areas that are water bodies, snow, glaciers, wetlands, exposed bedrock, and slopes <= 17 degrees.

Plot slope overlaid with mapped landslide types. Takes about a few minutes.

In [48]:

```
plt.figure('Landslides')
ls_mask1 = grid.at_node['landslides'] != 1.0
ls_mask2 = grid.at_node['landslides'] != 2.0
ls_mask3 = grid.at_node['landslides'] != 3.0
ls_mask4 = grid.at_node['landslides'] != 4.0
overlay_landslide1 = np.ma.array(grid.at_node['landslides'], mask=ls_mask1)
overlay_landslide2 = np.ma.array(grid.at_node['landslides'], mask=ls_mask2)
overlay_landslide3 = np.ma.array(grid.at_node['landslides'], mask=ls_mask3)
overlay_landslide4 = np.ma.array(grid.at_node['landslides'], mask=ls_mask4)
imshow_grid_at_node(grid, 'topographic__slope', cmap='pink',
grid_units=('coordinates', 'coordinates'), vmax=2.,
shrink=0.75, var_name='Slope', var_units='m/m',color_for_closed='white',)
imshow_grid_at_node(grid, overlay_landslide1, color_for_closed='None',
allow_colorbar=False, cmap='cool')
imshow_grid_at_node(grid, overlay_landslide2, color_for_closed='None',
allow_colorbar=False, cmap='autumn')
imshow_grid_at_node(grid, overlay_landslide3, color_for_closed='None',
allow_colorbar=False, cmap='winter')
imshow_grid_at_node(grid, overlay_landslide4, color_for_closed='None',
allow_colorbar=False,cmap='summer')
#plt.savefig('NOCA_Landslides_on_Slope.png')
plt.xlim(25000, 55000)
plt.ylim(30000, 55000)
```

Out[48]:

Legend to mapped landslides: blue - debris avalanches, cyan - falls/topples, red - debris torrents, and green - slumps/creeps

Plot of soil depth (m)

In [49]:

```
plt.figure('Soil Thickness')
imshow_grid_at_node(grid, 'soil__thickness', cmap='copper_r',
grid_units=('coordinates', 'coordinates'), shrink=0.75,
var_name='Soil Thickness', var_units='m', color_for_closed='white')
plt.xlim(25000, 55000)
plt.ylim(30000, 55000)
#plt.savefig('NOCA_SoilDepth.png')
```

Out[49]:

Plot probability of saturation

In [50]:

```
fig = plt.figure('Probability of Saturation')
xticks = np.arange(-0.1, 0.8, 0.4)
ax1 = fig.add_subplot(221)
ax1.xaxis.set_visible(True)
imshow_grid(grid, 'soil__probability_of_saturation',cmap='YlGnBu',
limits=((0), (1)),plot_name='Pre-Fire',
allow_colorbar=False,grid_units=('coordinates', 'coordinates'), color_for_closed='white')
plt.xlim(25000, 55000)
plt.ylim(30000, 55000)
ax2 = fig.add_subplot(222)
ax2.xaxis.set_visible(True)
ax2.yaxis.set_visible(False)
imshow_grid(grid_fire, 'soil__probability_of_saturation',cmap='YlGnBu',
limits=((0), (1)),plot_name='Post-Fire',grid_units=('coordinates', 'coordinates'), color_for_closed='white')
plt.xlim(25000, 55000)
plt.ylim(30000, 55000)
```

Out[50]:

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.

Plot probability of failure; Compare this with the elevation and slope maps.

In [51]:

```
fig = plt.figure('Probability of Failure')
xticks = np.arange(-0.1, 0.8, 0.4)
ax1 = fig.add_subplot(221)
#ax1.xaxis.set_visible(False)
imshow_grid(grid, 'landslide__probability_of_failure',cmap='OrRd',
limits=((0), (1)),plot_name='Pre-fire',
allow_colorbar=False,grid_units=('coordinates', 'coordinates'), color_for_closed='white')
plt.xlim(25000, 55000)
plt.ylim(30000, 55000)
ax2 = fig.add_subplot(222)
ax2.yaxis.set_visible(False)
imshow_grid(grid_fire, 'landslide__probability_of_failure',cmap='OrRd',
limits=((0), (1)),plot_name='Post-fire',grid_units=('coordinates', 'coordinates'), color_for_closed='white')
plt.xlim(25000, 55000)
plt.ylim(30000, 55000)
#plt.savefig('Probability_of_Failure_Original_Fire.png')
```

Out[51]:

The map of probability of failure shows higher probabilities in a post-fire scenario with reduced cohesion.

We can compare the range of probabilities before and after the simulated fire using a cumulative distribution plot.

In [52]:

```
X1 = np.sort(np.array(grid.at_node['landslide__probability_of_failure'][grid.core_nodes]))
Y1 = np.arange(1, len(X1)+1)/len(X1)
X2 = np.sort(np.array(grid_fire.at_node['landslide__probability_of_failure'][grid.core_nodes]))
Y2 = np.arange(1, len(X2)+1)/len(X2)
fig = plt.figure('Empirical CDF of Probability of Failure')
plt.plot(X1,Y1, marker=' ', color='black', linewidth=2, label='Pre-fire')
plt.plot(X2,Y2, marker=' ', color='red', linewidth=2, linestyle='dashed', label='Post-fire')
plt.xlabel('Probability of Failure')
plt.ylabel('Fraction of Area')
plt.legend()
```

Out[52]:

We can interpret this figure as a higher fraction of the landscape has a higher probability of failure (red) after fire than before fire (black).

To review the fields assigned to the grid, simply execute the following command.

In [53]:

```
grid.at_node
```

Out[53]:

Export data from model run: FS probability, mean Reletive wetness, probability of saturation

In [54]:

```
import pandas as pd
core_nodes = grid.core_nodes
data_extracted = {'Prob_fail_std': np.array(
grid.at_node['landslide__probability_of_failure'][grid.core_nodes]),
'mean_RW_std': np.array(grid.at_node['soil__mean_relative_wetness']
[grid.core_nodes]),'prob_sat_std': np.array(
grid.at_node['soil__probability_of_saturation'][grid.core_nodes]),
'Prob_fail_fire': np.array(
grid_fire.at_node['landslide__probability_of_failure'][grid_fire.core_nodes]),
'mean_RW_fire': np.array(grid_fire.at_node['soil__mean_relative_wetness']
[grid_fire.core_nodes]),'prob_sat_fire': np.array(
grid_fire.at_node['soil__probability_of_saturation'][grid_fire.core_nodes])}
headers = ['Prob_fail_std','mean_RW_std','prob_sat_std','Prob_fail_fire','mean_RW_fire','prob_sat_fire']
df = pd.DataFrame(data_extracted, index=core_nodes, columns=(headers))
df.to_csv('Landslide_std_fire.csv')
```

Make ascii files for raster creation in GIS

In [55]:

```
write_esri_ascii('prob_fail_std.txt',grid,names='landslide__probability_of_failure')
write_esri_ascii('mean_RW_std.txt',grid,names='soil__mean_relative_wetness')
write_esri_ascii('prob_sat_std.txt',grid,names='soil__probability_of_saturation')
write_esri_ascii('prob_fail_fire.txt',grid_fire,names='landslide__probability_of_failure')
write_esri_ascii('mean_RW_fire.txt',grid_fire,names='soil__mean_relative_wetness')
write_esri_ascii('prob_sat_fire.txt',grid_fire,names='soil__probability_of_saturation')
```

Out[55]:

How long did the code above take to run?

In [56]:

```
print('Elapsed time is %3.2f seconds' % (time.time() - st))
```

Using the `hs_utils`

library, the results of the Geoprocessing steps above can be saved back into HydroShare. First, define all of the required metadata for resource creation, i.e. *title*, *abstract*, *keywords*, *content files*. In addition, we must define the type of resource that will be created, in this case *genericresource*.

** Note:** Make sure you save the notebook at this point, so that all notebook changes will be saved into the new HydroShare resource.

** Option A** : define the resource from which this "NEW" content has been derived. This is one method for tracking resource provenance.

Create list of files to save to HydroShare. Verify location and names.

In [57]:

```
ThisNotebook='replicate_landslide_model_for_fire.ipynb' #check name for consistency
files=[os.path.join(homedir, ThisNotebook),
os.path.join(homedir, 'Landslide_std_fire.csv'),
os.path.join(homedir, 'prob_fail_std.txt'),
os.path.join(homedir, 'mean_RW_std.txt'),
os.path.join(homedir, 'prob_sat_std.txt'),
os.path.join(homedir, 'prob_fail_fire.txt'),
os.path.join(homedir, 'mean_RW_fire.txt'),
os.path.join(homedir, 'prob_sat_fire.txt')]
for f in files:
if not os.path.exists(f):
print('\n[ERROR] invalid path: %s' % f)
```

In [58]:

```
help(hs)
```

In [59]:

```
# Save the results back to HydroShare
# title for the new resource
title = 'Landslide Model run with Monte Carlo from NOCA Observatory - Thunder Creek with Fire'
# abstract for the new resource
abstract = 'This a reproducible demonstration of the landslide modeling results from eSurf paper: Strauch et al. (2018) '
# keywords for the new resource
keywords = ['landslide', 'climate', 'VIC','saturation','relative wetness','fire','Geohackweek']
# Hydroshare resource type
rtype = 'compositeresource'
# create the new resource
resource_id = hs.createHydroShareResource(abstract,
title,
keywords=keywords,
resource_type=rtype,
content_files=files,
public=False)
```

This Notebooks serves as content for:

Bandaragoda, C. J., A. Castronova, E. Istanbulluoglu, R. Strauch, S. S. Nudurupati, J. Phuong, J. M. Adams, et al. “Enabling Collaborative Numerical Modeling in Earth Sciences Using Knowledge Infrastructure.” Environmental Modelling & Software, April 24, 2019. https://doi.org/10.1016/j.envsoft.2019.03.020.

and

Bandaragoda, C., A. M. Castronova, J. Phuong, E. Istanbulluoglu, S. S. Nudurupati, R. Strauch, N. Lyons, K. Barnhart (2019). Enabling Collaborative Numerical Modeling in Earth Sciences using Knowledge Infrastructure: Landlab Notebooks, HydroShare, http://www.hydroshare.org/resource/fdc3a06e6ad842abacfa5b896df73a76

This notebook was developed from code written by Ronda Strauch as part of her Ph.D. disseration at the University of Washington.

Use or citation of this notebook should also reference:

Strauch R., Istanbulluoglu E., Nudurupati S.S., Bandaragoda C., Gasparini N.M., and G.E. Tucker (2018). A hydro-climatological approach to predicting regional landslide probability using Landlab. Earth Surf. Dynam., 6, 1–26. https://www.earth-surf-dynam.net/6/1/2018/

In [ ]:

```
```