<jupyter_start><jupyter_text>Economic infrastructure intersections  Table content:   1) Set-up - importing libraries and specifying data folder<jupyter_code>import os
import sys
import geopandas
import pandas
%matplotlib notebook
import matplotlib.pyplot as plt
data_folder = "/Users/lena/OneDrive - OnTheHub - The University of Oxford/University/Oxford/OX_2018_2019/stlucia/analysis/python/projects/nismod-int-hack/data/Hackathon/rv"<jupyter_output><empty_output><jupyter_text>Loading in all hazards and all economic sectors as shapefiles<jupyter_code>def load_hazard(data_folder, hazard_id):  
    hazard_path = os.path.join(
        data_folder, 'data', 'hazards', '{}.shp'.format(hazard_id))
    hazards = geopandas.read_file(hazard_path)
    
    if hazards.crs != {'init':'epsg:2006'}:
        hazards = hazards.to_crs({'init':'epsg:2006'})
    return hazards
def load_sector(data_folder, sector_id):  
    sector_path = os.path.join(
        data_folder, 'data', 'infrastructure', '{}.shp'.format(sector_id))
    sectors = geopandas.read_file(sector_path)
    
    if sectors.crs != {'init':'epsg:2006'}:
        sectors = sectors.to_crs({'init':'epsg:2006'})
    return sectors<jupyter_output><empty_output><jupyter_text>Loading in merged economic infrastructures as csv and admin<jupyter_code>#INPUT FILE
administrative_path = os.path.join(
    data_folder, 'data','admin', 'stluc_administrative.shp')

administrative = geopandas.read_file(administrative_path)

#OUTPUT FILE
econ_merged_path= os.path.join(
    data_folder, 'data', 'infrastructure', 'econ_merged_admin.csv')

intersections_econ_path = os.path.join(
    data_folder, 'results', 'intersections_econ_merged.csv')

econ_haz_cap_path = os.path.join(
    data_folder, 'results', 'econ_hazards_capacity.csv')

econ_haz_percent_exp_path = os.path.join(
    data_folder, 'results','econ_hazards_sector%_exp.csv')

econ_haz_count_path = os.path.join(
    data_folder, 'results','econ_hazards_sector_count_exp.csv')

econ_multi_haz_asset_path = os.path.join(
    data_folder, 'results','econ_multi_hazards_assets_exp.csv')


econ_haz_abs_path = os.path.join(
    data_folder, 'results','econ_hazards_sector_abs_exp.csv')

SDG_cap_path = os.path.join(
    data_folder, 'results','econ_hazards_SDG_exp.csv')

roads_haz_percent_exp_path = os.path.join(
    data_folder, 'results', 'roads_hazards_%_exp.csv')
sector_ids = ['wastewater', 'waste','electricity', 'freight', 'airports', 'port', 'water']
hazard_ids = ['1m_sea-level','4m_storm-surge','flashflooding','landslide_susceptibility']
#  Adding multihazard computation
def two_multi_hazards(row):
    if row['4m_storm-surge'] > 0:
        val = row['4m_storm-surge']
    elif row['flashflooding'] > 0:
        val = row['flashflooding']
    else:
        val = 0 
    return val

def three_multi_hazards(row):
    if row['4m_storm-surge'] > 0:
        val = row['4m_storm-surge']
    elif row['flashflooding'] > 0:
        val = row['flashflooding']
    elif row['landslide_susceptibility'] > 0:
        val = row['landslide_susceptibility']
    else:
        val = 0 
    return val<jupyter_output><empty_output><jupyter_text>Only for plotting (without having to rerun analysis)<jupyter_code>sector_exp = pandas.read_csv(econ_haz_percent_exp_path)
