Joplin Tornado Example#

Please note that you might need additional dependencies to run this notebook, such as geopandas and contextily

  • Tornadoes occur at a high frequency in the United States compared with other natural hazards such as earthquakes and tsunamis but have a substantially smaller footprint. Even a single high-intensity tornado can result in high casualty rates and catastrophic economic losses and social consequences, particularly for small to medium communities.

  • The city of Joplin, Missouri, USA, was hit by an EF-5 tornado on May 22, 2011. The National Institute of Standards and Technology (NIST) conducted a technical investigation of this devastating event which can be found at: https://nvlpubs.nist.gov/nistpubs/NCSTAR/NIST.NCSTAR.3.pdf . The Center for Risk-Based Community Resilience Planning simulated this event for buildings and the electrical power network of Joplin in IN-CORE. This Juypter Notebook provides an example of how to use IN-CORE.

  • The initial damage prediction utilized the tornado path, tornado fragility curves representative of a 19- archetype building dataset, and EPN datasets. Generic tornado paths are also available in IN-CORE, or a user defined tornado path is possible.

  • The functionality of the infrastructure was linked with a computable general equilibrium (CGE) economics model that computes specific community resilience metrics in economic terms. A population dislocation model provides resilience metrics related to socio-demographics such as population dislocation as a function of income or race/ethnicity.

  • This example demonstrates how users interact with the IN-CORE computational environment.

*This notebook was created by Lisa Wang, supervised by Professor John W. van de Lindt, with the help of the NCSA team (Jong Sung Lee, Chris Navarro, Diego Calderon, Chen Wang, Michal Ondrejcek, Gowtham Naraharisetty, and Yong Wook Kim). The population disclocation model was developed by Nathanael Rosenheim and the CGE model portion provided by Brad Hartman under the supervision of Professor Harvey Cutler.

Prerequisites#

The following modules are necessary to run this notebook. To ensure dependencies are correct, install all modules through conda.

Module

Version

Notes

pyIncore

=>1.3.0

see: https://incore.ncsa.illinois.edu/doc/incore/install_pyincore.html

pyIncore_viz

=>1.5.0

see: https://incore.ncsa.illinois.edu/doc/pyincore_viz/index.html

matplotlib

3.1.2

used for plotting results

contextily

0.99.0

used for context geo-tiles

descartes

1.1.0

used for visualization of geometric objects as matplotlib paths

Economic Computable General Equilibrium (CGE) Model uses the ipopt solver. The ipopt solver is installed with pyIncore, so the environment should have the solver in it. If you have a trouble to run CGE model later, please see the instruction at CGE section in this document.

1. Hazard, Data, and Fragility#

This section introduces the input for the infrastructure damage analysis including the tornado path, building dataset, and building fragility curves for tornado.

1.1. Tornado Path#

This figure shows the 2011 Joplin tornado path with EF zones provided in the legend. The tornado path represents the wind speed within the vortex (multi-vortex in the case of Joplin) that was estimated to have EF5 wind speeds of more than 200 mph, reducing to EF4 wind speeds as the areas move outward from the vortex, and eventually reaching EF1 zone. (Attary et al. 2018)

from pyincore import IncoreClient
from pyincore import FragilityService, MappingSet, Dataset, DataService
from pyincore_viz.plotutil import PlotUtil as frag_plot
from pyincore_viz.geoutil import GeoUtil as viz
from pyincore.models.fragilitycurveset import FragilityCurveSet

client = IncoreClient()  # The authentication service validates username/password and  
                         # approves users to access to other services."https://incore-dev.ncsa.illinois.edu"
client.clear_cache()
tornado_path_id = "5d07cbeab9219c065b080930"
viz.plot_tornado(tornado_path_id, client, basemap=True)

1.2. Building Archetypes#

The 19 archetype buildings are used herein to represent the Joplin community. The selected building archetypes consider building characteristics such as footprint area, roof structure, number of stories, and so on. (Attary et al. 2018, Memari et al. 2018)

bldg_dataset_id = "5dbc8478b9219c06dd242c0d"  # Load the Joplin building dataset.

bldg_dataset = Dataset.from_data_service(bldg_dataset_id, DataService(client))

viz.plot_map(bldg_dataset, column='archetype',category='True')

Joplin building archetypes shown in the Legend above:

1 and 5 - Residential wood building, 6 - Business and retail building, 7 - Light industrial building, 8 - Heavy industrial building, 9 - Elementary/middle school, 10 - High school, 11 - Fire/police station, 12 - Hospital, 13 - Community center/church, 14 - Goverment building, 15 - Large big-box, 16 - Small big-box, 17 - Mobile home, 18 - Shopping center, 19 - Office building

