from sunau import Au_write import pint # module to work with units; https://pypi.org/project/Pint/0.10/ import pandas as pd ureg = pint.UnitRegistry() import openpyxl import os # ********************************************************************************************** # Define custom units in pint ureg.define("dollar = [currency]") ureg.define("vehicle = [dimensionless]") ureg.define("person = [dimensionless]") ureg.define("accident = [dimensionless]") # ******************************************************************************************************* # Fill this section: # System id: system = "D" environment = "C5" # Possible values: C2, C3, C5 n_layers_equal_2 = "yes" # Number of coating layers in the system, 2 or 3. # 2 means no intermediate, just primer and top coat. surf_treat_SP6_yn = "no" surf_treat_SP10_yn = "yes" # Thickness of each layer: see Table 11 if system == "A": primer_thick_mils = 3 * ureg.mil # mils interm_thick_mils = 3 * ureg.mil # mils if n_layers_equal_2 == "yes": interm_thick_mils = 0 * ureg.mil # mils top_coat_thick_mils = 3 * ureg.mil # mils if system == "B": primer_thick_mils = 3 * ureg.mil # mils interm_thick_mils = 3 * ureg.mil # mils if n_layers_equal_2 == "yes": interm_thick_mils = 0 * ureg.mil # mils top_coat_thick_mils = 3 * ureg.mil # mils if system == "C": primer_thick_mils = 2 * ureg.mil # mils interm_thick_mils = 2 * ureg.mil # mils if n_layers_equal_2 == "yes": interm_thick_mils = 0 * ureg.mil # mils top_coat_thick_mils = 2 * ureg.mil # mils if system == "D": primer_thick_mils = 3 * ureg.mil # mils interm_thick_mils = 0 * ureg.mil # mils if n_layers_equal_2 == "yes": interm_thick_mils = 0 * ureg.mil # mils top_coat_thick_mils = 4 * ureg.mil # mils # Surface area to coat A_ft2 = 1e5 * ureg.ft**2 # surface area, ft2, Quebec Bridge 495,000 m2; 5328130 ft2 https://www.canada.ca/en/housing-infrastructure-communities/news/2021/04/did-you-know.html # Eq.2 # Coating cost during the 75-year service life of the bridge: # Attributing values to n1, n2, n3 for each case scenario if environment == "C2": if system == "A" and surf_treat_SP6_yn == "yes": n1 = 2 # number of times spot painting n2 = 1 # number of times over coating n3 = 1 # number of times removal and replacement if system == "A" and surf_treat_SP6_yn == "no": n1 = 1 # number of times spot painting n2 = 1 # number of times over coating n3 = 1 # number of times removal and replacement if system == "B" and surf_treat_SP6_yn == "yes": n1 = 1 # number of times spot painting n2 = 1 # number of times over coating n3 = 1 # number of times removal and replacement if system == "B" and surf_treat_SP6_yn == "no": n1 = 1 # number of times spot painting n2 = 1 # number of times over coating n3 = 1 # number of times removal and replacement if system == "C" and surf_treat_SP6_yn == "yes": n1 = 2 # number of times spot painting n2 = 2 # number of times over coating n3 = 2 # number of times removal and replacement if system == "C" and surf_treat_SP6_yn == "no": n1 = 2 # number of times spot painting n2 = 2 # number of times over coating n3 = 2 # number of times removal and replacement if system == "D" and surf_treat_SP6_yn == "yes": n1 = 1 # number of times spot painting n2 = 1 # number of times over coating n3 = 1 # number of times removal and replacement if system == "D" and surf_treat_SP6_yn == "no": n1 = 1 # number of times spot painting n2 = 1 # number of times over coating n3 = 1 # number of times removal and replacement if environment == "C3": if system == "A" and surf_treat_SP6_yn == "yes": n1 = 2 # number of times spot painting n2 = 2 # number of times over coating n3 = 2 # number of times removal and replacement if system == "A" and surf_treat_SP6_yn == "no": n1 = 2 # number of times spot painting n2 = 1 # number of times over coating n3 = 1 # number of times removal and replacement if system == "B" and surf_treat_SP6_yn == "yes": n1 = 2 # number of times spot painting n2 = 1 # number of times over coating n3 = 1 # number of times removal and replacement if system == "B" and surf_treat_SP6_yn == "no": n1 = 2 # number of times spot painting n2 = 1 # number of times over coating n3 = 1 # number of times removal and replacement if system == "C" and surf_treat_SP6_yn == "yes": n1 = 3 # number of times spot painting n2 = 3 # number of times over coating n3 = 3 # number of times removal and replacement if system == "C" and surf_treat_SP6_yn == "no": n1 = 3 # number of times spot painting n2 = 3 # number of times over coating n3 = 2 # number of times removal and replacement if system == "D" and surf_treat_SP6_yn == "yes": n1 = 2 # number of times spot painting n2 = 2 # number of times over coating n3 = 2 # number of times removal and replacement if system == "D" and surf_treat_SP6_yn == "no": n1 = 2 # number of times spot painting n2 = 2 # number of times over coating n3 = 2 # number of times removal and replacement if environment == "C5": if system == "A" and surf_treat_SP6_yn == "yes": n1 = 3 # number of times spot painting n2 = 3 # number of times over coating n3 = 2 # number of times removal and replacement if system == "A" and surf_treat_SP6_yn == "no": n1 = 3 # number of times spot painting n2 = 2 # number of times over coating n3 = 2 # number of times removal and replacement if system == "B" and surf_treat_SP6_yn == "yes": n1 = 3 # number of times spot painting n2 = 2 # number of times over coating n3 = 2 # number of times removal and replacement if system == "B" and surf_treat_SP6_yn == "no": n1 = 3 # number of times spot painting n2 = 2 # number of times over coating n3 = 2 # number of times removal and replacement if system == "C" and surf_treat_SP6_yn == "yes": n1 = 5 # number of times spot painting n2 = 5 # number of times over coating n3 = 4 # number of times removal and replacement if system == "C" and surf_treat_SP6_yn == "no": n1 = 4 # number of times spot painting n2 = 4 # number of times over coating n3 = 3 # number of times removal and replacement if system == "D" and surf_treat_SP6_yn == "yes": n1 = 3 # number of times spot painting n2 = 3 # number of times over coating n3 = 2 # number of times removal and replacement if system == "D" and surf_treat_SP6_yn == "no": n1 = 3 # number of times spot painting n2 = 3 # number of times over coating n3 = 2 # number of times removal and replacement # ******************************************************************************************************* # Data for the 4 systems considered in this paper if system == "A": # Organic zinc (OZ) /Epoxy/Polyurethane (PU) three-coat system primer_cost = 0.140 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil primer_spray_cost = 0.200 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil primer_voc_kg_per_m3 = 240 * ureg.kg/ureg.m**3 # 1 g/L = 1 kg/m3 primer_field_app_cost = 0.79 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 perc_solid_primer = 65 * ureg.percent interm_cost = 0.065 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil interm_spray_cost = 0.093 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil if n_layers_equal_2 == "yes": interm_cost = 0 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) interm_spray_cost = 0 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) interm_voc_kg_per_m3 = 210 * ureg.kg / ureg.m ** 3 if n_layers_equal_2 == "yes": interm_voc_kg_per_m3 = 0 * ureg.kg / ureg.m ** 3 # 1 g/L = 1 kg/m3 interm_field_app_cost = 0.67 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 if n_layers_equal_2 == "yes": interm_field_app_cost = 0 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 perc_solid_interm = 75 * ureg.percent top_coat_cost = 0.106 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil top_coat_spray_cost = 0.152 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil top_coat_voc_kg_per_m3 = 270 * ureg.kg/ureg.m**3 top_coat_field_app_cost = 0.72 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 perc_solid_top_coat = 65 * ureg.percent if system == "B": # Moisture-cured (MC) zinc rich polyurethane/MC urethane/MC urethane three-coat system primer_cost = 0.133 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil primer_spray_cost = 0.189 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil primer_voc_kg_per_m3 = 335 * ureg.kg/ureg.m**3 # 1 g/L = 1 kg/m3 primer_field_app_cost = 0.79 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 perc_solid_primer = 65 * ureg.percent interm_cost = 0.100 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil interm_spray_cost = 0.143 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil if n_layers_equal_2 == "yes": interm_cost = 0 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) interm_spray_cost = 0 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) interm_voc_kg_per_m3 = 335 * ureg.kg / ureg.m ** 3 if n_layers_equal_2 == "yes": interm_voc_kg_per_m3 = 0 * ureg.kg / ureg.m ** 3 # 1 g/L = 1 kg/m3 interm_field_app_cost = 0.72 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 if n_layers_equal_2 == "yes": interm_field_app_cost = 0 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 perc_solid_interm = 65 * ureg.percent top_coat_cost = 0.100 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil top_coat_spray_cost = 0.143 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil top_coat_voc_kg_per_m3 = 340 * ureg.kg/ureg.m**3 top_coat_field_app_cost = 0.72 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 perc_solid_top_coat = 65 * ureg.percent if system == "C": # Waterborne (WB) acrylic/WB acrylic/WB acrylic three-coat system primer_cost = 0.099 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil primer_spray_cost = 0.141 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil primer_voc_kg_per_m3 = 75 * ureg.kg/ureg.m**3 # 1 g/L = 1 kg/m3 primer_field_app_cost = 0.55 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 perc_solid_primer = 45 * ureg.percent interm_cost = 0.105 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil interm_spray_cost = 0.150 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil if n_layers_equal_2 == "yes": interm_cost = 0 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) interm_spray_cost = 0 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) interm_voc_kg_per_m3 = 75 * ureg.kg / ureg.m ** 3 if n_layers_equal_2 == "yes": interm_voc_kg_per_m3 = 0 * ureg.kg / ureg.m ** 3 # 1 g/L = 1 kg/m3 interm_field_app_cost = 0.55 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 if n_layers_equal_2 == "yes": interm_field_app_cost = 0 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 perc_solid_interm = 45 * ureg.percent top_coat_cost = 0.105 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil top_coat_spray_cost = 0.150 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil top_coat_voc_kg_per_m3 = 75 * ureg.kg/ureg.m**3 top_coat_field_app_cost = 0.55 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 perc_solid_top_coat = 45 * ureg.percent if system == "D": # Organic zinc (OZ)/Polyurethane two-coat system primer_cost = 0.140 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil primer_spray_cost = 0.200 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil primer_voc_kg_per_m3 = 240 * ureg.kg/ureg.m**3 # 1 g/L = 1 kg/m3 primer_field_app_cost = 0.79 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 perc_solid_primer = 65 * ureg.percent if n_layers_equal_2 == "yes": interm_cost = 0 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) interm_spray_cost = 0 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) if n_layers_equal_2 == "yes": interm_voc_kg_per_m3 = 0 * ureg.kg / ureg.m ** 3 # 1 g/L = 1 kg/m3 if n_layers_equal_2 == "yes": interm_field_app_cost = 0 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 top_coat_cost = 0.106 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil top_coat_spray_cost = 0.152 * (ureg.dollar / (ureg.mil * ureg.ft ** 2)) # $/ft2/mil top_coat_voc_kg_per_m3 = 270 * ureg.kg/ureg.m**3 top_coat_field_app_cost = 0.72 * (ureg.dollar / (ureg.ft ** 2)) # $/ft2 perc_solid_top_coat = 65 * ureg.percent if surf_treat_SP6_yn == "yes": surf_treat_cost = 1.80 * (ureg.dollar/(ureg.ft**2)) # $/ft2m if surf_treat_SP10_yn == "yes": surf_treat_cost = 2.47 * (ureg.dollar/(ureg.ft**2)) # $/ft2m surf_treat_cost_total = surf_treat_cost * A_ft2 original_cost = A_ft2 * (surf_treat_cost + primer_thick_mils * (primer_spray_cost) + interm_thick_mils * (interm_spray_cost) + top_coat_thick_mils * (top_coat_spray_cost) + primer_field_app_cost + interm_field_app_cost + top_coat_field_app_cost) print(original_cost) # Original costs need to be separated in 4 categories, as per Fig. 4 # 1. Material # 2. Surface preparation # 3. Coating application # 4. Access & HSE # 1. Material original_cost_material = A_ft2 * ((primer_thick_mils * primer_spray_cost) + (interm_thick_mils * interm_spray_cost) + (top_coat_thick_mils * top_coat_spray_cost)) # 2. Surface preparation original_cost_surf_treat = A_ft2 * surf_treat_cost # 3. Coating application original_cost_application = A_ft2 * (primer_field_app_cost + interm_field_app_cost + top_coat_field_app_cost) # 4. for HSE calculations, will be a constant # Refer to original_cost equation for variables id # Needed for access_hse_cost costs surf_treat_cost_SP10 = 2.47 * (ureg.dollar/(ureg.ft**2)) # $/ft2 --> reference for hse original_cost_SP10 = A_ft2 * (surf_treat_cost_SP10 + 3 * ureg.mil * (0.200 * (ureg.dollar / (ureg.mil * ureg.ft ** 2))) + 3 * ureg.mil * (0.093 * (ureg.dollar / (ureg.mil * ureg.ft ** 2))) + 3 * ureg.mil * (0.152 * (ureg.dollar / (ureg.mil * ureg.ft ** 2))) + 0.79 * (ureg.dollar / (ureg.ft ** 2)) + 0.67 * (ureg.dollar / (ureg.ft ** 2)) + 0.72 * (ureg.dollar / (ureg.ft ** 2))) print(original_cost_SP10) # For the whole 75-yr lifespan: # 1. Material overall_75yr_cost_material = (1 + 0.4*n1 + 0.7*n2 + 1.35*n3)*original_cost_material # includes original cost # 2. Surface preparation overall_75yr_surf_treat = (1 + 0.4*n1 + 0.7*n2 + 1.35*n3)*original_cost_surf_treat # includes original cost # 3. Coating application overall_75yr_cost_application = (1 + 0.4*n1 + 0.7*n2 + 1.35*n3)*original_cost_application # includes original cost # 4. Access & HSE overall_75yr_access_hse_cost = 3*(1 + 0.4*n1 + 0.7*n2 + 1.35*n3)*original_cost_SP10 # original_cost_SP10, baseline for access_hse_cost for ALL systems # Summary of costs for Fig.4 overall_75yr_cost_material overall_75yr_surf_treat overall_75yr_cost_application overall_75yr_access_hse_cost total_cost = (overall_75yr_cost_material + overall_75yr_surf_treat + overall_75yr_cost_application + overall_75yr_access_hse_cost) print(f"1. Material cost is {overall_75yr_cost_material:.2f}", (overall_75yr_cost_material/total_cost)) print(f"2. Surface preparation cost is {overall_75yr_surf_treat:.2f}", overall_75yr_surf_treat/total_cost) print(f"3. Coating application cost is {overall_75yr_cost_application:.2f}", overall_75yr_cost_application/total_cost) print(f"4. HSE cost is {overall_75yr_access_hse_cost:.2f}", overall_75yr_access_hse_cost/total_cost) coating_cost = (1 + 0.4*n1 + 0.7*n2 + 1.35*n3)*original_cost # includes original cost print(coating_cost) # Eq.3 access_hse_cost = 3*(1 + 0.4*n1 + 0.7*n2 + 1.35*n3)*original_cost_SP10 # original_cost_SP10, baseline for access_hse_cost for ALL systems print(access_hse_cost) # Eq.4 project_cost = coating_cost + access_hse_cost print(project_cost) field_app_cost_total = A_ft2 * (primer_field_app_cost + interm_field_app_cost + top_coat_field_app_cost) # ******************************************************************************************************* # VOC *************************************************************************************************** # Eq.5 A_m2 = A_ft2 / (10.764 * ureg.ft**2 / ureg.m**2) print(A_m2) primer_thick_m = primer_thick_mils / (39370 * ureg.mil/ureg.m) interm_thick_m = interm_thick_mils / (39370 * ureg.mil/ureg.m) top_coat_thick_m = top_coat_thick_mils / (39370 * ureg.mil/ureg.m) if n_layers_equal_2 == "yes": VOC_interm = 0 VOC_primer = (primer_thick_m/(perc_solid_primer/(100 * ureg.percent))) * A_m2 * primer_voc_kg_per_m3 # kg if n_layers_equal_2 == "no": VOC_interm = (interm_thick_m/(perc_solid_interm/(100 * ureg.percent))) * A_m2 * interm_voc_kg_per_m3 VOC_top_coat = (top_coat_thick_m/(perc_solid_top_coat/(100 * ureg.percent))) * A_m2 * top_coat_voc_kg_per_m3 # Eq.6 VOC_coat_sys_coating = VOC_primer + VOC_interm + VOC_top_coat print(VOC_coat_sys_coating) # Eq.7 # If 2 n_layers_equal_2 == "yes", VOC_interm is set to 0 above in the code. VOC_coat_sys_total = (1+n3)*VOC_coat_sys_coating + n2*(VOC_interm + VOC_top_coat) print(VOC_coat_sys_total) # ******************************************************************************************************* # Social cost ******************************************************************************************* # Extra vehicle operation cost ************************************************************************** # Variables for Eq.8 to Eq.15 ADT = 30000 * ureg.vehicle/ureg.day # Average daily traffic (ADT) Cr_car = 0.45 * ureg.dollar/ureg.km/ureg.vehicle # Car operation cost (Cr, car) Cr_truck = 1.50 * ureg.dollar/ureg.km/ureg.vehicle # Truck operation cost (Cr, car) D_l = 5 * ureg.km # Detour extra mileage (D_l) TT_p = 5 * ureg.percent # Percentage of trucks in ADT (TTP) FC_car = 13 * ureg.L/(100*ureg.km)/ureg.vehicle # Fuel consumption for cars FC_truck = 40 * ureg.L/(100*ureg.km)/ureg.vehicle # Fuel consumption for trucks C_gas = 0.90 * ureg.dollar/ureg.L C_diesel = 1.00 * ureg.dollar/ureg.L FCI = 25 * ureg.percent # Fuel consumption increase (FCI) P_t = 10 * ureg.percent # Percentage of the time for total closure d_r = 7 * ureg.day/ureg.person # Average recovery time for one person injured C_w = 19.30 * ureg.dollar/ureg.h # Average hourly income t_w = 7.5 * ureg.h/ureg.day # Daily working hours PPT_car = 1.50 * ureg.person/ureg.vehicle # Persons per trip for cars PPT_truck = 1.00 * ureg.person/ureg.vehicle # Persons per trip for trucks PPA = 1 * ureg.person/ureg.accident # Persons per accident speed_limit_reqular = 80 * ureg.km/ureg.h speed_limit_constr = 50 * ureg.km/ureg.h speed_limit_detour = 80 * ureg.km/ureg.h n_lane_regular = 4 n_lane_constr = 2 ANA = 0.02 * ureg.accident/ureg.vehicle/ureg.year # Average number of accidents r_injury = 0.28 # Rate of injuries RAI = 0.25 # Rate of accident increase r_damage = 0.71 # Rate of property damage in road accidents C_d = 4711 * ureg.dollar/ureg.accident # Average compensation value for property damage days_per_yr = 365 * ureg.day/ureg.year L = 4.57 * ureg.km # Maintenance affected distance dt_spot_painting = 30 * ureg.day # duration of maintenance activities - spot painting dt_over_coating = 90 * ureg.day # duration of maintenance activities - over coating dt_remove_replace = 180 * ureg.day # duration of maintenance activities - remove & replace hundred_pc = 100 * ureg.percent # Eq.8 extra_veh_oper_cost = ( n1 * ((Cr_car*(1-TT_p/hundred_pc) + Cr_truck*(TT_p/hundred_pc)) * D_l * (P_t/hundred_pc) * ADT * dt_spot_painting) + n2 * ((Cr_car*(1-TT_p/hundred_pc) + Cr_truck*(TT_p/hundred_pc)) * D_l * (P_t/hundred_pc) * ADT * dt_over_coating) + n3 * ((Cr_car*(1-TT_p/hundred_pc) + Cr_truck*(TT_p/hundred_pc)) * D_l * (P_t/hundred_pc) * ADT * dt_remove_replace) ) print(extra_veh_oper_cost) # same than Misagh # Eq.9 # Original coating process does not need to be added. Therefore: extra_fuel_cost = ( n1 * (((1-TT_p/hundred_pc)*FC_car*C_gas) + ((TT_p/hundred_pc)*FC_truck*C_diesel)) * # spot painting L * ADT * dt_spot_painting * (1-P_t/hundred_pc) * (FCI/hundred_pc) ) + ( n2 * (((1-TT_p/hundred_pc)*FC_car*C_gas) + ((TT_p/hundred_pc)*FC_truck*C_diesel)) * # over coating L * ADT * dt_over_coating * (1-P_t/hundred_pc) * (FCI/hundred_pc) ) + ( n3 * (((1-TT_p/hundred_pc)*FC_car*C_gas) + ((TT_p/hundred_pc)*FC_truck*C_diesel)) * # remove & replace L * ADT * dt_remove_replace * (1-P_t/hundred_pc) * (FCI/hundred_pc) ) print(extra_fuel_cost) # Eq.10 extra_veh_cost = extra_veh_oper_cost + extra_fuel_cost print(extra_veh_cost) # Eq.12 t_loss_detour = ( ((L+D_l)/speed_limit_detour) - (L/speed_limit_reqular) ) * (1/ureg.person) # *** added, / person, to get h/person print(t_loss_detour) # Eq.13 t_loss_lane_red = ( (n_lane_regular/n_lane_constr)*(L/speed_limit_constr) - (L/speed_limit_reqular) ) * (1/ureg.person) print(t_loss_lane_red) # Eq.11 # Original coating process does not need to be added. Therefore: travel_delay_cost = ( n1 * (PPT_car*(1-TT_p/hundred_pc) + PPT_truck*(TT_p/hundred_pc)) * ADT * dt_spot_painting * C_w * ((t_loss_detour*P_t/hundred_pc) + t_loss_lane_red*(1-P_t/hundred_pc)) ) + ( n2 * (PPT_car*(1-TT_p/hundred_pc) + PPT_truck*(TT_p/hundred_pc)) * ADT * dt_over_coating * C_w * ((t_loss_detour*P_t/hundred_pc) + t_loss_lane_red*(1-P_t/hundred_pc)) ) + ( n3 * (PPT_car*(1-TT_p/hundred_pc) + PPT_truck*(TT_p/hundred_pc)) * ADT * dt_remove_replace * C_w * ((t_loss_detour*P_t/hundred_pc) + t_loss_lane_red*(1-P_t/hundred_pc)) ) print(travel_delay_cost) # Eq.14 income_loss_cost = ( n1 * PPA * (ADT * dt_spot_painting) * (ANA * dt_spot_painting / days_per_yr) * r_injury * RAI * d_r * C_w * t_w ) + ( n2 * PPA * (ADT * dt_over_coating) * (ANA * dt_over_coating / days_per_yr) * r_injury * RAI * d_r * C_w * t_w ) + ( n3 * PPA * (ADT * dt_remove_replace) * (ANA * dt_remove_replace / days_per_yr) * r_injury * RAI * d_r * C_w * t_w ) print(income_loss_cost) # Eq.15 property_damage_cost = ( n1 * ADT * dt_spot_painting * (ANA * dt_spot_painting / days_per_yr) * r_damage * RAI * C_d ) + ( n2 * ADT * dt_over_coating * (ANA * dt_over_coating / days_per_yr) * r_damage * RAI * C_d ) + ( n3 * ADT * dt_remove_replace * (ANA * dt_remove_replace / days_per_yr) * r_damage * RAI * C_d ) print(property_damage_cost) # Eq.16 social_cost = extra_veh_cost + travel_delay_cost + income_loss_cost + property_damage_cost print(social_cost) # Summary print(f"1. The VOC content is {VOC_coat_sys_total:.2f} kg") print(f"2. The project cost is {project_cost/1e6:.2f} M$ (coating_cost, access_hse_cost)") print(f" 2.1. The coating cost is {coating_cost/1e6:.2f} M$") print(f" 2.2. The HSE cost is {access_hse_cost/1e6:.2f} M$") print(f"3. The social cost is {social_cost/1e6:.2f} M$ (extra_veh_cost, travel_delay_cost, income_loss_cost, property_damage_cost)") print(f" 3.1. The extra vehicle cost is {extra_veh_cost/1e6:.2f} M$") print(f" 3.2. The travel delay cost is {travel_delay_cost/1e6:.2f} M$") print(f" 3.3. The income loss is {income_loss_cost/1e6:.2f} M$") print(f" 3.4. The property damage cost is {property_damage_cost/1e6:.2f} M$") # *********************************************************************************************************** # Graphs **************************************************************************************************** # Graph 1 Pie chart showing the split of costs, as follows: (for Fig. 4) # # 1. Material # # 2. Surface preparation # # 3. Coating application # # 4. Access & HSE import matplotlib.pyplot as plt # Define the directory where you want to save the image save_directory = "C:/Users/gallantd/Documents/Industriels/Coating Project CONSTR/Paper Sept. 2024/Figures/version reviewed Febr2025/" # Change this to your desired folder # Ensure the directory exists os.makedirs(save_directory, exist_ok=True) # Prepare labels dynamically by converting pint values to million dollars labels = [ f'Material Cost ({(overall_75yr_cost_material / 1e6).magnitude:.1f} M$)', f'Surface Preparation ({(overall_75yr_surf_treat / 1e6).magnitude:.2f} M$)', f'Coating Application ({(overall_75yr_cost_application / 1e6).magnitude:.2f} M$)', f'Access & HSE Cost ({(overall_75yr_access_hse_cost / 1e6).magnitude:.1f} M$)' ] # Extract magnitudes for the pie chart values = [ (overall_75yr_cost_material / 1e6).magnitude, (overall_75yr_surf_treat / 1e6).magnitude, (overall_75yr_cost_application / 1e6).magnitude, (overall_75yr_access_hse_cost / 1e6).magnitude ] # Create the pie chart plt.figure(figsize=(8, 6)) wedges, texts, autotexts = plt.pie( values, autopct='%1.1f%%', # Format for the percentage startangle=140, pctdistance=1.2, # Standard distance for percentage text labeldistance=1.4, # Distance of labels from the center wedgeprops=dict(linewidth=1, edgecolor="black") # Add slice borders ) # Set font size for labels for text in texts: text.set_fontsize(14) # Increase label font size # Adjust percentage text to avoid overlaps for i, autotext in enumerate(autotexts): autotext.set_color('black') # Set percentage text color autotext.set_fontsize(14) autotext.set_horizontalalignment('center') autotext.set_verticalalignment('center') # Adjust the vertical position of overlapping percentages bbox = autotext.get_window_extent() for j, other in enumerate(autotexts): if i != j: # Compare with other percentage texts bbox_other = other.get_window_extent() if bbox.overlaps(bbox_other): # Check for overlap autotext.set_y(autotext.get_position()[1] + 0.05 * (j - i)) # Slightly offset overlapping texts # Add a legend at the bottom of the chart with a darker border legend = plt.legend( wedges, labels, title="Cost Categories", loc="upper center", # Position the legend at the top of bbox bbox_to_anchor=(0.5, -0.02), # Centered below the pie chart ncol=2, # Arrange legend items in two columns fontsize=12, # Increase font size for legend title_fontsize=14, # Increase font size for the legend title frameon=True, # Add a frame around the legend edgecolor='black', # Darker border color fancybox=False # Remove rounded corners ) # Adjust the border linewidth legend.get_frame().set_linewidth(1.5) # Make the border thicker # Adjust layout to make space for the legend plt.tight_layout() # Define file path file_path = os.path.join(save_directory, "systemD_project_cost_c5_sp10.jpg") # Save the figure as a JPG file plt.savefig(file_path, format="jpg", dpi=300) # Display the plot plt.show() # Close the plot to free memory plt.close() # Graph 2 **************************************************************** import matplotlib.pyplot as plt import pandas as pd import numpy as np import matplotlib.ticker as ticker # Define the directory where you want to save the image save_directory = "C:/Users/gallantd/Documents/Industriels/Coating Project CONSTR/Paper Sept. 2024/Figures/version reviewed Febr2025/" # Change this to your desired folder # Ensure the directory exists os.makedirs(save_directory, exist_ok=True) # Read the data from the Excel file file_path = "C:/Users/gallantd/Documents/Industriels/Coating Project CONSTR/Paper Sept. 2024/250221_3_info_graphs_Fig1_2_3_v_reviewed.xlsx" # Select which tab to read: #df = pd.read_excel(file_path, sheet_name = 'fig1_c2_sp6' ) #df = pd.read_excel(file_path, sheet_name = 'fig2_c3_sp10' ) df = pd.read_excel(file_path, sheet_name = 'fig3_c5_sp10' ) # Replace '/' with '\n' in the System column for multiline labels df["System"] = df["System"].str.replace("/", "/\n") # # Define the x positions for the bars x = np.arange(len(df["System"])) # Create a range for the x-axis positions bar_width = 0.35 # Width of the bars # Create the figure and axes fig, ax1 = plt.subplots() # Plot the bar plots for Project Cost and Social Cost side by side bars_project = ax1.bar(x - bar_width / 2, df["Project_Cost"], bar_width, label="Project Cost", alpha=0.7) bars_social = ax1.bar(x + bar_width / 2, df["Social_Cost"], bar_width, label="Social Cost", alpha=0.7) # Add numerical values on top of the bars # Add numerical values on top of the bars with 1 decimal place for bar in bars_project: height = bar.get_height() ax1.text(bar.get_x() + bar.get_width() / 2, height, f"{height:.1f}", ha="center", va="bottom") for bar in bars_social: height = bar.get_height() ax1.text(bar.get_x() + bar.get_width() / 2, height, f"{height:.0f}", ha="center", va="bottom") # Configure the first y-axis (left) for costs with a logarithmic scale ax1.set_yscale("log") # Set the y-axis to logarithmic scale ax1.set_ylim(1, 1200) # Set the limits for the left y-axis ax1.set_xlabel("Systems", fontsize=14, weight='bold') # Make x-axis label bold ax1.set_ylabel("Project & Social Costs (M$)", color="black", fontsize=14, weight='bold') # Make left y-axis label bold ax1.set_xticks(x) ax1.set_xticklabels(df["System"], fontsize=12) # Update x-axis labels to include multiline labels ax1.tick_params(axis="y", labelcolor="black", labelsize=11) # Increase the font size for y-axis tick labels ax1.legend(loc="upper left", fontsize=12) # Increase the font size for the legend # Create a twin y-axis for VOC ax2 = ax1.twinx() ax2.scatter(x, df["VOC"], color="red", label="VOC", zorder=5, s=100) # Add VOC values as labels next to red dots for i, voc in enumerate(df["VOC"]): ax2.text(x[i] + 0.03, voc - 350, f"{voc:.0f}", color="red", fontsize=10, ha="left", va="bottom") ax2.set_ylim(0, 5000) # Set the limits for the right y-axis # Set major ticks every 1000 and minor ticks every 250 ax2.yaxis.set_major_locator(ticker.MultipleLocator(500)) # Major ticks: 0, 1000, 2000, ..., 5000 ax2.yaxis.set_minor_locator(ticker.MultipleLocator(100)) # Minor ticks: 250, 500, 750, ..., 4750 ax2.set_ylabel("VOC (kg)", color="red", fontsize=14, weight='bold') # Make right y-axis label bold ax2.tick_params(axis="y", labelcolor="red", labelsize=11) # Add gridlines for readability ax1.grid(axis="y", linestyle="--", alpha=0.7) # Show the plot plt.tight_layout() # # Select which file path to save in #file_path_to_save = os.path.join(save_directory, "fig1_c2_sp6__bar_plot.jpg") #file_path_to_save = os.path.join(save_directory, "fig2_c3_sp10__bar_plot.jpg") file_path_to_save = os.path.join(save_directory, "fig3_c5_sp10__bar_plot.jpg") # Save the figure as a JPG file plt.savefig(file_path_to_save, format="jpg", dpi=300) plt.show() # Close the plot to free memory plt.close() # ******************************************************************************************************** # Graph 3 ************************************************************************************************ # Timeline for recoating ********************************************************************************* # Read the 241207_timeline_recoating_seq.xlsx file # Figures 1-2-3, replace below the ws = value, and the .jpg name for saving, at the end of the code import pandas as pd import matplotlib.pyplot as plt from openpyxl import load_workbook # File path for the Excel file file_path = 'C:/Users/gallantd/Documents/Industriels/Coating Project CONSTR/Paper Sept. 2024/241207_timeline_recoating_seq.xlsx' # Load the workbook and sheet wb = load_workbook(file_path, data_only=True) # Select 1 of the 3 ws: ws = wb['c2_sp6'] ws = wb['c3_sp10'] ws = wb['c5_sp10'] # Read the data into a DataFrame (no header in the sheet) data = pd.DataFrame([[cell.value for cell in row] for row in ws.iter_rows()]) data.reset_index(drop=True, inplace=True) # Reset the index to ensure consistency # Extract text formatting (bold, italic, underlined) formats = [] for row in ws.iter_rows(min_row=1, max_row=ws.max_row, min_col=2, max_col=ws.max_column): row_formats = [] for cell in row: if cell.font: if cell.font.bold: row_formats.append("bold") # Representing green vertical lines elif cell.font.italic: row_formats.append("italic") # Representing blue vertical lines elif cell.font.underline: row_formats.append("underline") # Representing red vertical lines else: row_formats.append(None) else: row_formats.append(None) formats.append(row_formats) # Customize figure size (width x height in inches) figure_width = 14 # Adjust width figure_height = 8 # Adjust height # Create the figure and axis fig, ax = plt.subplots(figsize=(figure_width, figure_height)) ax.set_xlim(0, 75) # Set x-axis range from 0 to 75 ax.set_xticks(range(0, 76, 5)) # Add ticks every 5 units #ax.set_ylim(-1, len(data)) # One row per system # Adjust the y-axis limits to reduce space above the upper coating system ax.set_ylim(-1, len(data) - 0.3) # Reduce space at the top # Add pale gray dotted vertical lines every 10 years for x in range(0, 76, 10): # Loop through every 10 years from 0 to 75 ax.axvline(x=x, color='lightgray', linestyle='dotted', linewidth=2) # Preprocess y-axis labels to add line breaks after every '/' system_labels = [str(row.iloc[0]).replace("/", "/\n") for _, row in data.iterrows()] # Enable y-axis ticks and labels ax.tick_params(axis="y", which="both", left=True, labelleft=True, right=False, labelright=False) ax.set_yticks(range(len(data))) ax.set_yticklabels(system_labels, fontsize=10) # Draw arrows as backgrounds for each timeline arrow_width = 0.4 # Set arrow width for i in range(len(data)): ax.arrow( 0, i, 75, 0, color='darkgray', alpha=0.3, width=arrow_width, head_width=0.6, head_length=2.5, length_includes_head=True ) # Manually add unique handles for the legend legend_handles = { "bold": plt.Line2D([0], [0], color="green", linewidth=4, label="Green Vertical Line (Bold)"), "italic": plt.Line2D([0], [0], color="blue", linewidth=4, label="Blue Vertical Line (Italic)"), "underline": plt.Line2D([0], [0], color="red", linewidth=4, label="Red Vertical Line (Underline)"), } # Plot timelines with vertical lines for i, row in data.iterrows(): times = row.iloc[1:] # Subsequent columns are the timeline values row_formats = formats[i] # Formats for this row for time, fmt in zip(times, row_formats): if pd.notna(time) and fmt: # Plot only if time is valid and format is defined try: time = float(time) # Ensure the time value is numeric if fmt == "bold": # Green vertical line ax.vlines(x=time, ymin=i - arrow_width / 2, ymax=i + arrow_width / 2, color="green", linewidth=8) elif fmt == "italic": # Blue vertical line ax.vlines(x=time, ymin=i - arrow_width / 2, ymax=i + arrow_width / 2, color="blue", linewidth=11) elif fmt == "underline": # Red vertical line ax.vlines(x=time, ymin=i - arrow_width / 2, ymax=i + arrow_width / 2, color="red", linewidth=15) except ValueError: print(f"Invalid time value: {time}") # Add an arrow pointing to the 75 years label ax.annotate( "75 years\nBridge \nEnd of Life", # Text to display xy=(75, -0.5), # Arrow tip (pointing location) xytext=(60, -0.9), # Text location arrowprops=dict( facecolor='black', # Arrow color arrowstyle="->,widthA=0.5,widthB=1.0", # Arrow width properties lw=3 # Line width ), fontsize=16, ha="center" ) # Manually add unique handles for the legend legend_handles = [ plt.Line2D([0], [0], color="green", linewidth=8, label="Spot\nPainting"), plt.Line2D([0], [0], color="blue", linewidth=11, label="Over\nCoating"), plt.Line2D([0], [0], color="red", linewidth=15, label="Remove &\nReplace") ] # Add legend outside of the graph on the right ax.legend( handles=legend_handles, loc="center left", # Position legend on the right side bbox_to_anchor=(1, 0.5), # Adjust position to be centered vertically outside the graph fontsize=18, frameon=True # Add a frame around the legend ) # Add labels and title ax.set_xlabel("Time (years)", fontsize=22) ax.tick_params(axis='x', labelsize=20) # Adjust label size for the x-axis ax.tick_params(axis='y', labelsize=20) # Adjust layout and show the plot plt.tight_layout() plt.show() # ********************************************************************************************************