merged_intersections = pandas.read_csv(econ_haz_cap_path)
econ_merged = pandas.read_csv(econ_merged_path)
all_intersections = pandas.read_csv(intersections_econ_path)
roads_exp_per = pandas.read_csv(roads_haz_percent_exp_path)
econ_merged = geopandas.GeoDataFrame(econ_merged)
econ_merged
len(econ_merged)
econ_merged.to_file(os.path.join(
    data_folder, 'data', 'infrastructure','econ_merged.shp'))<jupyter_output><empty_output><jupyter_text>Merge all economic assets, add administrative area to each sector<jupyter_code>def intersect_sector_with_admin(sectors, administrative, sector_id):
    """Intersect sectors with administrative
    """    
    intersections = geopandas.sjoin(sectors, administrative, how="inner", op='intersects')
    intersections = intersections[[
        'id_left',
        'capacity',
        'unit',
        'admin_name',
        'geometry'
    ]]
    intersections = intersections.rename(columns={
        'id_left': 'sector_id'
    })

    return intersections
econ_merged = []
for sector_id in sector_ids:
    sector_df = load_sector(data_folder, sector_id)
    sector_intersection = intersect_sector_with_admin(sector_df, administrative, sector_id)
    econ_merged.append(sector_intersection)

econ_merged = pandas.concat(econ_merged, axis=0)
econ_merged = econ_merged.fillna(0)
econ_merged.to_csv(econ_merged_path, index=False)

econ_merged.head()<jupyter_output><empty_output><jupyter_text>Intersections merged economics with each hazard as separate column<jupyter_code>def intersect_hazard_with_sector(hazards, sectors, hazard_id, sector_id):
    """Intersect sectors with hazards
    """
    ### INTERSECTIONS WITH HAZARDS        
    intersections_sector_haz = geopandas.sjoin(sectors, hazards, how="inner", op='intersects')

    intersections_sector_haz[hazard_id] = 1

    intersections_sector_haz['sector'] = sector_id

    intersections_sector_haz = intersections_sector_haz[[
        'sector',
        'id_left', 
        hazard_id
    ]]
    
    intersections_sector_haz = intersections_sector_haz.rename(columns={
        'id_left': 'sector_id'
        
    })
    intersections_sector_haz = intersections_sector_haz.drop_duplicates(subset='sector_id', keep='first')
    intersections_sector_haz = intersections_sector_haz.fillna(0)
    return intersections_sector_haz
all_intersections = []
for sector_id in sector_ids:
    sector_df = load_sector(data_folder, sector_id)
    
    for hazard_id in hazard_ids:
        hazard_df = load_hazard(data_folder, hazard_id)
        sector_intersection = intersect_hazard_with_sector(hazard_df, sector_df, hazard_id, sector_id)
        all_intersections.append(sector_intersection)
        
all_intersections = pandas.concat(all_intersections, axis=0, sort=False)
all_intersections = all_intersections.fillna(0)

all_intersections.to_csv(intersections_econ_path)
all_intersections.head()<jupyter_output><empty_output><jupyter_text>Point intersections for mapping<jupyter_code>sector_ids = ['wastewater_point', 'waste_point','electricity_point', 'freight', 'airports', 'port', 'water_point']
econ_merged = []
for sector_id in sector_ids:
    sector_df = load_sector(data_folder, sector_id)
    econ_merged.append(sector_df)

econ_merged = pandas.concat(econ_merged, axis=0)
econ_merged = econ_merged.fillna(0)
econ_merged = geopandas.GeoDataFrame(econ_merged)
econ_merged.to_file(os.path.join(data_folder,'results','econ_points_merge.shp'), index=False)

econ_merged.head()
for hazard_id in hazard_ids:
    hazard_df = load_hazard(data_folder, hazard_id)
    print(hazard_df)
    sector_intersection = geopandas.sjoin(econ_merged, hazard_df, how="inner", op='intersects')
    sector_intersection = geopandas.GeoDataFrame(all_intersections)
    sector_intersection.to_file(os.path.join(data_folder,'data','results','econ_merged_point_{}.shp'.format(hazard_id)))<jupyter_output>id  elevation                                           geometry