1.3. Fragility Curves#

There 19 fragility curves set for these 19 building types in the four damage states, which covers the entire range of wind speeds associated with tornadoes (Attary et al. 2018, Memari et al. 2018). Below we selected 3 types to plot as examples.

fragility_service = FragilityService(client)
mapping_id = "5d8942dfb9219c068ea795ea"
mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))

# plot fragility for the first 3 archetypes using pyincore viz method
for mapping in mapping_set.mappings[:3]:
    fragility_id = mapping.entry['Non-Retrofit Fragility ID Code']
    fragility_set = FragilityCurveSet(fragility_service.get_dfr3_set(fragility_id))
    plt = frag_plot.get_fragility_plot(fragility_set, start=20, end=80)
    plt.show()

2. Infrastructure Damage Analysis#

The models in this section implement the prediction of community-level damage to the buildings and Electric Power Network (EPN) for the 2011 Joplin tornado.

hazard_type = "tornado"
hazard_id = "5d07cbeab9219c065b080930"  # The simulated EF-5 tornado shows geographical locations and the range
                                        # of wind speed of tornado hazard in Joplin.

2.1. Buildings Damage#

from pyincore.analyses.buildingdamage import BuildingDamage  # Import building damage module integrated into pyIncore.
mapping_id = "5d8942dfb9219c068ea795ea"       # Create a mapping to assign tornado fragilities to 19 building archetypes.

fragility_service = FragilityService(client)  # loading fragility mapping
mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))

bldg_dmg = BuildingDamage(client) 

building_inv = bldg_dmg.load_remote_input_dataset("buildings", bldg_dataset_id)
bldg_dmg.set_input_dataset("dfr3_mapping_set", mapping_set)
result_name = "Joplin_bldg_dmg_result"  
bldg_dmg.set_parameter("result_name", result_name)
bldg_dmg.set_parameter("hazard_type", hazard_type)
bldg_dmg.set_parameter("hazard_id", hazard_id)
bldg_dmg.set_parameter("num_cpu", 4)  # Define the result name, etc., and choose the number of CPU locally 
                                      # to run the model parallelly.
bldg_dmg.run_analysis()  # Run the building damage module to get building damage results for Joplin in a .csv file format. 
                         # The building damage results herein are referred to fragilities at three damage states (moderate, 
                         # heavy, complete) for 28152 buildings located in Joplin that fall within the tornado path or not.
building_dmg_result = bldg_dmg.get_output_dataset('ds_result')

# Convert dataset to Pandas DataFrame
df_bldg_dmg = building_dmg_result.get_dataframe_from_csv()
# Display top 5 rows of output data
df_bldg_dmg.head()

2.2. Electrical Power Facility (EPF) - Substations Damage#

from pyincore.analyses.epfdamage.epfdamage import EpfDamage # Import epf damage module integrated into pyIncore.
epf_substations_id = "5d92355bb9219c06ae7e386a"  
mapping_id = "5d8a326fb9219c068ea798e7"  # Create a mapping to assign tornado fragilities to substations.

fragility_service = FragilityService(client)
mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))

epf_sub_dmg = EpfDamage(client)  

epf_sub_dmg.load_remote_input_dataset("epfs", epf_substations_id) 
                               # Load the Joplin substations dataset, which is a package of GIS files.
epf_sub_dmg.set_input_dataset("dfr3_mapping_set", mapping_set)
result_name = "Joplin_epf_substations_dmg_result"
epf_sub_dmg.set_parameter("result_name", result_name)
epf_sub_dmg.set_parameter("hazard_type", hazard_type)
epf_sub_dmg.set_parameter("hazard_id", hazard_id)
epf_sub_dmg.set_parameter("num_cpu", 4)
epf_sub_dmg.set_parameter("fragility_key", "substations")  # Define the result name, etc., and choose the number
                                                           # of CPU locally to run the model parallelly.
epf_sub_dmg.run_analysis()  # Run the EPF damage module to get substations damage results for Joplin in a .csv file format. The 
                            # substations damage results herein are referred to fragilities at four damage states (insignificant, 
                            # moderate, extensive, complete) for 18 substations located in Joplin that fall within the tornado 
                            # path or not.
substation_dmg_result = epf_sub_dmg.get_output_dataset('result')

df_sub_dmg = substation_dmg_result.get_dataframe_from_csv()
df_sub_dmg.head()

2.3. Electrical Power Facility (EPF) - Poles Damage#

epf_poles_id = "5d923daab9219c06ae84afb0"  # Load the Joplin poles dataset, which is a package of GIS files.
mapping_id = "5d8a326fb9219c068ea798e7"    # Create a mapping to assign tornado fragilities to poles.

