landslide_probability_20191202.py

#!/usr/env/python
print('You are not going crazy...')

"""Landlab component that simulates landslide probability of failure as well
as mean relative wetness and probability of saturation.

Relative wetness and factor-of-safety are based on the infinite slope
stability model driven by topographic and soils inputs and recharge provided
by user as inputs to the component. For each node, component simulates mean
relative wetness as well as the probability of saturation based on Monte Carlo
simulation of relative wetness where the probability is the number of
iterations with relative wetness >= 1.0 divided by the number of iterations.
Probability of failure for each node is also simulated in the Monte Carlo
simulation as the number of iterations with factor-of-safety <= 1.0
divided by the number of iterations.

.. codeauthor:: R.Strauch, E.Istanbulluoglu, & S.S.Nudurupati
University of Washington

Ref 1: Strauch et al., (2018) "A hydroclimatological approach to predicting regional landslide probability using Landlab", Earth Surf. Dynam., 6, 1-26 .

Ref 2: 'The Landlab LandslideProbability Component User Manual'

Created on Thu Aug 20, 2015
Last edit June 7, 2017

.. codeauthor:: R.Strauch, C.Bandaragoda
Seattle City Light, University of Washington

Ref 3: Strauch et al., (2020) "Slippery future 

Ref 4: Strauch, R., C. Raymond, C. Bandaragoda, N. Cristea (2019). Landslide Hazard Modeling in the Skagit Basin, HydroShare, http://www.hydroshare.org/resource/70d746c7da584ae6bd2f88deb5a4c188

Landslide Component Changes: Add depth to water table to optional input, alternative to recharge.
Created on May 7, 2019
Last edit Oct 22, 2019


"""

# %% Import Libraries
from landlab import Component
from landlab.utils.decorators import use_file_name_or_kwds
import numpy as np
import scipy.constants
from scipy import interpolate
from statsmodels.distributions.empirical_distribution import ECDF
import copy
import pandas as pd
# %% Instantiate Object