0        0          0  POLYGON ((513960.8484369325 1559718.924536677,...
1        1          1  POLYGON ((517676.475174851 1521324.000647028, ...
2        2          0  POLYGON ((513921.1223146291 1559682.98559937, ...
3        3          0  POLYGON ((520775.9107930948 1535575.979606579,...
4        4          1  POLYGON ((504712.2027580385 1520471.755477475,...
...    ...        ...                                                ...
3470  3470          0  POLYGON ((521305.2788649958 1532314.020539279,...
3471  3471          0  POLYGON ((517521.145317535 1556414.49932087, 5...
3472  3472          0  POLYGON ((499603.8596082596 1532216.20506565, ...
3473  3473          0  POLYGON ((520744.0103921712 1529583.372163905,...
3474  3474          0  POLYGON ((520835.7017037242 1539852.225847942,...

[3475 rows x 3 columns]<jupyter_text>Overview of all economic assets and whether they are exposed<jupyter_code>#import to not rerun all analyses
all_intersections = pandas.read_csv(intersections_econ_path)
all_intersections.head()
all_intersections = all_intersections.drop(columns= ['Unnamed: 0','Unnamed: 7'])
econ_merged.head()
merged_intersections = pandas.merge(
    econ_merged, all_intersections, how = 'outer', on=['sector','sector_id']) 
merged_intersections = merged_intersections.fillna(0)
merged_intersections.head()

merged_intersections.to_csv(econ_haz_count_path, index=False)<jupyter_output><empty_output><jupyter_text>Add capacity to each hazard<jupyter_code>merged_intersections
hazard_ids =['1m_sea-level', '4m_storm-surge', 'flashflooding', 'landslide_susceptibility']
for hazard_id in hazard_ids:
    merged_intersections[hazard_id] = merged_intersections[hazard_id]* merged_intersections['capacity']

merged_intersections.to_csv(econ_haz_cap_path, index=False)
merged_intersections.head()<jupyter_output><empty_output><jupyter_text>Add multiple hazard calculation<jupyter_code>merged_intersections = pandas.read_csv(econ_haz_cap_path)
list =['1m_sea-level', '4m_storm-surge', 'flashflooding', 'landslide_susceptibility', 'capacity']
#ensure all values are numeric
merged_intersections[list] = merged_intersections[list].apply(pandas.to_numeric)
# Merge information for each asset - grouping by 
# -> If sector_id is the same, add the columns for each of the hazards
merged_intersections = merged_intersections.groupby([
    'sector', 'sector_id', 'admin_name', 'unit', 'geometry'
    ])[
    list
    ].sum()
merged_intersections
#create a new hazard with combined hazard 
merged_intersections['storm_flash'] = merged_intersections.apply(two_multi_hazards, axis=1)
merged_intersections['storm_flash_landslide']= merged_intersections.apply(three_multi_hazards, axis=1)
merged_intersections
merged_intersections.to_csv(econ_multi_haz_asset_path)<jupyter_output><empty_output><jupyter_text>Calculate absolute exposed capacity per sector<jupyter_code>list =['1m_sea-level', '4m_storm-surge', 'flashflooding', 'landslide_susceptibility','storm_flash','storm_flash_landslide', 'capacity']
sector_exp = merged_intersections.groupby(['sector', 'unit'])[list].sum().reset_index()
sector_exp
sector_exp.to_csv(econ_haz_abs_path)<jupyter_output><empty_output><jupyter_text>Calculate % of exposed capacity per sector<jupyter_code>sector_exp_per = sector_exp.copy()
hazard_ids =['1m_sea-level', '4m_storm-surge', 'flashflooding', 'landslide_susceptibility', 'storm_flash', 'storm_flash_landslide']
for hazard_id in hazard_ids:
    sector_exp_per[hazard_id] = round((sector_exp_per[hazard_id] / sector_exp_per['capacity']*100))
    
sector_exp_per = sector_exp_per.fillna(0)



sector_exp_per.index = sector_exp_per.sector
#sector_exp_per = sector_exp_per.drop('sector', axis=1)

sector_exp_per.to_csv(econ_haz_percent_exp_path)
sector_exp_per.head()<jupyter_output><empty_output><jupyter_text>ADD ROADS<jupyter_code>#add roads as sector 
sector_exp_per = pandas.read_csv(econ_haz_percent_exp_path, index_col = 'sector')
roads_exp_per = pandas.read_csv(roads_haz_percent_exp_path, index_col = 'sector')
sector_exp_per = pandas.concat([sector_exp_per,roads_exp_per], axis=0,  sort=True)


roads_exp_per
sector_exp_per
sector_exp_per= sector_exp_per.replace('cargo','freight')
sector_exp_per 
sector_exp_per = sector_exp_per.drop(columns='sector.1')
sector_exp_per = sector_exp_per.fillna(0)
sector_exp_per.to_csv(econ_haz_percent_exp_path)<jupyter_output><empty_output><jupyter_text>Add SDG part<jupyter_code>SDG_direct_path = os.path.join(
        data_folder, 'data', 'SDG', 'SDG_influence_eco_soc_ecosystem_direct.csv')

SDG_direct = pandas.read_csv(SDG_direct_path, index_col ='Target')
SDG_direct.head()
SDG_direct = SDG_direct.sum()
SDG_direct_sec = pandas.DataFrame(SDG_direct)
SDG_direct_sec = SDG_direct_sec.rename(columns ={
    0:'SDG_direct'
})
SDG_direct_sec.index.name='sector'
SDG_direct_sec
SDG_exp = sector_exp_per.merge(SDG_direct_sec, on ='sector' )
SDG_exp
for hazard_id in hazard_ids:
    SDG_exp[hazard_id] = SDG_exp[hazard_id]* SDG_exp['SDG_direct']

SDG_exp.to_csv(SDG_cap_path, index=True)
SDG_exp<jupyter_output><empty_output><jupyter_text>PLOT ECONOMIC ASSETS<jupyter_code>#sector_exp = pandas.read_csv(econ_haz_percent_exp_path)
plot_econ_path = os.path.join(
        data_folder, 'figures', 'fig_economic_inf_stagged_hazards.svg')

plot_econ_SDG_path = os.path.join(
        data_folder, 'figures', 'fig_economic_SDG_hazards.svg')

SDG_cap_path = os.path.join(
    data_folder, 'results','econ_hazards_SDG_exp.csv')

SDG_direct_path = os.path.join(
        data_folder, 'data', 'SDG', 'SDG_influence_eco_soc_ecosystem_direct.csv')

SDG_direct = pandas.read_csv(SDG_direct_path)

SDG_exp = pandas.read_csv(SDG_cap_path)

sector_exp = pandas.read_csv(econ_haz_percent_exp_path)
#FOR SDG color=['red','darkred','yellow','lightcoral','darkgoldenrod','lightblue','blue']
def plot_to_ax(ax, title, hazard):
    ax = hazard.plot(ax=ax,kind='bar',use_index=True, legend = False,color=['red','green','yellow','orange','grey','lightblue','blue', 'black'], alpha=0.4)#x=sector_exp['sector']
    #
    ax.set_ylim(0,100)
    ax.set_title(title,fontweight = 'bold')
    ax.set_ylabel('Percentage of total capacity exposed')
    #ax.set_ylim(0, 100)
    ax.tick_params(axis='x', rotation = 90)
    ax.set_xlabel('economic infrastructures')
fig, axes = plt.subplots(nrows=2, ncols=2,sharex=True,sharey=True,figsize=(10,10), tight_layout = True)

plot_to_ax(axes[0,0], 'Sea-level rise (1 metre)', sector_exp['1m_sea-level'])
plot_to_ax(axes[0,1], 'Storm surge (4 metre)', sector_exp['4m_storm-surge'])
plot_to_ax(axes[1,0], 'Flashflooding (low, mod, high & very high susceptibility)', sector_exp['flashflooding'])
plot_to_ax(axes[1,1], 'Landslides (high susceptibility)', sector_exp['landslide_susceptibility'])
plt.savefig(plot_econ_path)
plt.show()   
plt.show()
fig, ax = plt.subplots(nrows=2, ncols=2,sharex=True,sharey=True,figsize=(10,10), tight_layout = True)

plot_to_ax(ax[0,0], 'Sea-level rise (1 metre)', sector_exp['1m_sea-level'])
plot_to_ax(ax[0,1], 'Storm surge (4 metre)', sector_exp['4m_storm-surge'])
plot_to_ax(ax[1,0], 'Flashflooding (low, mod, high & very high susceptibility)', sector_exp['flashflooding'])
plot_to_ax(ax[1,1], 'Landslides (high susceptibility)', sector_exp['landslide_susceptibility'])

ax2 = ax.twinx()
plot_to_ax2(ax2[0,0], 'Sea-level rise (1 metre)', SDG_exp['1m_sea-level'])
plot_to_ax2(ax2[0,1], 'Storm surge (4 metre)', SDG_exp['4m_storm-surge'])
plot_to_ax2(ax2[1,0], 'Flashflooding (low, mod, high & very high susceptibility)', SDG_exp['flashflooding'])
plot_to_ax2(ax2[1,1], 'Landslides (high susceptibility)', SDG_exp['landslide_susceptibility'])

plt.show()
plt.savefig(plot_econ_SDG_path)<jupyter_output><empty_output><jupyter_text>BY ADMIN CALCULATION  MERGE ADMIN WITH ALL ASSETS, CUMULATIVE RISK<jupyter_code>admin_econ_cum_exp_path = os.path.join(
    data_folder, 'results', 'admin_econ_cum_geo.csv')

admin_econ_cum_exp_disag_path = os.path.join(
    data_folder, 'results', 'admin_econ_cum_disag.csv')

admin_econ_exp_path = os.path.join(
    data_folder, 'results','admin_econ_sector_exp.csv')

merged_intersections = pandas.read_csv(econ_multi_haz_asset_path)
admin_econ_mean_exp_path= os.path.join(
    data_folder, 'results', 'admin_econ_mean.csv')

hazard_ids = ['1m_sea-level','4m_storm-surge','flashflooding','landslide_susceptibility', 'storm_flash', 'storm_flash_landslide']
merged_intersections
merged_intersections= merged_intersections.replace('cargo','freight')
merged_intersections
econ_admin = merged_intersections.groupby(['admin_name','sector', 'unit']).sum().reset_index()
econ_admin<jupyter_output><empty_output><jupyter_text>Exporting overview of percentage of sector exposed with admin info<jupyter_code>econ_admin_stagged = econ_admin.copy()
for hazard_id in hazard_ids:
    econ_admin_stagged[hazard_id] = round((econ_admin_stagged[hazard_id] / econ_admin_stagged.capacity)*100)
econ_admin_stagged = econ_admin_stagged.drop(columns={'unit', 'sector_id', 'capacity'})
econ_admin_stagged = econ_admin_stagged.set_index('admin_name')
econ_admin_stagged.to_csv(admin_econ_exp_path)
econ_admin_stagged<jupyter_output><empty_output><jupyter_text>Computing sector exposure for all different administrative areas<jupyter_code>def merge_admin(sector_id, administrative):
    sector = econ_admin[econ_admin['sector'] == sector_id]
    admin_econ_sector =  pandas.merge(administrative, sector, on ='admin_name', how='left')
    admin_econ_sector['sector'] = sector_id
    
    return admin_econ_sector<jupyter_output><empty_output><jupyter_text>Export each admin and sector<jupyter_code>for sector_id in sector_ids:
    admin_econ = merge_admin(sector_id, administrative)
    admin_econ = geopandas.GeoDataFrame(admin_econ)
    admin_econ.crs = {'init': 'epsg:2006'}
    admin_econ = admin_econ.reset_index()
    for hazard_id in hazard_ids:
        admin_econ[hazard_id] = round((admin_econ[hazard_id] / admin_econ['capacity'])*100)
    admin_econ = admin_econ.rename(columns={
        '1m_sea-level': 'sea-lev_1m',
        '4m_storm-surge': 'storm_s_4m',
        'flashflooding': 'flashflood', 
        'landslide_susceptibility':'landslide'
    })
    admin_econ = admin_econ.drop(columns=[ 'sector_id', 'index', 'id'])
    admin_econ = admin_econ.fillna(0)
    print(admin_econ)
    admin_econ_path = os.path.join(
        data_folder, 'results','admin_{}.shp'.format(sector_id))
    admin_econ.to_file(admin_econ_path)<jupyter_output>admin_name                                           geometry  \
0  Anse-la-Raye  POLYGON ((504489.1002518567 1543070.570047858,...   
1      Canaries  POLYGON ((504962.9645831303 1534981.799464189,...   
2      Castries  POLYGON ((504489.1002518567 1543070.570047858,...   
3      Choiseul  POLYGON ((507850.020189992 1529299.173142115, ...   
4       Dennery  POLYGON ((521289.1197935755 1543615.433000766,...   
5    Gros Islet  POLYGON ((514096.4336815997 1543239.529323319,...   
6       Laborie  POLYGON ((510296.0161504431 1518209.724388601,...   
7        Micoud  (POLYGON ((519386.4274498464 1533181.09109453,...   
8     Soufriere  POLYGON ((508887.6301650093 1530098.062839435,...   
9    Vieux Fort  (POLYGON ((507850.020189992 1529299.173142115,...   

       sector     unit  sea-lev_1m  storm_s_4m  flashflood  landslide  \
0  wastewater        0         0.0         0.0         0.0        0.0   
1  wastewater        0         0.0         0.0         0.0        0.0   
2  wastewa[...]<jupyter_text>[old_v072019]Combine administrative areas for each asset<jupyter_code>admin_econ_sector = []
for sector_id in sector_ids:
    admin_econ = merge_admin(sector_id, administrative)
    admin_econ_sector.append(admin_econ)
admin_econ_sector = pandas.concat(admin_econ_sector, axis=0)
admin_econ_sector_ind = admin_econ_sector.set_index('admin_name')
admin_econ_sector_ind
for hazard_id in hazard_ids:
    admin_econ_sector[hazard_id] = round((admin_econ_sector[hazard_id] / admin_econ_sector.capacity)*100)
    #sector_exp = rename_hazard_ids(hazard_id)
#sector_exp.index = sector_exp.sector
admin_econ_sector = admin_econ_sector.drop(['geometry','unit','sector_id','capacity','id'], axis=1)
admin_econ_sector_ind = admin_econ_sector.set_index('admin_name')
admin_econ_sector_ind
admin_econ_sector_ind.to_csv(admin_econ_cum_exp_disag_path)
admin_econ_cum_exp = admin_econ_sector.groupby(['admin_name']).sum().reset_index()

admin_econ_cum_exp
admin_econ_cum_exp.to_csv(admin_econ_cum_exp_path)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11809581934696368, 'fraction_numerical': 0.04557013528633913, 'mean_word_length': 5.6541455160744505, 'pattern_counts': {'":': 0, '<': 79, '<?xml version=': 0, '>': 85, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '359909', 'license': 'MIT', 'n_tokens_mistral': 8084, 'n_tokens_neox': 7199, 'n_words': 1077, 'path': 'src/intersection/hazards_econ.ipynb', 'repo_name': 'nismod/stlucia-adaptation', 'size': '19662'}