fragility_service = FragilityService(client)        # loading fragility mapping
mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))

epf_poles_dmg = EpfDamage(client) 

epf_poles_dmg.load_remote_input_dataset("epfs", epf_poles_id)

epf_poles_dmg.set_input_dataset("dfr3_mapping_set", mapping_set)
result_name = "Joplin_epf_poles_dmg_result"
epf_poles_dmg.set_parameter("result_name", result_name)
epf_poles_dmg.set_parameter("hazard_type", hazard_type)
epf_poles_dmg.set_parameter("hazard_id", hazard_id)
epf_poles_dmg.set_parameter("num_cpu", 4)
epf_poles_dmg.set_parameter("fragility_key", "poles")  # Define the result name, etc., and choose the number
                                                       # of CPU locally to run the model parallelly.
epf_poles_dmg.run_analysis()  # Run the EPF damage module to get poles damage results for Joplin in a .csv file 
                              # format. The poles damage results herein are referred to fragilities at four 
                              # damage states (insignificant, moderate, extensive, complete) for 23857 poles 
                              # located in Joplin that fall within the tornado path or not.
poles_dmg_result = epf_poles_dmg.get_output_dataset('result')


df_poles_dmg = poles_dmg_result.get_dataframe_from_csv()
df_poles_dmg.head()

3. Monte Carlo Simulation (MCS)#

Researchers can use Monte Carlo Simulation to estimate the probability of each building being in a particular damage state. This example uses 500 iterations to determine the failure probability of buildings reaching damage state 2, damage state 3, and damage state 4. Users can run 10000 samples or even more for a more accurate Monte Carlo Simulation to determine the building failure probabilities. Note that this takes several minutes and we are working on developing a more efficient algorithm.

num_samples = 500

3.1. MCS chaining with Joplin building damage#

from pyincore.analyses.montecarlofailureprobability import MonteCarloFailureProbability 
                                                  # Import Monte Carlo failure probability module integrated into pyIncore.
mc_bldg = MonteCarloFailureProbability(client)

mc_bldg.set_input_dataset("damage", building_dmg_result)  #  Load the Joplin building damage results dataset 
                                                         # generated from the previous model as an input
mc_bldg.set_parameter("num_cpu", 8)
mc_bldg.set_parameter("num_samples", num_samples)
mc_bldg.set_parameter("damage_interval_keys", ["DS_0", "DS_1", "DS_2", "DS_3"])
mc_bldg.set_parameter("failure_state_keys", ["DS_1", "DS_2", "DS_3"])

mc_bldg.set_parameter("result_name", "tornado_mc_failure_probability_buildings") # name of csv file with results
mc_bldg.run_analysis()  # Run the Monte Carlo Simulation module to obtain the building failure probabilities. The building failure
                       # probabilities herein only consider the physical damage without the interdependency.

building_failure_probability = mc_bldg.get_output_dataset('failure_probability')  # get buildings failure probabilities

df_bldg_fail = building_failure_probability.get_dataframe_from_csv()
df_bldg_fail.head()
building_damage_mcs_samples = mc_bldg.get_output_dataset('sample_failure_state')  # get buildings failure states

bdmcs = building_damage_mcs_samples.get_dataframe_from_csv()
bdmcs.head()

3.2. MCS chaining with Joplin EPF substations damage#

mc_sub = MonteCarloFailureProbability(client)

mc_sub.set_input_dataset("damage", substation_dmg_result)
mc_sub.set_parameter("num_cpu", 8)
mc_sub.set_parameter("num_samples", num_samples)
mc_sub.set_parameter("damage_interval_keys", ["DS_0", "DS_1", "DS_2", "DS_3", "DS_4"])
mc_sub.set_parameter("failure_state_keys", ["DS_1", "DS_2", "DS_3", "DS_4"])

mc_sub.set_parameter("result_name", "Joplin_mcs_substations_samples") # name of csv file with results
mc_sub.run_analysis()

substation_damage_mcs_samples = mc_sub.get_output_dataset('sample_failure_state')

sdmcs = substation_damage_mcs_samples.get_dataframe_from_csv()
sdmcs.head()

3.3. MCS chaining with Joplin EPF poles damage#

mc_pole = MonteCarloFailureProbability(client)

mc_pole.set_input_dataset("damage", poles_dmg_result)
mc_pole.set_parameter("num_cpu", 8)
mc_pole.set_parameter("num_samples", num_samples)
mc_pole.set_parameter("damage_interval_keys", ["DS_0", "DS_1", "DS_2", "DS_3", "DS_4"])
mc_pole.set_parameter("failure_state_keys", ["DS_1", "DS_2", "DS_3", "DS_4"])