class LandslideProbability(Component):
    """
    Landlab component designed to calculate probability of failure at
    each grid node based on the infinite slope stability model
    stability index (Factor of Safety).

    The driving force for failure is provided by the user in the form of
    groundwater recharge; four options for providing recharge are supported.
    The model uses topographic and soil characteristics provided as input
    by the user.

    The main method of the LandslideProbability class is
    calculate_landslide_probability(), which calculates the mean soil relative
    wetness, probability of soil saturation, and probability of failure at
    each node based on a Monte Carlo simulation.

    Construction::
    Option 1 - Uniform recharge
        LandslideProbability(grid, number_of_iterations=250,
        groundwater__recharge_distribution='uniform', 
        groundwater__recharge_min_value=5.,
        groundwater__recharge_max_value=121.)
    Option 2 - Lognormal recharge
        LandslideProbability(grid, number_of_iterations=250,
        groundwater__recharge_distribution='lognormal', 
        groundwater__recharge_mean=30.,
        groundwater__recharge_standard_deviation=0.25)
    Option 3 - Lognormal_spatial recharge
        LandslideProbability(grid, number_of_iterations=250,
        groundwater__recharge_distribution='lognormal_spatial', 
        groundwater__recharge_mean=np.random.randint(20, 120, grid_size),
        groundwater__recharge_standard_deviation=np.random.rand(grid_size))
    Option 4 - Data_driven_spatial recharge
        LandslideProbability(grid, number_of_iterations=250,
        groundwater__recharge_distribution='data_driven_spatial', 
        groundwater__recharge_HSD_inputs=[HSD_dict, HSD_id_dict,
        fract_dict])
    Option 5 - Uniform depth
        LandslideProbability(grid, number_of_iterations=250,
        groundwater__depth_distribution='uniform', 
        groundwater__depth_min_value=0., (EDIT)
        groundwater__depth_max_value=2.)
    Option 6 - Lognormal recharge
        LandslideProbability(grid, number_of_iterations=250,
        groundwater__depth_distribution='lognormal', 
        groundwater__depth_mean=0.5.,
        groundwater__depth_standard_deviation=0.1)
    Option 7 - Lognormal_spatial recharge
        LandslideProbability(grid, number_of_iterations=250,
        groundwater__depth_distribution='lognormal_spatial', 
        groundwater__depth_mean=np.random.randint(0, 2, grid_size),
        groundwater__depth_standard_deviation=np.random.rand(grid_size))        
    Option 8 - Data_driven_spatial depth to water
        LandslideProbability(grid, number_of_iterations=250,
        groundwater__depth_distribution='data_driven_spatial', 
        groundwater__depth_HSD_inputs=[HSD_dict])

    Parameters
    ----------
    grid: RasterModelGrid
        A raster grid.
    number_of_iterations: int, optional
        Number of iterations to run Monte Carlo simulation (default=250).
    groundwater__recharge_distribution: str, optional
        single word indicating recharge distribution, either 'uniform',
        'lognormal', 'lognormal_spatial,' or 'data_driven_spatial'.
         (default='uniform')
    groundwater__recharge_min_value: float, optional (mm/d)
        minimum groundwater recharge for 'uniform' (default=20.)
    groundwater__recharge_max_value: float, optional (mm/d)
        maximum groundwater recharge for 'uniform' (default=120.)
    groundwater__recharge_mean: float, optional (mm/d) 
        mean grounwater recharge for 'lognormal'
        and 'lognormal_spatial' (default=None)
    groundwater__recharge_standard_deviation: float, optional (mm/d)
        standard deviation of grounwater recharge for 'lognormal'
        and 'lognormal_spatial' (default=None)
    groundwater__recharge_HSD_inputs: list, optional
        list of 3 dictionaries in order (default=[]) - HSD_dict {Hydrologic
        Source Domain (HSD) keys: recharge numpy array values}, {node IDs keys:
        list of HSD_Id values}, HSD_fractions {node IDS keys: list of
        HSD fractions values} (none) - for more information refer to
        Ref 1 & Ref 2 mentioned above, as this set of inputs require
        rigorous pre-processing of data.
    groundwater__depth_distribution: str, optional
        single word indicating groundwater depth distribution, either 'uniform',
        'lognormal', 'lognormal_spatial,' or 'data_driven_spatial'.
         (default='uniform')
    groundwater__depth_min_value: float, optional (m)
        minimum groundwater depth for 'uniform' (default=0.)
    groundwater__depth_max_value: float, optional (m)
        maximum groundwater depth for 'uniform' (default=2.)
    groundwater__depth_mean: float, optional (m) 
        mean groundwater depth for 'lognormal'
        and 'lognormal_spatial' (default=None)
    groundwater__depth_standard_deviation: float, optional (m)
        standard deviation of grounwater depth for 'lognormal'
        and 'lognormal_spatial' (default=None)
    groundwater__depth_HSD_inputs: list, optional
        list of 3 dictionaries in order (default=[]) - HSD_dict {node IDs keys:depth numpy array values}
    g: float, optional (m/sec^2)
        acceleration due to gravity.
    seed: int, optional
        seed for random number generation. if seed is assigned any value
        other than the default value of zero, it will create different
        sequence. To create a certain sequence repititively, use the same
        value as input for seed.

    Examples
    ----------
    >>> from landlab import RasterModelGrid
    >>> from landlab.components.landslides import LandslideProbability
    >>> import numpy as np
    
    Create a grid on which to calculate landslide probability.
    
    >>> grid = RasterModelGrid((5, 4), spacing=(0.2, 0.2))
    
    Check the number of core nodes.
    
    >>> grid.number_of_core_nodes
    6
    
    The grid will need some input data. To check the names of the fields
    that provide the input to this component, use the *input_var_names*
    class property.
    
    >>> sorted(LandslideProbability.input_var_names)  # doctest: +NORMALIZE_WHITESPACE
    ['soil__density',
     'soil__internal_friction_angle',
     'soil__maximum_total_cohesion',
     'soil__minimum_total_cohesion',
     'soil__mode_total_cohesion',
     'soil__saturated_hydraulic_conductivity',
     'soil__thickness',
     'soil__transmissivity',
     'topographic__slope',
     'topographic__specific_contributing_area']
    
    Check the units for the fields.
    
    >>> LandslideProbability.var_units('topographic__specific_contributing_area')
    'm'
    
    Create an input field.
    
    >>> grid.at_node['topographic__slope'] = np.random.rand(grid.number_of_nodes)
    
    If you are not sure about one of the input or output variables, you can
    get help for specific variables.
    
    >>> LandslideProbability.var_help('soil__transmissivity')  # doctest: +NORMALIZE_WHITESPACE
    name: soil__transmissivity
    description:
      mode rate of water transmitted through a unit width of saturated
      soil - either provided or calculated with Ksat and soil depth
    units: m2/day
    at: node
    intent: in
    
    Additional required fields for component.
    
    >>> scatter_dat = np.random.randint(1, 10, grid.number_of_nodes)
    >>> grid.at_node['topographic__specific_contributing_area'] = np.sort(
    ...      np.random.randint(30, 900, grid.number_of_nodes))
    >>> grid.at_node['soil__transmissivity'] = np.sort(
    ...      np.random.randint(5, 20, grid.number_of_nodes), -1)
    >>> grid.at_node['soil__saturated_hydraulic_conductivity'] = np.sort(
    ...      np.random.randint(2, 10, grid.number_of_nodes), -1) 
   >>> grid.at_node['soil__mode_total_cohesion'] = np.sort(
    ...      np.random.randint(30, 900, grid.number_of_nodes))
    >>> grid.at_node['soil__minimum_total_cohesion'] = (
    ...      grid.at_node['soil__mode_total_cohesion'] - scatter_dat)
    >>> grid.at_node['soil__maximum_total_cohesion'] = (
    ...      grid.at_node['soil__mode_total_cohesion'] + scatter_dat)
    >>> grid.at_node['soil__internal_friction_angle'] = np.sort(
    ...      np.random.randint(26, 40, grid.number_of_nodes))
    >>> grid.at_node['soil__thickness'] = np.sort(
    ...      np.random.randint(1, 3, grid.number_of_nodes))
    >>> grid.at_node['soil__density'] = (2000. * np.ones(grid.number_of_nodes))
    
    Instantiate the 'LandslideProbability' component to work on this grid,
    and run it.
    
    >>> ls_prob = LandslideProbability(grid)
    >>> np.allclose(grid.at_node['landslide__probability_of_failure'], 0.)
    True
    
    Run the *calculate_landslide_probability* method to update output
    variables with grid
    
    >>> ls_prob.calculate_landslide_probability()
    
    Check the output variable names.
    
    >>> sorted(ls_prob.output_var_names) # doctest: +NORMALIZE_WHITESPACE
    ['landslide__probability_of_failure',
     'soil__mean_relative_wetness',
     'soil__mean_watertable_depth',
     'soil__probability_of_saturation']
    
    Check the output from the component, including array at one node.
    
    >>> np.allclose(grid.at_node['landslide__probability_of_failure'], 0.)
    False
    >>> core_nodes = ls_prob.grid.core_nodes
    """

    # component name
    _name = 'Landslide Probability'
    __version__ = '2.0'
    # component requires these values to do its calculation, get from driver
    _input_var_names = (
        'topographic__specific_contributing_area',
        'topographic__slope',
        'soil__transmissivity',
        'soil__saturated_hydraulic_conductivity',
        'soil__mode_total_cohesion',
        'soil__minimum_total_cohesion',
        'soil__maximum_total_cohesion',
        'soil__internal_friction_angle',
        'soil__density',
        'soil__thickness',
        )

    #  component creates these output values
    _output_var_names = (
        'landslide__probability_of_failure',
        'soil__mean_relative_wetness',
        'soil__mean_watertable_depth',
        'soil__probability_of_saturation',
        )

    # units for each parameter and output
    _var_units = {
        'topographic__specific_contributing_area': 'm',
        'topographic__slope': 'tan theta',
        'soil__transmissivity': 'm2/day',
        'soil__saturated_hydraulic_conductivity': 'm/day',
        'soil__mode_total_cohesion': 'Pa or kg/m-s2',
        'soil__minimum_total_cohesion': 'Pa or kg/m-s2',
        'soil__maximum_total_cohesion': 'Pa or kg/m-s2',
        'soil__internal_friction_angle': 'degrees',
        'soil__density': 'kg/m3',
        'soil__thickness': 'm',
        'soil__mean_relative_wetness': 'None',
        'landslide__probability_of_failure': 'None',
        'soil__probability_of_saturation': 'None',
        'soil__depth_to_groundwater': 'm',
        'soil__mean_watertable_depth': 'm',
        }

    # grid centering of each field and variable
    _var_mapping = {
        'topographic__specific_contributing_area': 'node',
        'topographic__slope': 'node',
        'soil__transmissivity': 'node',
        'soil__saturated_hydraulic_conductivity': 'node',
        'soil__mode_total_cohesion': 'node',
        'soil__minimum_total_cohesion': 'node',
        'soil__maximum_total_cohesion': 'node',
        'soil__internal_friction_angle': 'node',
        'soil__density': 'node',
        'soil__thickness': 'node',
        'soil__mean_relative_wetness': 'node',
        'landslide__probability_of_failure': 'node',
        'soil__probability_of_saturation': 'node',
        'soil__depth_to_groundwater': 'node',
        'soil__mean_watertable_depth': 'node',
        }

    # short description of each field
    _var_doc = {
        'topographic__specific_contributing_area':
            ('specific contributing (upslope area/cell face )' +
             ' that drains to node'),
        'topographic__slope':
        'slope of surface at node represented by tan theta',
        'soil__transmissivity':
            ('mode rate of water transmitted' +
             ' through a unit width of saturated soil - ' +
              'either provided or calculated with Ksat ' + 
               'and soil depth'),
        'soil__saturated_hydraulic_conductivity':
            ('mode rate of water transmitted' +
             ' through soil - provided if transmissivity ' +
              'is NOT provided to calculate tranmissivity ' + 
               ' with soil depth'),
        'soil__mode_total_cohesion':
            'mode of combined root and soil cohesion at node',
        'soil__minimum_total_cohesion':
            'minimum of combined root and soil cohesion at node',
        'soil__maximum_total_cohesion':
            'maximum of combined root and soil cohesion at node',
        'soil__internal_friction_angle':
            ('critical angle just before failure' +
             ' due to friction between particles'),
        'soil__density': 'wet bulk density of soil',
        'soil__thickness': 'soil depth to restrictive layer',
        'soil__mean_relative_wetness':
            ('Indicator of soil wetness;' +
             ' relative depth perched water table' +
             ' within the soil layer'),
        'landslide__probability_of_failure':
            ('number of times FS is <=1 out of number of' +
             ' iterations user selected'),
        'soil__probability_of_saturation':
            ('number of times relative wetness is >=1 out of' +
             ' number of iterations user selected'),
        'soil__mean_watertable_depth':
            ('mean distance to groundwater table from distribution' +
             ' of depths to groundwater or saturated soils'),
        }

    # Run Component
    @use_file_name_or_kwds
    def __init__(self, grid, number_of_iterations=250,
                 groundwater__recharge_distribution=None,
                 groundwater__recharge_min_value=20.,
                 groundwater__recharge_max_value=120.,
                 groundwater__recharge_mean=None,
                 groundwater__recharge_standard_deviation=None,
                 groundwater__recharge_HSD_inputs=[],
                 groundwater__depth_distribution=None,
                 groundwater__depth_min_value=0.,
                 groundwater__depth_max_value=3.,
                 groundwater__depth_mean=None,
                 groundwater__depth_standard_deviation=None,
                 groundwater__depth_HSD_inputs=[],
                 seed=0, **kwds):

        """
        Parameters
        ----------
        grid: RasterModelGrid
            A raster grid.
        number_of_iterations: int, optional
            Number of iterations to run Monte Carlo simulation (default=250).
        groundwater__recharge_distribution: str, optional
            single word indicating recharge distribution, either 'uniform',
            'lognormal', 'lognormal_spatial,' or 'data_driven_spatial'.
             (default=None)
        groundwater__recharge_min_value: float, optional (mm/d)
            minium groundwater recharge for 'uniform' (default=20.)
        groundwater__recharge_max_value: float, optional (mm/d)
            maximum groundwater recharge for 'uniform' (default=120.)
        groundwater__recharge_mean: float, optional (mm/d) 
            mean groundwater recharge for 'lognormal'
            and 'lognormal_spatial' (default=None)
        groundwater__recharge_standard_deviation: float, optional (mm/d)
            standard deviation of groundwater recharge for 'lognormal'
            and 'lognormal_spatial' (default=None)
        groundwater__recharge_HSD_inputs: list, optional
            list of 3 dictionaries in order (default=[]) - HSD_dict
            {Hydrologic Source Domain (HSD) keys: recharge numpy array values},
            {node IDs keys: list of HSD_Id values}, HSD_fractions
            {node IDS keys: list of HSD fractions values} (none) - for more
            information refer to Ref 1 & Ref 2 mentioned above, as this set of
            inputs require rigorous pre-processing of data.
        groundwater__depth_distribution: str, optional
            single word indicating depth to water table distribution, either 
            'uniform', 'lognormal', 'lognormal_spatial,' or 
            'data_driven_spatial'.
             (default=None)
        groundwater__depth_min_value: float, optional (m)
            minium groundwater depth to water table for 'uniform' (default=0.)
        groundwater__depth_max_value: float, optional (m)
            maximum groundwater depth for 'uniform' (default=2.)
        groundwater__depth_mean: float, optional (m) 
            mean groundwater depth to water table for 'lognormal'
            and 'lognormal_spatial' (default=None)
        groundwater__depth_standard_deviation: float, optional (m)
            standard deviation of groundwater depth to water table for 
            'lognormal' and 'lognormal_spatial' (default=None)
        groundwater__depth_HSD_inputs: list, optional
            list of 3 dictionaries in order (default=[]) - HSD_dict
            {Hydrologic Source Domain (HSD = node) keys: groundwater depth distribution numpy array values}
        g: float, optional (m/sec^2)
            acceleration due to gravity.
        seed: int, optional
            seed for random number generation. if seed is assigned any value
            other than the default value of zero, it will create different
            sequence. To create a certain sequence repititively, use the same
            value as input for seed.
        """

        # Initialize seeded random number generation        
        self._seed_generator(seed)

        super(LandslideProbability, self).__init__(grid)
        print('Input read for groundwater__recharge_distribution')
        print(groundwater__recharge_distribution)
        print('Input read for  groundwater__depth_distribution')
        print(groundwater__depth_distribution)
        
        # Store grid and parameters and do unit conversions
        self.n = int(number_of_iterations)
        self._g = kwds.get('g', scipy.constants.g)
        self.groundwater__recharge_distribution = (
            groundwater__recharge_distribution)
        self.groundwater__depth_distribution = (
            groundwater__depth_distribution)
      
        # Following code will deal with the input distribution and associated
        # parameters
    
        # Recharge Uniform distribution
        if self.groundwater__recharge_distribution == 'uniform':
            self._recharge_min = groundwater__recharge_min_value
            self._recharge_max = groundwater__recharge_max_value
            self._Re = np.random.uniform(self._recharge_min, self._recharge_max,
                                        size=self.n)
            self._Re /= 1000. # Convert mm to m
        # Recharge Lognormal Distribution - Uniform in space
        elif self.groundwater__recharge_distribution == 'lognormal':
            assert (groundwater__recharge_mean != None), (
                'Input mean of the distribution!')
            assert (groundwater__recharge_standard_deviation != None), (
                'Input standard deviation of the distribution!')
            self._recharge_mean = groundwater__recharge_mean
            self._recharge_stdev = groundwater__recharge_standard_deviation
            self._mu_lognormal = np.log((self._recharge_mean**2)/np.sqrt(
                self._recharge_stdev**2 + self._recharge_mean**2))
            self._sigma_lognormal = np.sqrt(np.log((self._recharge_stdev**2)/(
                self._recharge_mean**2)+1))
            self._Re = np.random.lognormal(self._mu_lognormal,
                                          self._sigma_lognormal, self.n)
            self._Re /= 1000. # Convert mm to m
        # Recharge Lognormal Distribution - Variable in space                                  
        elif self.groundwater__recharge_distribution == 'lognormal_spatial':
            assert (groundwater__recharge_mean.shape[0] == (
                self.grid.number_of_nodes)), (
                'Input array should be of the length of grid.number_of_nodes!')
            assert (groundwater__recharge_standard_deviation.shape[0] == (
                self.grid.number_of_nodes)), (
                'Input array should be of the length of grid.number_of_nodes!')
            self._recharge_mean = groundwater__recharge_mean
            self._recharge_stdev = groundwater__recharge_standard_deviation
        # Recharge Custom HSD inputs - Hydrologic Source Domain -> Model Domain
        elif self.groundwater__recharge_distribution == 'data_driven_spatial':
            self._HSD_dict = groundwater__recharge_HSD_inputs[0]
            self._HSD_id_dict = groundwater__recharge_HSD_inputs[1]
            self._fract_dict = groundwater__recharge_HSD_inputs[2]
            self._interpolate_HSD_dict()

        # Depth to water table - Uniform distribution
        if self.groundwater__depth_distribution == 'uniform':
            self._depth_min = groundwater__depth_min_value
            self._depth_max = groundwater__depth_max_value
            # Test printing 2019/11/04
            print(self._depth_max)
            print(self._depth_min)
            
            self._De = np.random.uniform(self._depth_min, self._depth_max,
                                        size=self.n)
            print('Depth to Water Table Distribution')
            print(self._De)
            print('Mean Depth to Water Table')
            print(np.mean(self._De))
            
        # Depth to water table - Lognormal Distribution - Uniform in space
        elif self.groundwater__depth_distribution == 'lognormal':
            assert (groundwater__depth_mean != None), (
                'Input mean of the distribution!')
            assert (groundwater__depth_standard_deviation != None), (
                'Input standard deviation of the distribution!')
            self._depth_mean = groundwater__depth_mean
            self._depth_stdev = groundwater__depth_standard_deviation
            self._mu_lognormal = np.log((self._depth_mean**2)/np.sqrt(
                self._depth_stdev**2 + self._depth_mean**2))
            self._sigma_lognormal = np.sqrt(np.log((self._depth_stdev**2)/(
                self._depth_mean**2)+1))
            self._De = np.random.lognormal(self._mu_lognormal,
                                          self._sigma_lognormal, self.n)
        # Depth to water table - Lognormal Distribution - Variable in space                                  
        elif self.groundwater__depth_distribution == 'lognormal_spatial':
            assert (groundwater__depth_mean.shape[0] == (
                self.grid.number_of_nodes)), (
                'Input array should be of the length of grid.number_of_nodes!')
            assert (groundwater__depth_standard_deviation.shape[0] == (
                self.grid.number_of_nodes)), (
                'Input array should be of the length of grid.number_of_nodes!')
            self._depth_mean = groundwater__depth_mean
            self._depth_stdev = groundwater__depth_standard_deviation 
            
        # Depth to water table Custom HSD inputs - Hydrologic Source Domain 
        # -> Model Domain
        elif self.groundwater__depth_distribution == 'data_driven_spatial':
            self._HSD_dict = groundwater__depth_HSD_inputs[0]
            
        # Check if all input fields are initialized
        for name in self._input_var_names:
            if name not in self.grid.at_node:
                self.grid.add_zeros(name, at=self._var_mapping[name],
                                    units=self._var_units[name])

        # Check if all output fields are initialized
        for name in self._output_var_names:
            if name not in self.grid.at_node:
                self.grid.add_zeros(name, at=self._var_mapping[name],
                                    units=self._var_units[name])
        
        # Create a switch to imply whether Ksat is provided.
        if np.all(
            self.grid.at_node['soil__saturated_hydraulic_conductivity'] == 0):
            self.Ksat_provided = 0  # False
        else:
            self.Ksat_provided = 1  # True
        # Create a switch to imply whether min/max cohesion is provided.
        if np.all(
            self.grid.at_node['soil__minimum_total_cohesion'] == 0):
            self.Cmin_provided = 0  # False
        else:
            self.Cmin_provided = 1  # True
            
        self._nodal_values = self.grid.at_node

        # Raise an error if no grid provided
        if self.grid is None:
            raise ValueError('You must now provide an existing grid!')

    def calculate_factor_of_safety(self, i):

        """
        Method calculates factor-of-safety stability index by using
        node specific parameters, creating distributions of these parameters,
        and calculating the index by sampling these distributions 'n' times.

        The index is calculated from the 'infinite slope stabilty
        factor-of-safety equation' in the format of Pack RT, Tarboton DG,
        and Goodwin CN (1998),The SINMAP approach to terrain stability mapping.

        Parameters
        ----------
        i: int
            index of core node ID.
        """

        # generate distributions to sample from to provide input parameters
        # currently triangle distribution using mode, min, & max
        self._a = self.grid.at_node[
            'topographic__specific_contributing_area'][i]
        self._theta = self.grid.at_node['topographic__slope'][i]
        self._Tmode = self.grid.at_node['soil__transmissivity'][i]
        self._Ksatmode = self.grid.at_node[
            'soil__saturated_hydraulic_conductivity'][i]
        self._Cmode = self.grid.at_node['soil__mode_total_cohesion'][i]
        self._Cmin = self.grid.at_node['soil__minimum_total_cohesion'][i]
        self._Cmax = self.grid.at_node['soil__maximum_total_cohesion'][i]
        self._phi_mode = self.grid.at_node['soil__internal_friction_angle'][i]
        self._rho = self.grid.at_node['soil__density'][i]
        self._hs_mode = self.grid.at_node['soil__thickness'][i]
        
        if self.groundwater__recharge_distribution == 'data_driven_spatial':
            self._calculate_HSD_recharge(i)
            self._Re /= 1000.  # mm->m
        elif self.groundwater__recharge_distribution == 'lognormal_spatial':
            mu_lognormal = np.log((self._recharge_mean[i]**2)/np.sqrt(
                self._recharge_stdev[i]**2 + self._recharge_mean[i]**2))
            sigma_lognormal = np.sqrt(np.log((self._recharge_stdev[i]**2)/(
                self._recharge_mean[i]**2)+1))
            self._Re = np.random.lognormal(mu_lognormal,
                                          sigma_lognormal, self.n)
            self._Re /= 1000. # Convert mm to m
            
               
        # Depth to water table distribution based on distribution type
        if self.groundwater__depth_distribution == 'data_driven_spatial':
            self._calculate_HSD_groundwater_depth(i)
            print('i')
            print(i)
            print('self._HSD_dict')
            print(self._HSD_dict)
            print('self._De')
            print(self._De)
        elif self.groundwater__depth_distribution == 'lognormal_spatial':  
            mu_lognormal = np.log((self._depth_mean[i]**2)/np.sqrt(
                self._depth_stdev[i]**2 + self._depth_mean[i]**2))
            sigma_lognormal = np.sqrt(np.log((self._depth_stdev[i]**2)/(
                self._depth_mean[i]**2)+1))
            self._De = np.random.lognormal(mu_lognormal,
                                          sigma_lognormal, self.n)


        
        # Create a switch to imply whether Recharge or depth
        if self.groundwater__depth_distribution is not None:
            # depth of water  
            self._soil__mean_watertable_depth = np.mean(self._De)
            if self.groundwater__depth_distribution == 'data_driven_spatial':
                print('self._HSD_dict')
                print(self._HSD_dict)
            
            hw_dist= self._hs_mode - self._De
            print('Depth of Water Distribution (SoilDepth-DTWdist)')
            print(hw_dist)
            print('Depth distribution is not none - count saturated cells')
            # calculate probability of saturation
            countr = 0
            for val in hw_dist:            # find how many RW values >= 1
                if val <= 0:
                    countr = countr + 1  # number with RW values (>=1)
            # probability: No. high RW values/total No. of values (n)
            self._soil__probability_of_saturation = countr/self.n
              
            print('Soil Thickness This Node')
            print(self._hs_mode)
                         
            hw_dist[np.where(hw_dist<0)] = self._hs_mode  #replace saturated cells to be the max soil thickness
            print('Depth of Water Distribution (SoilDepth-DTWdist) hw>hs correction')
            print(hw_dist)
              
            print('Soil thickness All')
            print(self.grid.at_node['soil__thickness'])
            print('Mean water table depth')
            print(self._soil__mean_watertable_depth)


        if self.groundwater__recharge_distribution == 'data_driven_spatial':
            self._calculate_HSD_recharge(i)
            self._Re /= 1000.  # mm->m
        elif self.groundwater__recharge_distribution == 'lognormal_spatial':
            mu_lognormal = np.log((self._recharge_mean[i]**2)/np.sqrt(
                self._recharge_stdev[i]**2 + self._recharge_mean[i]**2))
            sigma_lognormal = np.sqrt(np.log((self._recharge_stdev[i]**2)/(
                self._recharge_mean[i]**2)+1))
            self._Re = np.random.lognormal(mu_lognormal,
                                          sigma_lognormal, self.n)
            self._Re /= 1000. # Convert mm to m
            
            #dummy value
           # self._De = self._Re*1000 # dummy variable
            
        # Depth to water table distribution based on distribution type
        if self.groundwater__depth_distribution == 'data_driven_spatial':
            self._calculate_HSD_groundwater_depth(i)
            
        elif self.groundwater__depth_distribution == 'lognormal_spatial':  
            mu_lognormal = np.log((self._depth_mean[i]**2)/np.sqrt(
                self._depth_stdev[i]**2 + self._depth_mean[i]**2))
            sigma_lognormal = np.sqrt(np.log((self._depth_stdev[i]**2)/(
                self._depth_mean[i]**2)+1))
            self._De = np.random.lognormal(mu_lognormal,
                                          sigma_lognormal, self.n)


        # Cohesion
        if self.Cmin_provided:
            self._C = np.random.triangular(self._Cmin, self._Cmode,
                                          self._Cmax, size=self.n)        
        else:
            Cmin = self._Cmode-0.3*self._Cmode
            Cmax = self._Cmode+0.3*self._Cmode
            self._C = np.random.triangular(self._Cmin, self._Cmode,
                                      self._Cmax, size=self.n)
        # phi - internal angle of friction provided in degrees
        phi_min = self._phi_mode-0.18*self._phi_mode
        phi_max = self._phi_mode+0.32*self._phi_mode
        self._phi = np.random.triangular(phi_min, self._phi_mode,
                                        phi_max, size=self.n)
        # soil thickness
        #hs_min = min(0.005, self._hs_mode-0.3*self._hs_mode) # Alternative
        hs_min = self._hs_mode-0.3*self._hs_mode 
        hs_max = self._hs_mode+0.1*self._hs_mode
        self._hs = np.random.triangular(hs_min, self._hs_mode,
                                       hs_max, size=self.n)
        self._hs[self._hs <= 0.] = 0.005
        # Hydraulic conductivity (Ksat)
        if self.Ksat_provided:
            Ksatmin = self._Ksatmode-(0.3*self._Ksatmode)
            Ksatmax = self._Ksatmode+(0.1*self._Ksatmode)
            self._Ksat = np.random.triangular(Ksatmin, self._Ksatmode, Ksatmax,
                                             size=self.n)
            self._T = self._Ksat*self._hs
        else:
            # Transmissivity (T)
            Tmin = self._Tmode-(0.3*self._Tmode)
            Tmax = self._Tmode+(0.1*self._Tmode)
            self._T = np.random.triangular(Tmin, self._Tmode, Tmax, size=self.n)

        # calculate Factor of Safety for n number of times
        # calculate components of FS equation
        # dimensionless cohesion
        self._C_dim = self._C/(self._hs*self._rho*self._g)
        
        #print('Before Calculate rw')
        #print('Soil thickness')
        #print(self._hs_mode)    
        #print('Depth of Water Distribution')
        #print(hw_dist)    
        #print('Relative wetness Distribution')
        #print(hw_dist/self._hs_mode)

        # relative wetness
        if self.groundwater__recharge_distribution == 'uniform':
            print('recharge uniform loop')
            self._rel_wetness = ((self._Re)/self._T)*(self._a/np.sin(
            np.arctan(self._theta)))
        elif self.groundwater__recharge_distribution == 'lognormal':
            self._rel_wetness = ((self._Re)/self._T)*(self._a/np.sin(
            np.arctan(self._theta)))
        elif self.groundwater__recharge_distribution == 'lognormal_spatial':
            self._rel_wetness = ((self._Re)/self._T)*(self._a/np.sin(
            np.arctan(self._theta)))
        elif self.groundwater__recharge_distribution == 'data_driven_spatial':
            self._rel_wetness = ((self._Re)/self._T)*(self._a/np.sin(
            np.arctan(self._theta)))
        elif self.groundwater__depth_distribution == 'uniform':
            self._rel_wetness = ((hw_dist)/self._hs_mode)   
            print('depth uniform loop')
        elif self.groundwater__depth_distribution == 'lognormal':
            self._rel_wetness = ((hw_dist)/self._hs_mode)       
        elif self.groundwater__depth_distribution == 'lognormal_spatial':
            self._rel_wetness = ((hw_dist)/self._hs_mode)          
        elif self.groundwater__depth_distribution == 'data_driven_spatial':
            self._rel_wetness = ((hw_dist)/self._hs_mode)
            
        print('self._rel_wetness')
        print(self._rel_wetness)
        
        
        # calculate probability of saturation
        if self.groundwater__recharge_distribution is not None:
            print('recharge  distribution is not none - count saturated cells')
            print('Relative wetness from recharge')
            print(self._rel_wetness)
            countr = 0
            for val in self._rel_wetness:            # find how many RW values >= 1
                if val >= 1.0:
                    countr = countr + 1  # number with RW values (>=1)
                
        # probability: No. high RW values/total No. of values (n)
        self._soil__probability_of_saturation = countr/self.n
        # Maximum Rel_wetness = 1.0        
        np.place(self._rel_wetness, self._rel_wetness > 1, 1.0)
        self._soil__mean_relative_wetness = np.mean(self._rel_wetness)
        print('self._soil__mean_relative_wetness')
        print(self._soil__mean_relative_wetness)
        print('self._rel_wetness')
        print(self._rel_wetness)
        print('self._phi')
        print(self._phi)
        Y = np.tan(np.radians(self._phi))*(1 - (self._rel_wetness*0.5))
        # convert from degrees; 0.5 = water to soil density ratio
        # calculate Factor-of-safety
        self._FS = (self._C_dim/np.sin(np.arctan(self._theta))) + (
            np.cos(np.arctan(self._theta)) *
            (Y/np.sin(np.arctan(self._theta))))
        count = 0
        for val in self._FS:                   # find how many FS values <= 1
            if val <= 1.0:
                count = count + 1   # number with unstable FS values (<=1)
        # probability: No. unstable values/total No. of values (n)
        self._landslide__probability_of_failure = np.array(count/self.n)

    def calculate_landslide_probability(self, **kwds):

        """
        Method creates arrays for output variables then loops through all
        the core nodes to run the method 'calculate_factor_of_safety.'
        Output parameters probability of failure, mean relative wetness,
        and probability of saturation are assigned as fields to nodes. 

        """

        # Create arrays for data with -9999 as default to store output
        self.mean_Relative_Wetness = np.full(self.grid.number_of_nodes, -9999.)
        self.prob_fail = np.full(self.grid.number_of_nodes, 0)
        self.prob_sat = np.full(self.grid.number_of_nodes, 0) 
        self.mean_watertable_depth = np.full(self.grid.number_of_nodes, -9999.)
        #self._hs_mode = np.full(self.grid.number_of_nodes, -9999.) 
        #self._De = np.full(self.grid.number_of_nodes, -9999.) 
        
        # Run factor of safety Monte Carlo for all core nodes in domain
        # i refers to each core node id
        for i in self.grid.core_nodes:
            self.calculate_factor_of_safety(i)
            # Populate storage arrays with calculated values
            self.mean_Relative_Wetness[i] = self._soil__mean_relative_wetness
            self.prob_fail[i] = self._landslide__probability_of_failure
            self.prob_sat[i] = self._soil__probability_of_saturation
            if self.groundwater__depth_distribution is not None:
                self.mean_watertable_depth[i]=self._soil__mean_watertable_depth

            
        # Values can't be negative
        self.mean_Relative_Wetness[
            self.mean_Relative_Wetness < 0.] = 0.
        self.prob_fail[self.prob_fail < 0.] = 0.
        
       
        
        # assign output fields to nodes
        self.grid.at_node['soil__mean_relative_wetness'] = (self.mean_Relative_Wetness)
        self.grid.at_node['landslide__probability_of_failure'] = self.prob_fail
        self.grid.at_node['soil__probability_of_saturation'] = self.prob_sat
        
        if self.groundwater__depth_distribution is not None:
            self.mean_watertable_depth[
                self.mean_watertable_depth < 0.] = 0.
            self.grid.at_node['soil__mean_watertable_depth'] = (
                self.mean_watertable_depth)
            
 
    def _seed_generator(self, seed=0):
        """Seed the random-number generator. This method will create the same
        sequence again by re-seeding with the same value (default value is
        zero). To create a sequence other than the default, assign non-zero
        value for seed.
        """
        np.random.seed(seed)


    def _interpolate_HSD_dict(self):
        """This method uses a non-parametric approach to expand the input
        recharge array to the length of number of iterations. Output is
        a new dictionary of interpolated recharge for each HSD id.
        """
        HSD_dict = copy.deepcopy(self._HSD_dict)
        # First generate interpolated Re for each HSD grid
        Yrand = np.sort(np.random.rand(self.n))
        # n random numbers (0 to 1) in a column
        for vkey in HSD_dict.keys():
            if isinstance(HSD_dict[vkey], int):
                continue       # loop back up if value is integer (e.g. -9999)
            Re_temp = HSD_dict[vkey]	 # an array of annual Re for 1 HSD grid
            Fx = ECDF(Re_temp)  # instantiate to get probabilities with Re
            Fx_ = Fx(Re_temp)    # probability array associated with Re data
            # interpolate function based on recharge data & probability
            f = interpolate.interp1d(Fx_, Re_temp, bounds_error=False,
                                     fill_value=min(Re_temp))
            # array of Re interpolated from Yrand probabilities (n count)
            Re_interpolated = f(Yrand)
            # replace values in HSD_dict with interpolated Re
            HSD_dict[vkey] = Re_interpolated

            self._interpolated_HSD_dict = HSD_dict
            
    def _interpolate_HSD_array(self):
        
        Yrand = np.sort(np.random.rand(self.n))
        Fx = ECDF(self._hwdist.values[0])
        Fx_ = Fx(self._hwdist.values[0])
        print(Fx_)
        f = interpolate.interp1d(Fx_, self._hwdist.values[0], bounds_error=False,
                                     fill_value=min(self._hwdist.values[0]))

        self._interp_De = f(Yrand)
        print('Input Depth Annual Year Max Data = 3 years')
        print(self._hwdist.values[0])
        print('Input Depth Annual Year Max Interpolated Data = n years')
        print(self._interp_De)


    def _calculate_HSD_recharge(self, i):
        """This method calculates the resultant recharge at node i of the
        model domain, using recharge of contributing HSD ids and the areal
        fractions of upstream contributing HSD ids. Output is a numpy array
        of recharge at node i.        
        """
        store_Re = np.zeros(self.n)
        HSD_id_list = self._HSD_id_dict[i]
        fract_list = self._fract_dict[i]
        for j in range(0, len(HSD_id_list)):
            Re_temp = self._interpolated_HSD_dict[HSD_id_list[j]]
            fract_temp = fract_list[j]
            Re_adj = (Re_temp*fract_temp)
            store_Re = np.vstack((store_Re, np.array(Re_adj)))
        self._Re = np.sum(store_Re, 0)
        
    def _calculate_HSD_groundwater_depth(self, i):
        """This method calculates the resultant groundwater depth at node i of the
        model domain, using depth to water table. Output is a numpy array
        of depth to groundwater at node i.        
        """
        
        groundwater__depth_HSD_inputs=self._HSD_dict[i]
        self._hwdist = pd.DataFrame.from_dict(groundwater__depth_HSD_inputs, orient='index')
        
        print('In calculate HSD loop - _hwdist.values[0]')
        print(self._hwdist.values[0])
        
        self._interpolate_HSD_array()
        
        print('In calculate HSD loop -interpolate _hwdist.values[0]')
        print(self._interp_De)
        
        self._De=self._interp_De
        
        print('interpolated Depth dist FOS input')
        print(self._De)