mc_pole.set_parameter("result_name", "Joplin_mcs_poles_samples") # name of csv file with results
mc_pole.run_analysis()

pole_damage_mcs_samples = mc_pole.get_output_dataset('sample_failure_state')

pdmcs = pole_damage_mcs_samples.get_dataframe_from_csv()
pdmcs.head()

3.4. Building Damage Spatial Distribution Results#

Joplin building failure probability map. Legend shows failure probability from 0 to 1 (total) failure.

# getting geodataframework of building dataset and merge with output
bldg_gdf = bldg_dataset.get_dataframe_from_shapefile()
bldg_fail_gdf = bldg_gdf.merge(df_bldg_fail, on='guid')
viz.plot_gdf_map(bldg_fail_gdf, column='failure_probability')

4. Infrastructure Functionality Analysis#

The functionality analysis module below can be used to calculate building functionality probabilities considering two situations: buildings are in at least a damage state 2 or greater or buildings are not damaged but electric power is not available to the building. Whether buildings can receive electrical power is assumed to depend on the interdependency between buildings and substations, and between buildings and poles in close proximity. If both the nearest pole to the building and the substation where buildings belong to its service area are functional, buildings are considered to be able to receive electric power.

4.1. Building Functionality Probability#

from pyincore.analyses.buildingfunctionality import BuildingFunctionality
                                            # Import building functionality module integrated into pyIncore.
bldg_func = BuildingFunctionality(client)

# Load the datasets of building samples, substation samples, and pole samples. All of the samples are randomly 
# generated using the Monte Carlo Simulation module introduced in the previous section.
bldg_func.set_input_dataset("building_damage_mcs_samples", building_damage_mcs_samples)
bldg_func.set_input_dataset("substations_damage_mcs_samples", substation_damage_mcs_samples)
bldg_func.set_input_dataset("poles_damage_mcs_samples", pole_damage_mcs_samples)

# Load the dataset of the interdependency table between buildings and substations, and between buildings and poles.
interdependency_id = "5dcf4a34b9219ca5e4118312"
bldg_func.load_remote_input_dataset("interdependency_dictionary", interdependency_id)

bldg_func.set_parameter("result_name", "Joplin_mcs_functionality_probability")
bldg_func.run_analysis()  # Run the module to obtain the building functionality probabilities.

bldg_func_samples = bldg_func.get_output_dataset('functionality_samples')
df_bldg_samples = bldg_func_samples.get_dataframe_from_csv()
df_bldg_samples.head()
bldg_func_probability = bldg_func.get_output_dataset('functionality_probability')
df_bldg_func = bldg_func_probability.get_dataframe_from_csv()
df_bldg_func = df_bldg_func.rename(columns={"building_guid": "guid"})
df_bldg_func.head()

4.2. Building Functionality Spatial Distribution Results#

Joplin building functionality map.

# merge building geodataframe with output
bldg_func_gdf = bldg_gdf.merge(df_bldg_func, on='guid')
viz.plot_gdf_map(bldg_func_gdf, column='probability')

4.3. Physical Service Resilience Metrics#

After finding the damage level for each component (buildings, substations, and poles) based on the components’ fragility curves, their intrinsic failure status is expressed as a binary format with either failed (0) or not-failed (1). Then the failure status of all the buildings is updated by considering their dependencies with the corresponding electric power facilities. Each component generates 500 samples randomly with their failure status determined, and the percentage of building functional and nonfunctional could be calculated using the updated status of the same sample (such as sample #1) for all the buildings, illustrated herein as an example.

The percent shown in the horizontal bar below refers to the number of nonfunctional/damaged buildings divided by all buildings in Joplin. For example, the mean percentage of building nonfunctional is around 52%, and the standard deviation percentage shown as an error bar. The percentage is two times more than the percent of building in damage is because three substations fall within the tornado path. Consequently, there is a significant probability for buildings located in their service area cannot be available to access the electric power even though some buildings themselves are still safe and operational in the immediate response of the tornado hazard. It is worthy to note that a random tornado path analysis will produce a larger standard deviation of building damage and building functionality.

Calculate and show all three infrustructure damages and resulting total building funcionality.

import pandas as pd
import numpy as np

# buildings damage
bdm = bdmcs["failure"].str.split(pat=",")
# strings to integers, change fail to 1
bdm = 1 - bdm.apply(pd.to_numeric)
# convert panda frame to numpy ndarray and do the statistics
bdm = bdm.to_numpy()
# convert list of np arrays into single one
bdm = np.stack(bdm)
# percent
bdm_guid = 100 * bdm.mean(axis=0)
bdm_mean = np.mean(bdm_guid)
bdm_std = np.std(bdm_guid)

# substations damage
sdm = sdmcs["failure"].str.split(pat=",")
sdm = 1 - sdm.apply(pd.to_numeric)
sdm = sdm.to_numpy()
sdm = np.stack(sdm)
sdm_guid = 100 * sdm.mean(axis=0)
sdm_mean = np.mean(sdm_guid)
sdm_std = np.std(sdm_guid)

# poles damage
pdm = pdmcs["failure"].str.split(pat=",")
pdm = 1 - pdm.apply(pd.to_numeric)
pdm = pdm.to_numpy()
pdm = np.stack(pdm)
pdm_guid = 100 * pdm.mean(axis=0)
pdm_mean = np.mean(pdm_guid)
pdm_std = np.std(pdm_guid)

# buildings combined functionality
bfm = df_bldg_samples["samples"].str.split(pat=",")
bfm = 1 - bfm.apply(pd.to_numeric)
bfm = bfm.to_numpy()
bfm = np.stack(bfm)
bfm_guid = 100 * bfm.mean(axis=0)
bfm_mean = np.mean(bfm_guid)
bfm_std = np.std(bfm_guid)
# dataframe from statistical results
infra = ('Percent of building nonfunctional', 'Percent of building in damage', 'Percent of substation in damage',\
         'Percent of pole in damage')
res_data = {"infrastructure": infra,
            "percentage": [np.round(bfm_mean, 2), np.round(bdm_mean, 2), 
                           np.round(sdm_mean, 2), np.round(pdm_mean, 2)],
            "std": [np.round(bfm_std, 2), np.round(bdm_std, 2), 
                    np.round(sdm_std, 2), np.round(pdm_std, 2)]}
func_data = pd.DataFrame(res_data)
print(func_data)
plt.rcdefaults()
fig, ax = plt.subplots()

y_pos = np.arange(len(infra))
performance = res_data['percentage']

for i, v in enumerate(performance):
    ax.text(v + 2, i, str(v), color='black', fontweight='bold')
    
ax.barh(y_pos, performance, xerr=res_data['std'], align='center', height=0.6)
ax.set_yticks(y_pos)
ax.set_yticklabels(infra)
ax.invert_yaxis()
ax.set_xlabel('Failures to Infrastructure')
ax.set_title('Physical Service Resilience Metrics Percentage Change',fontsize = 12,weight='bold')
plt.xlim(0, 100)
plt.show()  # Physical service resilience metrics percentage change immediately after the simulated EF-5 tornado.

5. Economic analysis#

5.1. Computable General Equilibrium (CGE) Model#

A computable general equilibrium (CGE) model is based on fundamental economic principles. A CGE model uses multiple data sources to model the interactions of households, firms and relevant government entities as they contribute to economic activity. The model is based on (1) utility-maximizing households that supply labor and capital, using the proceeds to pay for goods and services (both locally produced and imported) and taxes; (2) the production sector, with perfectly competitive, profit-maximizing firms using intermediate inputs, capital, land and labor to produce goods and services for both domestic consumption and export; (3) the government sector that collects taxes and uses tax revenues in order to finance the provision of public services; and (4) the rest of the world.

5.1.2. Capital stock shock analysis#

The first step is to aggregate building functionality states and calculates total capital shock losses per sector. We use results of building failure probabilities obtained in previous MCS chaining with Joplin building damage section.

from pyincore.analyses.capitalshocks import CapitalShocks
# Joplin building to sector mapping table
building_to_sectors_id = "5f202d674620b643d787a5e7"
# Create Capital shocks analysis
capital_shocks = CapitalShocks(client)

# Load remote datasets
capital_shocks.load_remote_input_dataset("buildings_to_sectors", building_to_sectors_id)
# Set datasets
# Joplin building inventory
capital_shocks.load_remote_input_dataset("buildings", bldg_dataset_id)
# Joplin building failure probability
capital_shocks.set_input_dataset("failure_probability", building_failure_probability)

capital_shocks.set_parameter("result_name", "sector_shocks") # name of csv file with results
# Run capital shocks analysis
capital_shocks.run_analysis()

sector_shocks_result = capital_shocks.get_output_dataset("sector_shocks")
sector_shocks_result.get_dataframe_from_csv()

5.1.3. CGE Model#

from pyincore.analyses.joplincge import JoplinCGEModel

Each of the objects with string IDs below refer to datasets used as inputs in the CGE model.

# SAM
# Social accounting matrix (SAM) contains data for firms, households and government which are organized 
# in a way to represent the interactions of all three entities in a typical economy
SAM = "5dd85ae7b9219c06d4da8de4"

# CAPITAL COMP
# BB is a matrix which describes how investment in physical infrastructure is transformed into functioning capital such as commercial and residential buildings. 
# These data are collected from the Bureau of Economic Analysis (BEA).
BB = "5dc1e620b9219c06dd2f473a"

# MISC TABLES
IOUT = "5dc1e6d8b9219c06dd2f475e"  # This is matrix that describes the transfer of tax revenue collected by the local government to help finance local government expenditures. 
MISC = "5dc1e736b9219c06dd2f4782"  # This is the name of a file that contains data for commercial sector employment and physical capital. It also contains data for the number of households and working households in the economy. 
MISCH = "5dc1e7b5b9219c06dd2f47a6"  # A file that contains elasticities for the supply of labor with respect to paying income taxes.
LANDCAP = "5dc1e810b9219c06dd2f47ca"  # Contains information regarding elasticity values for the response of changes in the price of physical capital with respect to the supply of investment.
EMPLOY = "5dc1e85ab9219c06dd2f47ee"  # Table name containing data for commercial sector employment.
IGTD = "5dc1e895b9219c06dd2f4812"  # This variable represents a matrix describing the transfer of taxes collected to a variable which permits governments to spend the tax revenue on workers and intermediate inputs.
TAUFF = "5dc1e8eeb9219c06dd2f4836"  # Represents social security tax rates
JOBCR = "5dc1e962b9219c06dd2f487e"  # This is a matrix describing the supply of workers coming from each household group in the economy. 
OUTCR = "5dc1e9aeb9219c06dd2f48bc"  # This a matrix describing the number of workers who live in Joplin but commute outside of town to work.
sector_shocks = "5f21d40d4620b643d78bb4c2"# This is the aggregation of building functionality states to capital shocks per sector.
# Create Joplin CGE Model
joplin_cge = JoplinCGEModel(client)
# Load analysis input datasets
joplin_cge.load_remote_input_dataset("SAM", SAM)
joplin_cge.load_remote_input_dataset("BB", BB)
joplin_cge.load_remote_input_dataset("IOUT", IOUT)
joplin_cge.load_remote_input_dataset("MISC", MISC)
joplin_cge.load_remote_input_dataset("MISCH", MISCH)
joplin_cge.load_remote_input_dataset("LANDCAP", LANDCAP)
joplin_cge.load_remote_input_dataset("EMPLOY", EMPLOY)
joplin_cge.load_remote_input_dataset("IGTD", IGTD)
joplin_cge.load_remote_input_dataset("TAUFF", TAUFF)
joplin_cge.load_remote_input_dataset("JOBCR", JOBCR)
joplin_cge.load_remote_input_dataset("OUTCR", OUTCR)

# Set analysis input dataset from previous Capital stock shock analysis
joplin_cge.set_input_dataset("sector_shocks", sector_shocks_result)

The ipopt solver is installed with pyIncore, so the environment should have the solver in it, but the solver_path variable lets the user specify a version elsewhere when desired. The variable points to the absolute path for where the ipopt solver is installed. To know the absolute path of the ipopt solver in your local machine, open a terminal window, activate your pyIncore environment, and run the following command according to your OS:

  • on Mac: open a terminal window, activate your conda environment and type “which ipopt”.

  • on Windows: open a command terminal, activate your conda environment and type “where ipopt”.

and make sure that the folder containing “ipopt” is in PATH environment of your OS.

# Set analysis parameters
joplin_cge.set_parameter("model_iterations", 1)
# Run Joplin CGE model analysis
joplin_cge.run_analysis()
# Retrieve result datasets
domestic_supply = joplin_cge.get_output_dataset('domestic-supply')
ds = domestic_supply.get_dataframe_from_csv()
ds.head(6)
gross_income = joplin_cge.get_output_dataset('gross-income')
gi = gross_income.get_dataframe_from_csv()
gi.head()
pre_factor_demand = joplin_cge.get_output_dataset('pre-disaster-factor-demand')
pre_fd = pre_factor_demand.get_dataframe_from_csv()
pre_fd.head()
post_factor_demand = joplin_cge.get_output_dataset('post-disaster-factor-demand')
pos_fd = post_factor_demand.get_dataframe_from_csv()
pos_fd.head()
household_count = joplin_cge.get_output_dataset('household-count')
hc = household_count.get_dataframe_from_csv()
hc.head()

5.2. Economic Resilience Metrics#

5.2.1. Percent Reduction of Domestic Supply#

plt.rcdefaults()
fig_ds, ax_ds = plt.subplots()

sectors_ds = ('Goods', 'Trades', 'Others', 'Housing service 1', 'Housing service 2', 'Housing service 3')
y_pos_ds = np.arange(len(sectors_ds))
perf_ds = (1 - ds['DSL'] / ds['DS0']) * 100  # percentage

for i, v in enumerate(round(perf_ds, 2)):
    ax_ds.text(v + 2, i, str(v), color='black', fontweight='bold')
ax_ds.barh(y_pos_ds, perf_ds, align='center', height=0.8)
ax_ds.set_yticks(y_pos_ds)
ax_ds.set_yticklabels(sectors_ds)
ax_ds.invert_yaxis()
ax_ds.set_title('Economic Resilience Metrics Percentage Change', fontsize=12, weight='bold')
ax_ds.set_xlabel('Percent reduction of domestic supply (%)')
plt.xlim(0, 50)
fig_ds.tight_layout()
plt.show()

5.2.2. Factor demand percent reduction before and after disaster#

# Labor groups
data_pre = {"Sectors": ["L1", "L2", "L3"],
            "GOODS": [449.006038, 3282.049564, 3013.051661],
            "TRADE": [2729.055379, 4360.095686, 1851.044407],
            "OTHER": [6306.146353, 12396.308192, 5444.146479]}
pre_fd = pd.DataFrame(data_pre)
data_pos = {"Sectors": ["L1", "L2", "L3"],
            "GOODS": [438.913649, 3195.254485, 2919.538642],
            "TRADE": [2586.192593, 4115.076833, 1738.785852],
            "OTHER": [6026.389334, 11798.285962, 5157.077844]}
pos_fd = pd.DataFrame(data_pos)

perf_gd = (1 - pos_fd['GOODS'] / pre_fd['GOODS']) * 100  # percentage
perf_tr = (1 - pos_fd['TRADE'] / pre_fd['TRADE']) * 100
perf_ot = (1 - pos_fd['OTHER'] / pre_fd['OTHER']) * 100

# Labor groups
labor_groups = ('L1', 'L2', 'L3')
goods = perf_gd.values.round(2).tolist()
trade = perf_tr.values.round(2).tolist()
other = perf_ot.values.round(2).tolist()

y_pos_fd = np.arange(len(labor_groups))

plt.rcdefaults()
df_all = pd.DataFrame({'goods': goods, 'trade': trade, 'other': other}, index = labor_groups)
ax_fd = df_all.plot.barh(align='center')

for i, v in enumerate(round(perf_gd, 2)):
    ax_fd.text(v + 0.3, i - 0.15, str(v), color='black', fontweight='bold')
for i, v in enumerate(round(perf_tr, 2)):
    ax_fd.text(v + 0.3, i + 0.05, str(v), color='black', fontweight='bold')
for i, v in enumerate(round(perf_ot, 2)):
    ax_fd.text(v + 0.3, i + 0.25, str(v), color='black', fontweight='bold')
ax_fd.set_yticks(y_pos_fd)
ax_fd.set_yticklabels(labor_groups)
ax_fd.invert_yaxis()
ax_fd.set_xlabel('Factor demand reduction (%)')
ax_fd.set_title('Factor demand after disaster', fontsize=12, weight='bold')
plt.xlim(0, 10)
plt.show()

5.2.3. Percent Reduction of Employment#

plt.rcdefaults()
fig_hc, ax_hc = plt.subplots()
household_hc = ('HH1 (less than $15,000)', 'HH2 ($15,000 to 24,999)', 'HH3 ($25,000 to 74,999)',
                'HH4 ($75,000 to 99,999)', 'HH5 ($100,000 or more)')
y_pos = np.arange(len(household_hc))
perf_hc = (1 - hc['HHL'] / hc['HH0']) * 100  # percentage

for i, v in enumerate(round(perf_hc, 2)):
    ax_hc.text(v + 0.3, i, str(v), color='black', fontweight='bold')
ax_hc.barh(y_pos, perf_hc, align='center', height=0.6)
ax_hc.set_yticks(y_pos)
ax_hc.set_yticklabels(household_hc)
ax_hc.invert_yaxis()
ax_hc.set_title('Economic Resilience Metrics Percentage Change', fontsize=12, weight='bold')
ax_hc.set_xlabel('Percent reduction of employment (%)')
plt.xlim(0, 15)
fig_hc.tight_layout()
plt.show()

6. Sociological Analysis#

This section introduces the sociology-based analysis which may be integrated with hazard and building damage analysis. In this example the damage to buildings drives the population dislocation algorithm at the household level.

6.1. Population Dislocation Model#

The population dislocation model depends first on the allocation of detailed housing unit and household characteristic data to each residential building within the community. The allocation of detailed household characteristic data provides an estimate of the number of people that live within each structure. The dislocation model is based on household surveys, for a more detailed summary of the dislocation model please refer to Rosenheim, Guidotti, Gardoni, and Peacock 2019.

The model predicts which households within the community will be dislocated from their homes immediately following the hazard event.

from pyincore.analyses.populationdislocation import PopulationDislocation, PopulationDislocationUtil
housing_unit_alloc = "5dc1c196b9219c06dd2e3f0b"
bg_data = "5d4c9545b9219c0689b2358a"
value_loss = "60354810e379f22e16560dbd"
# Create Population dislocation analysis
pop_dis = PopulationDislocation(client)
# Load analysis input datasets
pop_dis.load_remote_input_dataset("housing_unit_allocation", housing_unit_alloc)
pop_dis.load_remote_input_dataset("block_group_data", bg_data)
pop_dis.load_remote_input_dataset("value_loss_param", value_loss)

pop_dis.set_input_dataset("building_dmg", building_dmg_result)  # Load the Joplin building damage results dataset 
                                                                # generated from the previous model as an input
result_name = "pop-dislocation-results"
seed = 1111

pop_dis.set_parameter("result_name", result_name)
pop_dis.set_parameter("seed", seed)
pop_dis.run_analysis()

# Retrieve result dataset
result = pop_dis.get_output_dataset("result")

# Convert dataset to Pandas DataFrame
df = result.get_dataframe_from_csv(low_memory=False)
df.head()

6.2. Sociological Resilience Metrics#

6.2.1. Percent of Population Dislocation#

The PopulationDislocation model predicts the probability of dislocation using a logistic regression model based on household level surveys from Hurricane Andrew in 1992, which was a major wind event. The logistic regression equation provides a probability of dislocation based on the value loss and the other factors. For more details on the population dislocation algorithm and the logistic regression model see Rosenheim et al 2019.

To determine if a household dislocate a random value was uniformly sampled between 0 and 1: if this value was lower than the probability of dislocation, then the household was set to dislocate. The results below provide one example of how the output can be explored by population subgroups. The model is designed to run as part of a larger Monte Carlo Simulation.

result = pop_dis.get_output_dataset("result")

df = result.get_dataframe_from_csv(low_memory=False)
# Select Housing units allocated to buildings
table_title = "Total population by tenure status and dislocation."
pd.crosstab(df['ownershp'], df['dislocated'], df['numprec'], aggfunc = sum,
            dropna=True, margins=True, margins_name="Total").style.set_caption(table_title)
table_title = "Percent population by tenure status and dislocation."
pd.crosstab(df['ownershp'], df['dislocated'], df['numprec'], aggfunc = sum, normalize='index',
            dropna=True, margins=True, margins_name="Total").round(4)*100

The variable ownershp represents the tenure status of households. When ownershp = 1 the housing unit is owner-occupied. When the ownershp = 2 the housing unit is renter-occupied. The results above provide the percent of the population predicted to dislocate based on tenure status. When the variable dislocated = True the household is predicted to dislocate.

7. References#

Attary, N., van de Lindt, J. W., Mahmoud, H., Smith, S., Navarro, C. M., Kim, Y. W., & Lee, J. S. (2018). Hindcasting community-level building damage for the 2011 Joplin EF5 tornado. Natural Hazards, 93(3), 1295-1316.

Attary, N., van de Lindt, J. W., Mahmoud, H., & Smith, S. (2018). Hindcasting Community-Level Damage to the Interdependent Buildings and Electric Power Network after the 2011 Joplin, Missouri, Tornado. Natural Hazards Review, 20(1), 04018027.

Memari, M., Attary, N., Masoomi, H., Mahmoud, H., van de Lindt, J. W., Pilkington, S. F., & Ameri, M. R. (2018). Minimal building fragility portfolio for damage assessment of communities subjected to tornadoes. Journal of Structural Engineering, 144(7), 04018072.

Guidotti, R., Gardoni, P., & Rosenheim, N. (2019). Integration of physical infrastructure and social systems in communities’ reliability and resilience analysis. Reliability Engineering & System Safety, 185, 476-492. https://doi.org/10.1016/j.ress.2019.01.008

Rosenheim, N., Guidotti, R., Gardoni, P., & Peacock, W. (2019). Integration of detailed household and housing unit characteristic data with critical infrastructure for post-hazard resilience modeling. Sustainable and Resilient Infrastructure. https://doi.org/10.1080/23789689.2019.1681821