#import necessary modules import os import csv import numpy as np from qgis.core import QgsProject from qgis.core import QgsVectorLayer from PyQt5.QtCore import QVariant from qgis.core import QgsField, QgsExpression, QgsFeature from qgis.core import QgsVectorFileWriter #Add paths to input layers here #get the path to needed shapefiles path_to_AADFlayer="" path_to_boundaries="" path_to_pop="" path_to_popdens="" path_to_inputmatrix="" #Format for the input layers AADFlayer=QgsVectorLayer(path_to_AADFlayer,"AADF_test","ogr") WECA_boundaries=QgsVectorLayer(path_to_boundaries,"WECA","ogr") #Add paths to output layers here path_to_emissionsresults="" path_to_hhresults="" path_to_eiresults="" #read input matrix from csv file (1st - UoB PC | 2nd - Personal PC) and store values as specfic variables with open(path_to_inputmatrix, 'r') as f: reader=csv.reader(f) #read rows in csv file rows=list(reader) #store rows in a "master" list #Import emissions input values #store lists for EF and FW NOx_EF=[float(rows[2][col]) for col in [1,2,3,4,5,6,7,8,9,10]] PM25exhaust_EF=[float(rows[3][col]) for col in [1,2,3,4,5,6,7,8,9,10]] S0_PM25tbr_EF=[float(rows[4][col]) for col in [1,2,3,4,5,6,7,8,9,10]] Sx_PM25tbr_EF=[float(rows[5][col]) for col in [1,2,3,4,5,6,7,8,9,10]] #Sx_PM25_tbr - will change for electrification CO2_EF=[float(rows[6][col]) for col in [1,2,3,4,5,6,7,8,9,10]] FW_Motorway=[float(rows[7][col]) for col in [1,2,3,4,5,6,7,8,9,10]] FW_Urban=[float(rows[8][col]) for col in [1,2,3,4,5,6,7,8,9,10]] FW_Rural=[float(rows[9][col]) for col in [1,2,3,4,5,6,7,8,9,10]] u_e=[float(rows[10][col]) for col in [1,2,3,4,5,6,7,8,9,10]] r_e=[float(rows[11][col]) for col in [1,2,3,4,5,6,7,8,9,10]] reg_e=[float(rows[12][col]) for col in [1,2,3,4,5,6,7,8,9,10]] uc=[float(rows[13][col]) for col in [1,2,3,4,5,6,7,8,9,10]] rc=[float(rows[14][col]) for col in [1,2,3,4,5,6,7,8,9,10]] reg_c=[float(rows[15][col]) for col in [1,2,3,4,5,6,7,8,9,10]] T0ms_percent=[float(rows[16][col]) for col in [1,2,3,4,5,6,7,8,9,10]] urban_mschange_percent=[float(rows[17][col]) for col in [1,2,3,4,5,6,7,8,9,10]] rural_mschange_percent=[float(rows[18][col]) for col in [1,2,3,4,5,6,7,8,9,10]] #Import total trip information S0year=[float(rows[1][col]) for col in [13]] Sxyear=[float(rows[2][col]) for col in [13]] health_d=[float(rows[3][col]) for col in [13]] standard_d=[float(rows[4][col]) for col in [13]] urban_cycl_tripsperson=[float(rows[5][col]) for col in [13]] urban_passcar_tripsperson=[float(rows[6][col]) for col in [13]] urban_motorcycle_tripsperson=[float(rows[7][col]) for col in [13]] urban_bus_tripsperson=[float(rows[8][col]) for col in [13]] urban_taxi_tripsperson=[float(rows[9][col]) for col in [13]] rtown_cycl_tripsperson=[float(rows[10][col]) for col in [13]] rtown_passcar_tripsperson=[float(rows[11][col]) for col in [13]] rtown_motorcycle_tripsperson=[float(rows[12][col]) for col in [13]] rtown_bus_tripsperson=[float(rows[13][col]) for col in [13]] rtown_taxi_tripsperson=[float(rows[14][col]) for col in [13]] rvillage_cycl_tripsperson=[float(rows[15][col]) for col in [13]] rvillage_passcar_tripsperson=[float(rows[16][col]) for col in [13]] rvillage_motorcycle_tripsperson=[float(rows[17][col]) for col in [13]] rvillage_bus_tripsperson=[float(rows[18][col]) for col in [13]] rvillage_taxi_tripsperson=[float(rows[19][col]) for col in [13]] passcar_occup=[float(rows[20][col]) for col in [13]] bus_occup=[float(rows[21][col]) for col in [13]] taxi_occup=[float(rows[22][col]) for col in [13]] wpercent_cycling=[float(rows[23][col]) for col in [13]] wpercent_motorcycle=[float(rows[24][col]) for col in [13]] wpercent_passcar=[float(rows[25][col]) for col in [13]] wpercent_bus=[float(rows[26][col]) for col in [13]] wpercent_taxi=[float(rows[27][col]) for col in [13]] #Import hh input values annual_cycling=[float(rows[1][col]) for col in [17]] trip_length=[float(rows[2][col]) for col in [17]] local_pop=[float(rows[3][col]) for col in [17]] local_popdens=[float(rows[4][col]) for col in [17]] local_fatalities=[float(rows[5][col]) for col in [17]] local_accidents=[float(rows[6][col]) for col in [17]] cycling_speed=[float(rows[7][col]) for col in [17]] ref_cycling_vol=[float(rows[8][col]) for col in [17]] dr_cycling=[float(rows[9][col]) for col in [17]] dr_motorcycles=[float(rows[10][col]) for col in [17]] dr_passcars=[float(rows[11][col]) for col in [17]] dr_bus=[float(rows[12][col]) for col in [17]] seriousr_cycling=[float(rows[13][col]) for col in [17]] seriousr_motorcycles=[float(rows[14][col]) for col in [17]] seriousr_passcars=[float(rows[15][col]) for col in [17]] seriousr_bus=[float(rows[16][col]) for col in [17]] slightr_cycling=[float(rows[17][col]) for col in [17]] slightr_motorcycles=[float(rows[18][col]) for col in [17]] slightr_passcars=[float(rows[19][col]) for col in [17]] slightr_bus=[float(rows[20][col]) for col in [17]] S0_vkm_cycling=[float(rows[21][col]) for col in [17]] S0_vkm_motorcycles=[float(rows[22][col]) for col in [17]] S0_vkm_passcars=[float(rows[23][col]) for col in [17]] S0_vkm_bus=[float(rows[24][col]) for col in [17]] HR_PM25=[float(rows[25][col]) for col in [17]] HR_NOx=[float(rows[26][col]) for col in [17]] RR_cycling=[float(rows[27][col]) for col in [17]] #Import economic impacts - environmental and social - input values NOxGBP=[float(rows[2][col]) for col in [21]] PM25GBP=[float(rows[3][col]) for col in [21]] CO2GBP=[float(rows[4][col]) for col in [21]] life_expect=[float(rows[5][col]) for col in [21]] YLL_mGBP=[float(rows[6][col]) for col in [21]] fatal_net_mGBP=[float(rows[7][col]) for col in [21]] fatal_WTP_mGBP=[float(rows[8][col]) for col in [21]] fatal_medical_mGBP=[float(rows[9][col]) for col in [21]] fatal_police_mGBP=[float(rows[10][col]) for col in [21]] fatal_property_mGBP=[float(rows[11][col]) for col in [21]] fatal_insurance_mGBP=[float(rows[12][col]) for col in [21]] serious_net_mGBP=[float(rows[13][col]) for col in [21]] serious_WTP_mGBP=[float(rows[14][col]) for col in [21]] serious_medical_mGBP=[float(rows[15][col]) for col in [21]] serious_police_mGBP=[float(rows[16][col]) for col in [21]] serious_property_mGBP=[float(rows[17][col]) for col in [21]] serious_insurance_mGBP=[float(rows[18][col]) for col in [21]] slight_net_mGBP=[float(rows[19][col]) for col in [21]] slight_WTP_mGBP=[float(rows[20][col]) for col in [21]] slight_medical_mGBP=[float(rows[21][col]) for col in [21]] slight_police_mGBP=[float(rows[22][col]) for col in [21]] slight_property_mGBP=[float(rows[23][col]) for col in [21]] slight_insurance_mGBP=[float(rows[24][col]) for col in [21]] #Import economic impacts - user and provider - input values VAT=[float(rows[2][col]) for col in [24]] e_VAT=[float(rows[3][col]) for col in [24]] motor_chargeS0=[float(rows[4][col]) for col in [24]] motor_chargeSx=[float(rows[5][col]) for col in [24]] no_household_pop=[float(rows[6][col]) for col in [24]] bus_ticket_S0=[float(rows[7][col]) for col in [24]] bus_ticket_Sx=[float(rows[8][col]) for col in [24]] taxi_cost_S0=[float(rows[9][col]) for col in [24]] taxi_cost_Sx=[float(rows[10][col]) for col in [24]] workh_cycling=[float(rows[11][col]) for col in [24]] workh_passcar=[float(rows[12][col]) for col in [24]] workh_PT=[float(rows[13][col]) for col in [24]] nworkh_cycling=[float(rows[14][col]) for col in [24]] nworkh_passcar=[float(rows[15][col]) for col in [24]] nworkh_PT=[float(rows[16][col]) for col in [24]] S0_cycling_wGBP=[float(rows[17][col]) for col in [24]] S0_passcar_wGBP=[float(rows[18][col]) for col in [24]] S0_motorcycle_wGBP=[float(rows[19][col]) for col in [24]] S0_bus_wGBP=[float(rows[20][col]) for col in [24]] S0_taxi_wGBP=[float(rows[21][col]) for col in [24]] Sx_cycling_wGBP=[float(rows[22][col]) for col in [24]] Sx_passcar_wGBP=[float(rows[23][col]) for col in [24]] Sx_motorcycle_wGBP=[float(rows[24][col]) for col in [24]] Sx_bus_wGBP=[float(rows[25][col]) for col in [24]] Sx_taxi_wGBP=[float(rows[26][col]) for col in [24]] S0_nwGBP=[float(rows[27][col]) for col in [24]] Sx_nwGBP=[float(rows[28][col]) for col in [24]] #Import economic impacts - operating costs - input values car_fueltx=[float(rows[2][col]) for col in [27]] bus_fueltx=[float(rows[3][col]) for col in [27]] motorcycle_fueltx=[float(rows[4][col]) for col in [27]] fop_S0w_passcar=[float(rows[5][col]) for col in [27]] fop_S0w_epasscar=[float(rows[6][col]) for col in [27]] fop_S0w_bus=[float(rows[7][col]) for col in [27]] fop_S0w_ebus=[float(rows[8][col]) for col in [27]] fop_S0nw_passcar=[float(rows[9][col]) for col in [27]] fop_S0nw_epasscar=[float(rows[10][col]) for col in [27]] fop_Sxw_passcar=[float(rows[11][col]) for col in [27]] fop_Sxw_epasscar=[float(rows[12][col]) for col in [27]] fop_Sxw_bus=[float(rows[13][col]) for col in [27]] fop_Sxw_ebus=[float(rows[14][col]) for col in [27]] fop_Sxnw_passcar=[float(rows[15][col]) for col in [27]] fop_Sxnw_epasscar=[float(rows[16][col]) for col in [27]] nfop_S0wnw_passcar=[float(rows[17][col]) for col in [27]] nfop_Sxwnw_passcar=[float(rows[18][col]) for col in [27]] ############################################################# 1 Emissions ######################################################## #check if input layer is valid if not AADFlayer.isValid(): print("Layer failed to load") else: QgsProject.instance().addMapLayer(AADFlayer) AADFlayer.startEditing() #if layer is valid start editing #add field(s) for emissions layer_provider = AADFlayer.dataProvider() layer_provider.addAttributes([QgsField("S0_AADF_NOx", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0_AADF_PM25", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0_AADF_CO2", QVariant.Double)]) layer_provider.addAttributes([QgsField("Sx_AADF_NOx", QVariant.Double)]) layer_provider.addAttributes([QgsField("Sx_AADF_PM25", QVariant.Double)]) layer_provider.addAttributes([QgsField("Sx_AADF_CO2", QVariant.Double)]) AADFlayer.updateFields() #update the fields to the layer #start calculations (AADF x FleetWeight x EF)and round with 6 decimal places for feature in AADFlayer.getFeatures(): if feature['Classification']=='Motorway': #Motorways feature['S0_AADF_NOx']= float(round(((feature['passcars'])*((FW_Motorway[0]*NOx_EF[0])+(FW_Motorway[1]*NOx_EF[1]))) +((feature['lgvs'])*((FW_Motorway[2]*NOx_EF[2])+(FW_Motorway[3]*NOx_EF[3]))) +((feature['hgvs'])*((FW_Motorway[4]*NOx_EF[4])+(FW_Motorway[5]*NOx_EF[5]))) +(feature['motorcycles']*NOx_EF[6]) +(feature['bus']*NOx_EF[7]),6)) feature['S0_AADF_PM25'] = float(round((((feature['passcars'])*((FW_Motorway[0]*PM25exhaust_EF[0])+(FW_Motorway[1]*PM25exhaust_EF[1])))+(feature['passcars']*S0_PM25tbr_EF[1])) +(((feature['lgvs'])*((FW_Motorway[2]*PM25exhaust_EF[2])+(FW_Motorway[3]*PM25exhaust_EF[3])))+(feature['lgvs']*S0_PM25tbr_EF[2])) +(((feature['hgvs'])*((FW_Motorway[4]*PM25exhaust_EF[4])+(FW_Motorway[5]*PM25exhaust_EF[5])))+((feature['hgvs'])*((FW_Motorway[4]*S0_PM25tbr_EF[4])+(FW_Motorway[5] * S0_PM25tbr_EF[5])))) +(feature['motorcycles'])*((PM25exhaust_EF[6]+S0_PM25tbr_EF[6])) +(feature['bus'])*((PM25exhaust_EF[7]+S0_PM25tbr_EF[7])),6)) feature['S0_AADF_CO2']= float(round(((feature['passcars'])*((FW_Motorway[0]*CO2_EF[0])+(FW_Motorway[1] * CO2_EF[1]))) +((feature['lgvs'])*((FW_Motorway[2]*CO2_EF[2])+(FW_Motorway[3]*CO2_EF[3]))) +((feature['hgvs'])*((FW_Motorway[4]*CO2_EF[4])+(FW_Motorway[5]*CO2_EF[5]))) +(feature['motorcycles']*CO2_EF[6]) +(feature['bus'] *CO2_EF[7]),6)) feature['Sx_AADF_NOx']= float(round(((uc[0]*feature['passcars']*(1-u_e[0]))*((FW_Motorway[0]*NOx_EF[0])+(FW_Motorway[1]*NOx_EF[1]))) +((uc[2]*feature['lgvs']*(1-u_e[2]))*((FW_Motorway[2]*NOx_EF[2])+(FW_Motorway[3]*NOx_EF[3]))) +((uc[4]*feature['hgvs']*(1-u_e[4]))*((FW_Motorway[4]*NOx_EF[4])+(FW_Motorway[5]*NOx_EF[5]))) +(uc[6]*feature['motorcycles']*(1-u_e[6])*NOx_EF[6]) +(uc[7]*feature['bus']*(1-u_e[7])*NOx_EF[7]),6)) feature['Sx_AADF_PM25'] = float(round((((uc[0]*feature['passcars']*(1-u_e[0]))*((FW_Motorway[0]*PM25exhaust_EF[0])+(FW_Motorway[1]*PM25exhaust_EF[1])+(S0_PM25tbr_EF[1])))+(uc[0]*u_e[0]*feature['passcars']*Sx_PM25tbr_EF[1])) +(((uc[2]*feature['lgvs']*(1-u_e[2]))*((FW_Motorway[2]*PM25exhaust_EF[2])+(FW_Motorway[3]*PM25exhaust_EF[3])+S0_PM25tbr_EF[2]))+(uc[2]*u_e[2]*feature['lgvs']*Sx_PM25tbr_EF[2])) +(((uc[4]*feature['hgvs']*(1-u_e[4]))*((FW_Motorway[4]*PM25exhaust_EF[4])+(FW_Motorway[5]*PM25exhaust_EF[5])+S0_PM25tbr_EF[4]))+((uc[4]*u_e[4]*feature['hgvs'])*((FW_Motorway[4]*Sx_PM25tbr_EF[4])+(FW_Motorway[5]*Sx_PM25tbr_EF[5])))) +(((uc[6]*feature['motorcycles']*(1-u_e[6]))*(PM25exhaust_EF[6]+S0_PM25tbr_EF[6]))+(uc[6]*feature['motorcycles']*(u_e[6])*Sx_PM25tbr_EF[6])) +((((uc[7]*feature['bus'])*(1-u_e[7]))*(PM25exhaust_EF[7]+S0_PM25tbr_EF[7]))+((uc[7]*feature['bus'])*(u_e[7])*Sx_PM25tbr_EF[7])),6)) feature['Sx_AADF_CO2']= float(round(((uc[0]*feature['passcars']*(1-u_e[0]))*((FW_Motorway[0]*CO2_EF[0])+(FW_Motorway[1]*CO2_EF[1]))) +((uc[2]*feature['lgvs']*(1-u_e[2]))*((FW_Motorway[2]*CO2_EF[2])+(FW_Motorway[3]*CO2_EF[3]))) +((uc[4]*feature['hgvs']*(1-u_e[4]))*((FW_Motorway[4]*CO2_EF[4])+(FW_Motorway[5]*CO2_EF[5]))) +(uc[6]*feature['motorcycles']*(1-u_e[6])*CO2_EF[6]) +(uc[7]*feature['bus']*(1-u_e[7])*CO2_EF[7]),6)) elif feature['Classification']=='Urban': #Urban roads feature['S0_AADF_NOx']= float(round(((feature['passcars'])*((FW_Urban[0]*NOx_EF[0])+(FW_Urban[1]* NOx_EF[1]))) +((feature['lgvs'])*((FW_Urban[2]*NOx_EF[2])+(FW_Urban[3]*NOx_EF[3]))) +((feature['hgvs'])*((FW_Urban[4]*NOx_EF[4])+(FW_Urban[5]*NOx_EF[5]))) +(feature['motorcycles']*NOx_EF[6]) +(feature['bus']*NOx_EF[7]),6)) feature['S0_AADF_PM25'] = float(round ((((feature['passcars'])*((FW_Urban[0]*PM25exhaust_EF[0])+(FW_Urban[1]* PM25exhaust_EF[1])))+(feature['passcars']*S0_PM25tbr_EF[1])) +(((feature['lgvs'])*((FW_Urban[2]*PM25exhaust_EF[2])+(FW_Urban[3]*PM25exhaust_EF[3])))+(feature['lgvs']*S0_PM25tbr_EF[2])) +(((feature['hgvs'])*((FW_Urban[4]*PM25exhaust_EF[4])+(FW_Urban[5]*PM25exhaust_EF[5])))+((feature['hgvs'])*((FW_Urban[4]*S0_PM25tbr_EF[4])+(FW_Urban[5] * S0_PM25tbr_EF[5])))) +(feature['motorcycles'])*((PM25exhaust_EF[6]+S0_PM25tbr_EF[6])) +(feature['bus'])*((PM25exhaust_EF[7]+S0_PM25tbr_EF[7])),6)) feature['S0_AADF_CO2']= float(round (((feature['passcars'])*((FW_Urban[0]*CO2_EF[0])+(FW_Urban[1] * CO2_EF[1]))) +((feature['lgvs'])*((FW_Urban[2]*CO2_EF[2])+(FW_Urban[3]*CO2_EF[3]))) +((feature['hgvs'])*((FW_Urban[4]*CO2_EF[4])+(FW_Urban[5]*CO2_EF[5]))) +(feature['motorcycles']*CO2_EF[6]) +(feature['bus'] *CO2_EF[7]),6)) feature['Sx_AADF_NOx']= float(round(((uc[0]*feature['passcars']*(1-u_e[0]))*((FW_Urban[0]*NOx_EF[0])+(FW_Urban[1]*NOx_EF[1]))) +((uc[2]*feature['lgvs']*(1-u_e[2]))*((FW_Urban[2]*NOx_EF[2])+(FW_Urban[3]*NOx_EF[3]))) +((uc[4]*feature['hgvs']*(1-u_e[4]))*((FW_Urban[4]*NOx_EF[4])+(FW_Urban[5]*NOx_EF[5]))) +(uc[6]*feature['motorcycles']*(1-u_e[6])*NOx_EF[6]) +(uc[7]*feature['bus']*(1-u_e[7])*NOx_EF[7]),6)) feature['Sx_AADF_PM25'] = float(round((((uc[0]*feature['passcars']*(1-u_e[0]))*((FW_Urban[0]*PM25exhaust_EF[0])+(FW_Urban[1]*PM25exhaust_EF[1])+(S0_PM25tbr_EF[1])))+(uc[0]*u_e[0]*feature['passcars']*Sx_PM25tbr_EF[1])) +(((uc[2]*feature['lgvs']*(1-u_e[2]))*((FW_Urban[2]*PM25exhaust_EF[2])+(FW_Urban[3]*PM25exhaust_EF[3])+S0_PM25tbr_EF[2]))+(uc[2]*u_e[2]*feature['lgvs']*Sx_PM25tbr_EF[2])) +(((uc[4]*feature['hgvs']*(1-u_e[4]))*((FW_Urban[4]*PM25exhaust_EF[4])+(FW_Urban[5]*PM25exhaust_EF[5])+S0_PM25tbr_EF[4]))+((uc[4]*u_e[4]*feature['hgvs'])*((FW_Urban[4]*Sx_PM25tbr_EF[4])+(FW_Urban[5]*Sx_PM25tbr_EF[5])))) +(((uc[6]*feature['motorcycles']*(1-u_e[6]))*(PM25exhaust_EF[6]+S0_PM25tbr_EF[6]))+(uc[6]*feature['motorcycles']*(u_e[6])*Sx_PM25tbr_EF[6])) +((((uc[7]*feature['bus'])*(1-u_e[7]))*(PM25exhaust_EF[7]+S0_PM25tbr_EF[7]))+((uc[7]*feature['bus'])*(u_e[7])*Sx_PM25tbr_EF[7])),6)) feature['Sx_AADF_CO2']= float(round(((uc[0]*feature['passcars']*(1-u_e[0]))*((FW_Urban[0]*CO2_EF[0])+(FW_Urban[1]*CO2_EF[1]))) +((uc[2]*feature['lgvs']*(1-u_e[2]))*((FW_Urban[2]*CO2_EF[2])+(FW_Urban[3]*CO2_EF[3]))) +((uc[4]*feature['hgvs']*(1-u_e[4]))*((FW_Urban[4]*CO2_EF[4])+(FW_Urban[5]*CO2_EF[5]))) +(uc[6]*feature['motorcycles']*(1-u_e[6])*CO2_EF[6]) +(uc[7]*feature['bus']*(1-u_e[7]) *CO2_EF[7]),6)) elif feature['Classification']=='Rural': #Rural roads feature['S0_AADF_NOx']= float(round(((feature['passcars'])*((FW_Rural[0]*NOx_EF[0])+(FW_Rural[1] * NOx_EF[1]))) +((feature['lgvs'])*((FW_Rural[2]*NOx_EF[2])+(FW_Rural[3]*NOx_EF[3]))) +((feature['hgvs'])*((FW_Rural[4]*NOx_EF[4])+(FW_Rural[5]*NOx_EF[5]))) +(feature['motorcycles']*NOx_EF[6]) +(feature['bus']*NOx_EF[7]),6)) feature['S0_AADF_PM25'] = float(round((((feature['passcars'])*((FW_Rural[0]*PM25exhaust_EF[0])+(FW_Rural[1]* PM25exhaust_EF[1])))+(feature['passcars']*S0_PM25tbr_EF[1])) +(((feature['lgvs'])*((FW_Rural[2]*PM25exhaust_EF[2])+(FW_Rural[3]*PM25exhaust_EF[3])))+(feature['lgvs']*S0_PM25tbr_EF[2])) +(((feature['hgvs'])*((FW_Rural[4]*PM25exhaust_EF[4])+(FW_Rural[5]*PM25exhaust_EF[5])))+((feature['hgvs'])*((FW_Rural[4]*S0_PM25tbr_EF[4])+(FW_Rural[5] * S0_PM25tbr_EF[5])))) +(feature['motorcycles'])*((PM25exhaust_EF[6]+S0_PM25tbr_EF[6])) +(feature['bus'])*((PM25exhaust_EF[7]+S0_PM25tbr_EF[7])),5)) feature['S0_AADF_CO2']= float(round(((feature['passcars'])*((FW_Rural[0]*CO2_EF[0])+(FW_Rural[1]*CO2_EF[1]))) +((feature['lgvs'])*((FW_Rural[2]*CO2_EF[2])+(FW_Rural[3]*CO2_EF[3]))) +((feature['hgvs'])*((FW_Rural[4]*CO2_EF[4])+(FW_Rural[5]*CO2_EF[5]))) +(feature['motorcycles']*CO2_EF[6]) +(feature['bus'] *CO2_EF[7]),6)) feature['Sx_AADF_NOx']= float(round(((rc[0]*feature['passcars']*(1-r_e[0]))*((FW_Rural[0]*NOx_EF[0])+(FW_Rural[1]*NOx_EF[1]))) +((rc[2]*feature['lgvs']*(1-r_e[2]))*((FW_Rural[2]*NOx_EF[2])+(FW_Rural[3]*NOx_EF[3]))) +((rc[4]*feature['hgvs']*(1-r_e[4]))*((FW_Rural[4]*NOx_EF[4])+(FW_Rural[5]*NOx_EF[5]))) +(rc[6]*feature['motorcycles']*(1-r_e[6])*NOx_EF[6]) +(rc[7]*feature['bus']*(1-r_e[7])*NOx_EF[7]),6)) feature['Sx_AADF_PM25'] = float(round((((rc[0]*feature['passcars']*(1-r_e[0]))*((FW_Rural[0]*PM25exhaust_EF[0])+(FW_Rural[1]*PM25exhaust_EF[1])+(S0_PM25tbr_EF[1])))+(rc[0]*r_e[0]*feature['passcars']*Sx_PM25tbr_EF[1])) +(((rc[2]*feature['lgvs']*(1-r_e[2]))*((FW_Rural[2]*PM25exhaust_EF[2])+(FW_Rural[3]*PM25exhaust_EF[3])+S0_PM25tbr_EF[2]))+(rc[2]*r_e[2]*feature['lgvs']*Sx_PM25tbr_EF[2])) +(((rc[4]*feature['hgvs']*(1-r_e[4]))*((FW_Rural[4]*PM25exhaust_EF[4])+(FW_Rural[5]*PM25exhaust_EF[5])+S0_PM25tbr_EF[4]))+((rc[4]*r_e[4]*feature['hgvs'])*((FW_Rural[4]*Sx_PM25tbr_EF[4])+(FW_Rural[5]*Sx_PM25tbr_EF[5])))) +(((rc[6]*feature['motorcycles']*(1-r_e[6]))*(PM25exhaust_EF[6]+S0_PM25tbr_EF[6]))+(rc[6]*feature['motorcycles']*(r_e[6])*Sx_PM25tbr_EF[6])) +((((rc[7]*feature['bus'])*(1-r_e[7]))*(PM25exhaust_EF[7]+S0_PM25tbr_EF[7]))+((rc[7]*feature['bus'])*(r_e[7])*Sx_PM25tbr_EF[7])),6)) feature['Sx_AADF_CO2']= float(round(((rc[0]*feature['passcars']*(1-r_e[0]))*((FW_Rural[0]*CO2_EF[0])+(FW_Rural[1]*CO2_EF[1]))) +((rc[2]*feature['lgvs']*(1-r_e[2]))*((FW_Rural[2]*CO2_EF[2])+(FW_Rural[3]*CO2_EF[3]))) +((rc[4]*feature['hgvs']*(1-r_e[4]))*((FW_Rural[4]*CO2_EF[4])+(FW_Rural[5]*CO2_EF[5]))) +(rc[6]*feature['motorcycles']*(1-r_e[6])*CO2_EF[6]) +(rc[7]*feature['bus']*(1-r_e[7])*CO2_EF[7]),6)) AADFlayer.updateFeature(feature) #update layer changes AADFlayer.commitChanges() #record layer changes and stop editing #IDW Interpolation of AADFxEFxFW #INTERPOLATION_DATA format interpolation_dataS0AADFNOx = '{}::~::0::~::{}::~::0'.format(AADFlayer.source(), AADFlayer.fields().indexFromName('S0_AADF_NOx')) interpolation_dataS0AADFPM25 = '{}::~::0::~::{}::~::0'.format(AADFlayer.source(), AADFlayer.fields().indexFromName('S0_AADF_PM25')) interpolation_dataS0AADFCO2 = '{}::~::0::~::{}::~::0'.format(AADFlayer.source(), AADFlayer.fields().indexFromName('S0_AADF_CO2')) interpolation_dataSxAADFNOx = '{}::~::0::~::{}::~::0'.format(AADFlayer.source(), AADFlayer.fields().indexFromName('Sx_AADF_NOx')) interpolation_dataSxAADFPM25 = '{}::~::0::~::{}::~::0'.format(AADFlayer.source(), AADFlayer.fields().indexFromName('Sx_AADF_PM25')) interpolation_dataSxAADFCO2 = '{}::~::0::~::{}::~::0'.format(AADFlayer.source(), AADFlayer.fields().indexFromName('Sx_AADF_CO2')) #IDW function for emissions #S0 resultS0AADFNOx=processing.run("qgis:idwinterpolation", {'INTERPOLATION_DATA':interpolation_dataS0AADFNOx, 'DISTANCE_COEFFICIENT':3, 'EXTENT':'346395.8385,382607.1013,152770.0196,198319.9009 [EPSG:27700]', 'PIXEL_SIZE':1000, 'OUTPUT':'TEMPORARY_OUTPUT'}) S0_AADF_NOx=resultS0AADFNOx['OUTPUT'] resultS0AADFPM25=processing.run("qgis:idwinterpolation", {'INTERPOLATION_DATA':interpolation_dataS0AADFPM25, 'DISTANCE_COEFFICIENT':3, 'EXTENT':'346395.8385,382607.1013,152770.0196,198319.9009 [EPSG:27700]', 'PIXEL_SIZE':1000, 'OUTPUT':'TEMPORARY_OUTPUT'}) S0_AADF_PM25=resultS0AADFPM25['OUTPUT'] resultS0AADFCO2=processing.run("qgis:idwinterpolation", {'INTERPOLATION_DATA':interpolation_dataS0AADFCO2, 'DISTANCE_COEFFICIENT':3, 'EXTENT':'346395.8385,382607.1013,152770.0196,198319.9009 [EPSG:27700]', 'PIXEL_SIZE':1000, 'OUTPUT':'TEMPORARY_OUTPUT'}) S0_AADF_CO2=resultS0AADFCO2['OUTPUT'] #Sx resultSxAADFNOx=processing.run("qgis:idwinterpolation", {'INTERPOLATION_DATA':interpolation_dataSxAADFNOx, 'DISTANCE_COEFFICIENT':3, 'EXTENT':'346395.8385,382607.1013,152770.0196,198319.9009 [EPSG:27700]', 'PIXEL_SIZE':1000, 'OUTPUT':'TEMPORARY_OUTPUT'}) Sx_AADF_NOx=resultSxAADFNOx['OUTPUT'] resultSxAADFPM25=processing.run("qgis:idwinterpolation", {'INTERPOLATION_DATA':interpolation_dataSxAADFPM25, 'DISTANCE_COEFFICIENT':3, 'EXTENT':'346395.8385,382607.1013,152770.0196,198319.9009 [EPSG:27700]', 'PIXEL_SIZE':1000, 'OUTPUT':'TEMPORARY_OUTPUT'}) Sx_AADF_PM25=resultSxAADFPM25['OUTPUT'] resultSxAADFCO2=processing.run("qgis:idwinterpolation", {'INTERPOLATION_DATA':interpolation_dataSxAADFCO2, 'DISTANCE_COEFFICIENT':3, 'EXTENT':'346395.8385,382607.1013,152770.0196,198319.9009 [EPSG:27700]', 'PIXEL_SIZE':1000, 'OUTPUT':'TEMPORARY_OUTPUT'}) Sx_AADF_CO2=resultSxAADFCO2['OUTPUT'] #Prepare S0_AADF grid S0_AADF_CO2_points=processing.run("native:pixelstopoints", {'INPUT_RASTER':S0_AADF_CO2, 'RASTER_BAND':1, 'FIELD_NAME':'S0_AADF_CO2', 'OUTPUT':'TEMPORARY_OUTPUT'}) S0_AADF_CO2_points_output=S0_AADF_CO2_points['OUTPUT'] #Add coordinates S0_AADF_CO2_coord=processing.run("qgis:exportaddgeometrycolumns", {'INPUT':S0_AADF_CO2_points_output, 'CALC_METHOD':0, 'OUTPUT':'TEMPORARY_OUTPUT'}) S0_AADF_CO2_coord_output=S0_AADF_CO2_coord['OUTPUT'] #Buffer S0_AADF_CO2_grid=processing.run("native:buffer", {'INPUT':S0_AADF_CO2_coord_output, 'DISTANCE':500,'SEGMENTS':5, 'END_CAP_STYLE':2,'JOIN_STYLE':0, 'MITER_LIMIT':2,'DISSOLVE':False, 'OUTPUT':'TEMPORARY_OUTPUT'}) S0_AADFCO2_grid=S0_AADF_CO2_grid['OUTPUT'] #Intersection S0_AADF_CO2_int=processing.run("native:intersection", {'INPUT':S0_AADFCO2_grid, 'OVERLAY':WECA_boundaries, 'INPUT_FIELDS':['fid','xcoord','ycoord','S0_AADF_CO2'], 'OVERLAY_FIELDS':['NAME'], 'OVERLAY_FIELDS_PREFIX':'', 'OUTPUT':'TEMPORARY_OUTPUT', 'GRID_SIZE':None}) S0_AADFCO2grid_output=S0_AADF_CO2_int['OUTPUT'] #Add Pop., Pop.density and other AADFs #Pop density and Pop. S0_grid_pop=processing.run("native:zonalstatisticsfb",{'INPUT':S0_AADFCO2grid_output, 'INPUT_RASTER':path_to_pop, 'RASTER_BAND':1, 'COLUMN_PREFIX':'Pop.', 'STATISTICS':[1], 'OUTPUT':'TEMPORARY_OUTPUT'}) S0_pop_output=S0_grid_pop['OUTPUT'] S0_grid_popdens=processing.run("native:zonalstatisticsfb",{'INPUT':S0_pop_output, 'INPUT_RASTER':path_to_popdens, 'RASTER_BAND':1, 'COLUMN_PREFIX':'PopDens', 'STATISTICS':[2], 'OUTPUT':'TEMPORARY_OUTPUT'}) S0_AADF_pop=S0_grid_popdens['OUTPUT'] #S0 S0_AADF_PM25_ZS=processing.run("native:zonalstatisticsfb", {'INPUT':S0_AADF_pop, 'INPUT_RASTER':S0_AADF_PM25, 'RASTER_BAND':1, 'COLUMN_PREFIX':'S0_AADF_PM25', 'STATISTICS':[2], 'OUTPUT':'TEMPORARY_OUTPUT'}) S0_AADF_PM25_grid=S0_AADF_PM25_ZS['OUTPUT'] S0_AADF_NOx_ZS=processing.run("native:zonalstatisticsfb", {'INPUT':S0_AADF_PM25_grid, 'INPUT_RASTER':S0_AADF_NOx, 'RASTER_BAND':1, 'COLUMN_PREFIX':'S0_AADF_NOx', 'STATISTICS':[2], 'OUTPUT':'TEMPORARY_OUTPUT'}) S0_AADF_NOx_grid=S0_AADF_NOx_ZS['OUTPUT'] #Sx Sx_AADF_CO2_ZS=processing.run("native:zonalstatisticsfb", {'INPUT':S0_AADF_NOx_grid, 'INPUT_RASTER':Sx_AADF_CO2, 'RASTER_BAND':1, 'COLUMN_PREFIX':'Sx_AADF_CO2', 'STATISTICS':[2], 'OUTPUT':'TEMPORARY_OUTPUT'}) Sx_AADF_CO2_grid=Sx_AADF_CO2_ZS['OUTPUT'] Sx_AADF_PM25_ZS=processing.run("native:zonalstatisticsfb", {'INPUT':Sx_AADF_CO2_grid, 'INPUT_RASTER':Sx_AADF_PM25, 'RASTER_BAND':1, 'COLUMN_PREFIX':'Sx_AADF_PM25', 'STATISTICS':[2], 'OUTPUT':'TEMPORARY_OUTPUT'}) Sx_AADF_PM25_grid=Sx_AADF_PM25_ZS['OUTPUT'] Sx_AADF_NOx_ZS=processing.run("native:zonalstatisticsfb", {'INPUT':Sx_AADF_PM25_grid, 'INPUT_RASTER':resultSxAADFNOx['OUTPUT'], 'RASTER_BAND':1, 'COLUMN_PREFIX':'Sx_AADF_NOx', 'STATISTICS':[2], 'OUTPUT':'TEMPORARY_OUTPUT'}) Sx_AADF_NOx_grid=Sx_AADF_NOx_ZS['OUTPUT'] #Edit AADFxEFxFW layer to calculate final emissions Sx_AADF_NOx_grid.startEditing() #if layer is valid start editing #delete NULL fields in Pop. feature Sx_AADF_NOx_grid.selectByExpression ('"Pop.sum" is NULL or "PopDensmean" is NULL') for feature in Sx_AADF_NOx_grid.selectedFeatures(): Sx_AADF_NOx_grid.deleteFeature(feature.id()) #delete features when Pop.sum = 0 Sx_AADF_NOx_grid.selectByExpression ('"Pop.sum"=0') for feature in Sx_AADF_NOx_grid.selectedFeatures(): Sx_AADF_NOx_grid.deleteFeature(feature.id()) #add field(s) layer_provider = Sx_AADF_NOx_grid.dataProvider() layer_provider.addAttributes([QgsField("lower_S0_NOx", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0_NOx", QVariant.Double)]) layer_provider.addAttributes([QgsField("upper_S0_NOx", QVariant.Double)]) layer_provider.addAttributes([QgsField("lower_S0_PM25", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0_PM25", QVariant.Double)]) layer_provider.addAttributes([QgsField("upper_S0_PM25", QVariant.Double)]) layer_provider.addAttributes([QgsField("lower_S0_CO2", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0_CO2", QVariant.Double)]) layer_provider.addAttributes([QgsField("upper_S0_CO2", QVariant.Double)]) layer_provider.addAttributes([QgsField("lower_NOx", QVariant.Double)]) layer_provider.addAttributes([QgsField("Sx_NOx", QVariant.Double)]) layer_provider.addAttributes([QgsField("upper_NOx", QVariant.Double)]) layer_provider.addAttributes([QgsField("lower_PM25", QVariant.Double)]) layer_provider.addAttributes([QgsField("Sx_PM25", QVariant.Double)]) layer_provider.addAttributes([QgsField("upper_PM25", QVariant.Double)]) layer_provider.addAttributes([QgsField("lower_CO2", QVariant.Double)]) layer_provider.addAttributes([QgsField("Sx_CO2", QVariant.Double)]) layer_provider.addAttributes([QgsField("upper_CO2", QVariant.Double)]) Sx_AADF_NOx_grid.updateFields() #update the fields to the layer #Final emissions using equations from regression equations for feature in Sx_AADF_NOx_grid.getFeatures(): #for the base scenario (S0) feature['S0_NOx']=float(round(8.640721244 + (0.000918*feature['S0_AADF_NOxmean'])+(0.175504*feature['PopDensmean']),6)) feature['S0_PM25']=float(round(7.737099 + (0.001385*feature['S0_AADF_PM25mean'])+(0.022317*feature['PopDensmean']),6)) feature['S0_CO2']=float(round((0.199916688*feature['S0_AADF_CO2'])+(14.69986462*feature['PopDensmean']),6)) #for scenario Sx feature['Sx_NOx']=float(round(8.640721244 + (0.000918*feature['Sx_AADF_NOxmean'])+(0.175504*feature['PopDensmean']),6)) feature['Sx_PM25']=float(round(7.737099 + (0.001385*feature['Sx_AADF_PM25mean'])+(0.022317*feature['PopDensmean']),6)) feature['Sx_CO2']=float(round((0.199916688*feature['Sx_AADF_CO2mean'])+(14.69986462*feature['PopDensmean']),6)) #lower 95% feature['lower_S0_NOx']=float(round(8.640721244 + (0.000824727*feature['S0_AADF_NOxmean'])+(0.165638391 *feature['PopDensmean']),6)) feature['lower_NOx']=float(round(8.640721244 + (0.000824727*feature['Sx_AADF_NOxmean'])+(0.165638391*feature['PopDensmean']),6)) feature['lower_S0_PM25']=float(round(7.692567129 + (0.001272817*feature['S0_AADF_PM25mean'])+(0.020968707*feature['PopDensmean']),6)) feature['lower_PM25']=float(round(7.692567129 + (0.001272817*feature['Sx_AADF_PM25mean'])+(0.020968707*feature['PopDensmean']),6)) feature['lower_S0_CO2']=float(round((0.183575664*feature['S0_AADF_CO2'])+(13.05462129*feature['PopDensmean']),6)) feature['lower_CO2']=float(round((0.183575664*feature['Sx_AADF_CO2mean'])+(13.05462129*feature['PopDensmean']),6)) #upper 95% feature['upper_S0_NOx']=float(round(9.361432756 + (0.001011273*feature['S0_AADF_NOxmean'])+(0.185369609*feature['PopDensmean']),6)) feature['upper_NOx']=float(round(9.361432756 + (0.001011273*feature['Sx_AADF_NOxmean'])+(0.185369609*feature['PopDensmean']),6)) feature['upper_S0_PM25']=float(round(7.781630871 + (0.001497183*feature['S0_AADF_PM25mean'])+(0.023665293*feature['PopDensmean']),6)) feature['upper_PM25']=float(round(7.781630871 + (0.001497183*feature['Sx_AADF_PM25mean'])+(0.023665293*feature['PopDensmean']),6)) feature['upper_S0_CO2']=float(round((0.216257711*feature['S0_AADF_CO2'])+(16.34510795*feature['PopDensmean']),6)) feature['upper_CO2']=float(round((0.216257711*feature['Sx_AADF_CO2mean'])+(16.34510795*feature['PopDensmean']),6)) Sx_AADF_NOx_grid.updateFeature(feature) #update layer changes Sx_AADF_NOx_grid.commitChanges() #record layer changes and stop editing #Clean layer Sx_AADF_NOx_grid.startEditing() #start editing Sx_AADF_NOx_grid.dataProvider().deleteAttributes([2,5,6,7,8,9]) Sx_AADF_NOx_grid.commitChanges() #record layer changes and stop editing #Add emissions results to QGIS emissions_layername="emissions" #define new name for layer Sx_AADF_NOx_grid.setName(emissions_layername) QgsProject.instance().addMapLayer(Sx_AADF_NOx_grid) #Add to QGIS #Duplicate layer for HH hhlayer = Sx_AADF_NOx_grid.clone() #duplicate emissions hh_layername="human_healh" #define new name for layer hhlayer.setName(hh_layername) QgsProject.instance().addMapLayer(hhlayer) #Add to QGIS ############################################################# 2 Human Health ######################################################## #Edit layer and add fields hhlayer.startEditing() #if layer is valid start editing #add field(s) layer_provider = hhlayer.dataProvider() layer_provider.addAttributes([QgsField("S0_mortality", QVariant.Double)]) layer_provider.addAttributes([QgsField("NOx_S0x", QVariant.Double)]) layer_provider.addAttributes([QgsField("PM25_S0x", QVariant.Double)]) layer_provider.addAttributes([QgsField("CO2_S0x", QVariant.Double)]) layer_provider.addAttributes([QgsField("Sx0_AP_mortality", QVariant.Double)]) layer_provider.addAttributes([QgsField("mins_cycling", QVariant.Double)]) layer_provider.addAttributes([QgsField("Sx0_cycl_mortality", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0_accidents", QVariant.Double)]) layer_provider.addAttributes([QgsField("Sx0_fatalacc", QVariant.Double)]) layer_provider.addAttributes([QgsField("Sx0_seriousacc", QVariant.Double)]) layer_provider.addAttributes([QgsField("Sx0_slightacc", QVariant.Double)]) layer_provider.addAttributes([QgsField("Sx0_mortality", QVariant.Double)]) hhlayer.updateFields() #update the fields to the layer #start calculations (for new features)and round with 3 decimal places for feature in hhlayer.getFeatures(): #Calculate mortality in each 1x1km square feature['S0_mortality']=float(local_fatalities[0]*(feature['Pop.sum']/local_pop[0])) feature['S0_accidents']=float(local_accidents[0]*(feature['Pop.sum']/local_pop[0])) #Air Pollution and mortality feature['NOx_S0x']=feature['S0_NOx']-feature['Sx_NOx'] feature['PM25_S0x']=feature['S0_PM25']-feature['Sx_PM25'] feature['CO2_S0x']=feature['S0_CO2']-feature['Sx_CO2'] feature['Sx0_AP_mortality']=float(((np.exp((np.log(HR_NOx[0])*(feature['Sx_NOx']-feature['S0_NOx']))+(np.log(HR_PM25[0])*(feature['Sx_PM25']-feature['S0_PM25']))))*feature['S0_mortality'])-feature['S0_mortality']) popdens_ratio=float(feature['PopDensmean']/local_popdens[0]) #Differentiating between urban and rural areas if 0.34 <= feature['PopDensmean']<=7.19: #rural areas #cycling benefits + 95% CI local_cycling_Sx=float(((annual_cycling[0]*rc[8]*trip_length[0]*60)/(52*local_pop[0]*cycling_speed[0]))*popdens_ratio) local_cycling_S0=float(((annual_cycling[0]*trip_length[0]*60)/(52*local_pop[0]*cycling_speed[0]))*popdens_ratio) feature['mins_cycling']=(local_cycling_Sx-local_cycling_S0)/ref_cycling_vol[0] feature['Sx0_cycl_mortality']=float(((np.exp((np.log(RR_cycling[0]))*(feature['mins_cycling'])))*feature['S0_mortality'])-feature['S0_mortality']) #accidental mortality fx_f0_cycling=float(round((dr_cycling[0]*S0_vkm_cycling[0]*rc[8]*0.000000001)-(dr_cycling[0]*S0_vkm_cycling[0]*0.000000001),3)) fx_f0_motorcycles=float(round((dr_motorcycles[0]*S0_vkm_motorcycles[0]*rc[6]*0.000000001)-(dr_motorcycles[0]*S0_vkm_motorcycles[0]*0.000000001),3)) fx_f0_passcars=float(round((dr_passcars[0]*S0_vkm_passcars[0]*rc[0]*0.000000001)-(dr_passcars[0]*S0_vkm_passcars[0]*0.000000001),3)) fx_f0_bus=float(round((dr_bus[0]*S0_vkm_bus[0]*rc[7]*0.000000001)-(dr_bus[0]*S0_vkm_bus[0]*0.000000001),3)) rr_facc=float((local_fatalities[0]+((fx_f0_cycling)+(fx_f0_motorcycles)+(fx_f0_passcars)+(fx_f0_bus)))/local_fatalities[0]) feature['Sx0_fatalacc']=float(feature['S0_mortality']*((rr_facc-1)/rr_facc)) #serious accidents serx_ser0_cycling=float(round((seriousr_cycling[0]*S0_vkm_cycling[0]*rc[8]*0.000000001)-(seriousr_cycling[0]*S0_vkm_cycling[0]*0.000000001),3)) serx_ser0_motorcycles=float(round((seriousr_motorcycles[0]*S0_vkm_motorcycles[0]*rc[6]*0.000000001)-(seriousr_motorcycles[0]*S0_vkm_motorcycles[0]*0.000000001),3)) serx_ser0_passcars=float(round((seriousr_passcars[0]*S0_vkm_passcars[0]*rc[0]*0.000000001)-(seriousr_passcars[0]*S0_vkm_passcars[0]*0.000000001),3)) serx_ser0_bus=float(round((seriousr_bus[0]*S0_vkm_bus[0]*rc[7]*0.000000001)-(seriousr_bus[0]*S0_vkm_bus[0]*0.000000001),3)) rr_seracc=float((local_accidents[0]+((serx_ser0_cycling)+(serx_ser0_motorcycles)+(serx_ser0_passcars)+(serx_ser0_bus)))/local_accidents[0]) feature['Sx0_seriousacc']=float(feature['S0_accidents']*((rr_seracc-1)/rr_seracc)) #slight accidents slhx_slh0_cycling=float(round((slightr_cycling[0]*S0_vkm_cycling[0]*rc[8]*0.000000001)-(slightr_cycling[0]*S0_vkm_cycling[0]*0.000000001),3)) slhx_slh0_motorcycles=float(round((slightr_motorcycles[0]*S0_vkm_motorcycles[0]*rc[6]*0.000000001)-(slightr_motorcycles[0]*S0_vkm_motorcycles[0]*0.000000001),3)) slhx_slh0_passcars=float(round((slightr_passcars[0]*S0_vkm_passcars[0]*rc[0]*0.000000001)-(slightr_passcars[0]*S0_vkm_passcars[0]*0.000000001),3)) slhx_slh0_bus=float(round((slightr_bus[0]*S0_vkm_bus[0]*rc[7]*0.000000001)-(slightr_bus[0]*S0_vkm_bus[0]*0.000000001),3)) rr_slhacc=float((local_accidents[0]+((slhx_slh0_cycling)+(slhx_slh0_motorcycles)+(slhx_slh0_passcars)+(slhx_slh0_bus)))/local_accidents[0]) feature['Sx0_slightacc']=float(feature['S0_accidents']*((rr_slhacc-1)/rr_slhacc)) elif feature['PopDensmean']>7.19: #urban areas #cycling benefits local_cycling_Sx=float(((annual_cycling[0]*uc[8]*trip_length[0]*60)/(52*local_pop[0]*cycling_speed[0]))*popdens_ratio) local_cycling_S0=float(((annual_cycling[0]*trip_length[0]*60)/(52*local_pop[0]*cycling_speed[0]))*popdens_ratio) feature['mins_cycling']=(local_cycling_Sx-local_cycling_S0)/ref_cycling_vol[0] feature['Sx0_cycl_mortality']=float(((np.exp((np.log(RR_cycling[0]))*(feature['mins_cycling'])))*feature['S0_mortality'])-feature['S0_mortality']) #accidental mortality fx_f0_cycling=float(round((dr_cycling[0]*S0_vkm_cycling[0]*uc[8]*0.000000001)-(dr_cycling[0]*S0_vkm_cycling[0]*0.000000001),3)) fx_f0_motorcycles=float(round((dr_motorcycles[0]*S0_vkm_motorcycles[0]*uc[6]*0.000000001)-(dr_motorcycles[0]*S0_vkm_motorcycles[0]*0.000000001),3)) fx_f0_passcars=float(round((dr_passcars[0]*S0_vkm_passcars[0]*uc[0]*0.000000001)-(dr_passcars[0]*S0_vkm_passcars[0]*0.000000001),3)) fx_f0_bus=float(round((dr_bus[0]*S0_vkm_bus[0]*uc[7]*0.000000001)-(dr_bus[0]*S0_vkm_bus[0]*0.000000001),3)) rr_facc=float((local_fatalities[0]+((fx_f0_cycling)+(fx_f0_motorcycles)+(fx_f0_passcars)+(fx_f0_bus)))/local_fatalities[0]) feature['Sx0_fatalacc']=float(feature['S0_mortality']*((rr_facc-1)/rr_facc)) #serious accidents serx_ser0_cycling=float(round((seriousr_cycling[0]*S0_vkm_cycling[0]*uc[8]*0.000000001)-(seriousr_cycling[0]*S0_vkm_cycling[0]*0.000000001),3)) serx_ser0_motorcycles=float(round((seriousr_motorcycles[0]*S0_vkm_motorcycles[0]*uc[6]*0.000000001)-(seriousr_motorcycles[0]*S0_vkm_motorcycles[0]*0.000000001),3)) serx_ser0_passcars=float(round((seriousr_passcars[0]*S0_vkm_passcars[0]*uc[0]*0.000000001)-(seriousr_passcars[0]*S0_vkm_passcars[0]*0.000000001),3)) serx_ser0_bus=float(round((seriousr_bus[0]*S0_vkm_bus[0]*uc[7]*0.000000001)-(seriousr_bus[0]*S0_vkm_bus[0]*0.000000001),3)) rr_seracc=float((local_accidents[0]+((serx_ser0_cycling)+(serx_ser0_motorcycles)+(serx_ser0_passcars)+(serx_ser0_bus)))/local_accidents[0]) feature['Sx0_seriousacc']=float(feature['S0_accidents']*((rr_seracc-1)/rr_seracc)) #slight accidents slhx_slh0_cycling=float(round((slightr_cycling[0]*S0_vkm_cycling[0]*uc[8]*0.000000001)-(slightr_cycling[0]*S0_vkm_cycling[0]*0.000000001),3)) slhx_slh0_motorcycles=float(round((slightr_motorcycles[0]*S0_vkm_motorcycles[0]*uc[6]*0.000000001)-(slightr_motorcycles[0]*S0_vkm_motorcycles[0]*0.000000001),3)) slhx_slh0_passcars=float(round((slightr_passcars[0]*S0_vkm_passcars[0]*uc[0]*0.000000001)-(slightr_passcars[0]*S0_vkm_passcars[0]*0.000000001),3)) slhx_slh0_bus=float(round((slightr_bus[0]*S0_vkm_bus[0]*uc[7]*0.000000001)-(slightr_bus[0]*S0_vkm_bus[0]*0.000000001),3)) rr_slhacc=float((local_accidents[0]+((slhx_slh0_cycling)+(slhx_slh0_motorcycles)+(slhx_slh0_passcars)+(slhx_slh0_bus)))/local_accidents[0]) feature['Sx0_slightacc']=float(feature['S0_accidents']*((rr_slhacc-1)/rr_slhacc)) #Calculate final Sx mortality feature['Sx0_mortality']=float(feature['Sx0_AP_mortality']+feature['Sx0_cycl_mortality']+feature['Sx0_fatalacc']) hhlayer.updateFeature(feature) #update layer changes hhlayer.commitChanges() #record layer changes and stop editing #Clean human health layer hhlayer.startEditing() #start editing hhlayer.dataProvider().deleteAttributes([4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]) #delete emissions fields hhlayer.commitChanges() #record layer changes and stop editing #Duplicate layer for economic impacts economic_impacts= hhlayer.clone() #duplicate human health ei_layername="economic_impacts" #define new name for layer economic_impacts.setName(ei_layername) QgsProject.instance().addMapLayer(economic_impacts) #Add to QGIS ############################################################# 3 Economic Impacts ######################################################## #Edit layer and add fields economic_impacts.startEditing() #if layer is valid start editing #add field(s) layer_provider = economic_impacts.dataProvider() layer_provider.addAttributes([QgsField("S0x_NOx_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0x_PM25_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0x_CO2_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0x_cyc_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0x_acnt_netout_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0x_acnt_medical_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0x_acnt_WTP_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0x_acnt_police_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0x_acnt_property_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("S0x_acnt_insurance_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("vehT0cycl", QVariant.Double)]) layer_provider.addAttributes([QgsField("vehT0passcar", QVariant.Double)]) layer_provider.addAttributes([QgsField("vehT0motorcycl", QVariant.Double)]) layer_provider.addAttributes([QgsField("vehT0bus", QVariant.Double)]) layer_provider.addAttributes([QgsField("vehT0taxi", QVariant.Double)]) layer_provider.addAttributes([QgsField("vehTxcycl", QVariant.Double)]) layer_provider.addAttributes([QgsField("vehTxpasscar", QVariant.Double)]) layer_provider.addAttributes([QgsField("vehTxmotorcycl", QVariant.Double)]) layer_provider.addAttributes([QgsField("vehTxbus", QVariant.Double)]) layer_provider.addAttributes([QgsField("vehTxtaxi", QVariant.Double)]) layer_provider.addAttributes([QgsField("user_charges_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("PTprovider_revenue_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("travel_time_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("PT_operational_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("private_operational_m£", QVariant.Double)]) layer_provider.addAttributes([QgsField("NPV_m£", QVariant.Double)]) economic_impacts.updateFields() #update the fields to the layer #start calculations (for new features)and round with 3 decimal places for feature in economic_impacts.getFeatures(): #price discounts health_disc=float((1/(1+(health_d[0]/100)))**(Sxyear[0]-S0year[0])) standard_disc=float((1/(1+(standard_d[0]/100)))**(Sxyear[0]-S0year[0])) #Environmental_social #Calculate AP costs/benefits feature['S0x_NOx_m£']=float((feature['NOx_S0x']*feature['Pop.sum']*NOxGBP[0])*health_disc) feature['S0x_PM25_m£']=float((feature['PM25_S0x']*feature['Pop.sum']*PM25GBP[0])*health_disc) feature['S0x_CO2_m£']=float(((feature['CO2_S0x']*3.67*CO2GBP[0])/feature['Pop.sum'])*health_disc) #Calculate active travel costs/benefits feature['S0x_cyc_m£']=float(((-feature['Sx0_cycl_mortality'])*life_expect[0]*YLL_mGBP[0])*health_disc) #Calculate accidental categories' costs feature['S0x_acnt_netout_m£']=float(((-feature['Sx0_fatalacc']*fatal_net_mGBP[0])+(-feature['Sx0_seriousacc']*serious_net_mGBP[0])+(-feature['Sx0_slightacc']*slight_net_mGBP[0]))*standard_disc) feature['S0x_acnt_medical_m£']=float(((-feature['Sx0_fatalacc']*fatal_medical_mGBP[0])+(-feature['Sx0_seriousacc']*serious_medical_mGBP[0])+(-feature['Sx0_slightacc']*slight_medical_mGBP[0]))*standard_disc) feature['S0x_acnt_WTP_m£']=float(((-feature['Sx0_fatalacc']*fatal_WTP_mGBP[0])+(-feature['Sx0_seriousacc']*serious_WTP_mGBP[0])+(-feature['Sx0_slightacc']*slight_WTP_mGBP[0]))*health_disc) feature['S0x_acnt_police_m£']=float(((-feature['Sx0_fatalacc']*fatal_police_mGBP[0])+(-feature['Sx0_seriousacc']*serious_police_mGBP[0])+(-feature['Sx0_slightacc']*slight_police_mGBP[0]))*standard_disc) feature['S0x_acnt_property_m£']=float(((-feature['Sx0_fatalacc']*fatal_property_mGBP[0])+(-feature['Sx0_seriousacc']*serious_property_mGBP[0])+(-feature['Sx0_slightacc']*slight_property_mGBP[0]))*standard_disc) feature['S0x_acnt_insurance_m£']=float(((-feature['Sx0_fatalacc']*fatal_insurance_mGBP[0])+(-feature['Sx0_seriousacc']*serious_insurance_mGBP[0])+(-feature['Sx0_slightacc']*slight_insurance_mGBP[0]))*standard_disc) #User and provider impacts #Trip distribution by km travelled if 0.34 <= feature['PopDensmean']<=1.51: #Rural village and dispersed #Calculate distribution ratios feature['vehT0cycl']=float((rvillage_cycl_tripsperson[0]*feature['Pop.sum'])) feature['vehT0passcar']=float((rvillage_passcar_tripsperson[0]*feature['Pop.sum'])/(passcar_occup[0])) feature['vehT0motorcycl']=float((rvillage_motorcycle_tripsperson[0]*feature['Pop.sum'])) feature['vehT0bus']=float((rvillage_bus_tripsperson[0]*feature['Pop.sum'])/(bus_occup[0])) feature['vehT0taxi']=float((rvillage_taxi_tripsperson[0]*feature['Pop.sum'])/(taxi_occup[0])) #Calculate Sx trips feature['vehTxcycl']=float(feature['vehT0cycl']*((T0ms_percent[8]+rural_mschange_percent[8])/T0ms_percent[8])) feature['vehTxpasscar']=float(feature['vehT0passcar']*((T0ms_percent[0]+rural_mschange_percent[0])/T0ms_percent[0])) feature['vehTxmotorcycl']=float(feature['vehT0motorcycl']*((T0ms_percent[6]+rural_mschange_percent[6])/T0ms_percent[6])) feature['vehTxbus']=float(feature['vehT0bus']*((T0ms_percent[7]+rural_mschange_percent[7])/T0ms_percent[7])) feature['vehTxtaxi']=float(feature['vehT0taxi']*((T0ms_percent[9]+rural_mschange_percent[9])/T0ms_percent[9])) if 1.51 < feature['PopDensmean']<=7.19: #Rural towns #Calculate distribution ratios feature['vehT0cycl']=float((rtown_cycl_tripsperson[0]*feature['Pop.sum'])) feature['vehT0passcar']=float((rtown_passcar_tripsperson[0]*feature['Pop.sum'])/(passcar_occup[0])) feature['vehT0motorcycl']=float((rtown_motorcycle_tripsperson[0]*feature['Pop.sum'])) feature['vehT0bus']=float((rtown_bus_tripsperson[0]*feature['Pop.sum'])/(bus_occup[0])) feature['vehT0taxi']=float((rtown_taxi_tripsperson[0]*feature['Pop.sum'])/(taxi_occup[0])) #Calculate Sx trips feature['vehTxcycl']=float(feature['vehT0cycl']*((T0ms_percent[8]+rural_mschange_percent[8])/T0ms_percent[8])) feature['vehTxpasscar']=float(feature['vehT0passcar']*((T0ms_percent[0]+rural_mschange_percent[0])/T0ms_percent[0])) feature['vehTxmotorcycl']=float(feature['vehT0motorcycl']*((T0ms_percent[6]+rural_mschange_percent[6])/T0ms_percent[6])) feature['vehTxbus']=float(feature['vehT0bus']*((T0ms_percent[7]+rural_mschange_percent[7])/T0ms_percent[7])) feature['vehTxtaxi']=float(feature['vehT0taxi']*((T0ms_percent[9]+rural_mschange_percent[9])/T0ms_percent[9])) elif feature['PopDensmean']>7.19: #Urban town #Calculate distribution ratios feature['vehT0cycl']=float((urban_cycl_tripsperson[0]*feature['Pop.sum'])) feature['vehT0passcar']=float((urban_passcar_tripsperson[0]*feature['Pop.sum'])/(passcar_occup[0])) feature['vehT0motorcycl']=float((urban_motorcycle_tripsperson[0]*feature['Pop.sum'])) feature['vehT0bus']=float((urban_bus_tripsperson[0]*feature['Pop.sum'])/(bus_occup[0])) feature['vehT0taxi']=float((urban_taxi_tripsperson[0]*feature['Pop.sum'])/(taxi_occup[0])) #Calculate Sx trips feature['vehTxcycl']=float(feature['vehT0cycl']*((T0ms_percent[8]+urban_mschange_percent[8])/T0ms_percent[8])) feature['vehTxpasscar']=float(feature['vehT0passcar']*((T0ms_percent[0]+urban_mschange_percent[0])/T0ms_percent[0])) feature['vehTxmotorcycl']=float(feature['vehT0motorcycl']*((T0ms_percent[6]+urban_mschange_percent[6])/T0ms_percent[6])) feature['vehTxbus']=float(feature['vehT0bus']*((T0ms_percent[7]+urban_mschange_percent[7])/T0ms_percent[7])) feature['vehTxtaxi']=float(feature['vehT0taxi']*((T0ms_percent[9]+urban_mschange_percent[9])/T0ms_percent[9])) #Calculations for economic impacts + distribution work_usercharges=float((1/2)*(1+VAT[0])*((((feature['vehT0passcar']*wpercent_passcar[0])+(feature['vehTxpasscar']*wpercent_passcar[0])))*(motor_chargeS0[0]-motor_chargeSx[0]))) non_work_usercharges=float((1/2)*((((feature['vehT0passcar']*(1-wpercent_passcar[0]))+(feature['vehTxpasscar']*(1-wpercent_passcar[0]))))*(motor_chargeS0[0]-motor_chargeSx[0]))) feature['user_charges_m£']=float(((work_usercharges+non_work_usercharges)*(10**-6))*standard_disc) #Calculate PT provider revenues feature['PTprovider_revenue_m£']=float((((bus_ticket_Sx[0]*feature['vehTxbus'])-(bus_ticket_S0[0]*feature['vehT0bus'])) +((taxi_cost_Sx[0]*feature['vehTxtaxi'])-(taxi_cost_S0[0]*feature['vehT0taxi'])))*(10**-6)*standard_disc) #Calculating travel time work_traveltime=float((1/2)*(1+VAT[0])*(((wpercent_cycling[0]*(feature['vehT0cycl']+feature['vehTxcycl']))*(workh_cycling[0]*(S0_cycling_wGBP[0]-Sx_cycling_wGBP[0]))) +((wpercent_motorcycle[0]*(feature['vehT0motorcycl']+feature['vehTxmotorcycl']))*(workh_passcar[0]*(S0_motorcycle_wGBP[0]-Sx_motorcycle_wGBP[0]))) +((wpercent_passcar[0]*(feature['vehT0passcar']+feature['vehTxpasscar']))*(workh_passcar[0]*(S0_passcar_wGBP[0]-Sx_passcar_wGBP[0]))) +((wpercent_bus[0]*(feature['vehT0bus']+feature['vehTxbus']))*(workh_PT[0]*(S0_bus_wGBP[0]-Sx_bus_wGBP[0]))) +((wpercent_taxi[0]*(feature['vehT0taxi']+feature['vehTxtaxi']))*(workh_PT[0]*(S0_taxi_wGBP[0]-Sx_taxi_wGBP[0]))))) nwork_traveltime=float((1/2)*((((1-wpercent_cycling[0])*(feature['vehT0cycl']+feature['vehTxcycl']))*(nworkh_cycling[0]*(S0_nwGBP[0]-Sx_nwGBP[0]))) +(((1-wpercent_motorcycle[0])*(feature['vehT0motorcycl']+feature['vehTxmotorcycl']))*(nworkh_passcar[0]*(S0_nwGBP[0]-Sx_nwGBP[0]))) +(((1-wpercent_passcar[0])*(feature['vehT0passcar']+feature['vehTxpasscar']))*(nworkh_passcar[0]*(S0_nwGBP[0]-Sx_nwGBP[0]))) +(((1-wpercent_bus[0])*(feature['vehT0bus']+feature['vehTxbus']))*(nworkh_PT[0]*(S0_nwGBP[0]-Sx_nwGBP[0]))) +(((1-wpercent_taxi[0])*(feature['vehT0taxi']+feature['vehTxtaxi']))*(nworkh_PT[0]*(S0_nwGBP[0]-Sx_nwGBP[0]))))) feature['travel_time_m£']=float((nwork_traveltime+work_traveltime)*(10**-6)*standard_disc) #Calculating operational costs feature['PT_operational_m£']=float((((1/2)*(1+bus_fueltx[0])*(((((feature['vehT0bus']*wpercent_bus[0])+((feature['vehTxbus']*wpercent_bus[0])*(1-reg_e[7]))))*(fop_S0w_bus[0]-fop_Sxw_bus[0])))) +((1/2)*(1+e_VAT[0])*((((feature['vehTxbus']*wpercent_bus[0])*reg_e[7]))*(fop_S0w_ebus[0]-fop_Sxw_ebus[0]))) +((1/2)*(1+car_fueltx[0])*((((feature['vehT0taxi']*wpercent_taxi[0])+((feature['vehTxtaxi']*wpercent_taxi[0])*(1-reg_e[9]))))*(fop_S0w_passcar[0]+nfop_S0wnw_passcar[0]-fop_Sxw_passcar[0]-nfop_Sxwnw_passcar[0]))) +((1/2)*(1+e_VAT[0])*((((feature['vehTxtaxi']*wpercent_taxi[0])*reg_e[9]))*(fop_S0w_epasscar[0]+nfop_S0wnw_passcar[0]-fop_Sxw_epasscar[0]-nfop_Sxwnw_passcar[0]))))*(10**-6)*standard_disc) wprivate_operational=float(((1/2)*(1+motorcycle_fueltx[0])*((((feature['vehT0motorcycl']*wpercent_motorcycle[0])+(feature['vehTxmotorcycl']*wpercent_motorcycle[0])))*(fop_S0w_passcar[0]+nfop_S0wnw_passcar[0]-fop_Sxw_passcar[0]-nfop_Sxwnw_passcar[0]))) +((1/2)*(1+car_fueltx[0])*((((feature['vehT0passcar']*wpercent_passcar[0])+((feature['vehTxpasscar']*wpercent_passcar[0])*(1-reg_e[0]))))*(fop_S0w_passcar[0]+nfop_S0wnw_passcar[0]-fop_Sxw_passcar[0]-nfop_Sxwnw_passcar[0]))) +((1/2)*(1+e_VAT[0])*((((feature['vehTxpasscar']*wpercent_passcar[0])*reg_e[0]))*(fop_S0w_epasscar[0]+nfop_S0wnw_passcar[0]-fop_Sxw_epasscar[0]-nfop_Sxwnw_passcar[0])))) nwprivate_operational=float((((1/2)*((((feature['vehT0motorcycl']*(1-wpercent_motorcycle[0]))+(feature['vehTxmotorcycl']*(1-wpercent_motorcycle[0]))))*(fop_S0nw_passcar[0]-fop_Sxnw_passcar[0])))-(((feature['vehTxmotorcycl']*(1-wpercent_motorcycle[0]))*nfop_Sxwnw_passcar[0])-((feature['vehT0motorcycl']*(1-wpercent_motorcycle[0]))*nfop_S0wnw_passcar[0]))) +(((1/2)*((((feature['vehT0passcar']*(1-wpercent_passcar[0]))+((feature['vehTxpasscar']*(1-wpercent_passcar[0]))*(1-reg_e[0]))))*(fop_S0nw_passcar[0]-fop_Sxnw_passcar[0])))-(((feature['vehTxpasscar']*(1-wpercent_passcar[0]))*(1-reg_e[0])*nfop_Sxwnw_passcar[0])-((feature['vehT0passcar']*(1-wpercent_passcar[0]))*nfop_S0wnw_passcar[0]))) +(((1/2)*((((feature['vehTxpasscar']*(1-wpercent_passcar[0]))*reg_e[0]))*(fop_S0nw_epasscar[0]-fop_Sxnw_epasscar[0])))-(((feature['vehTxpasscar']*(1-wpercent_passcar[0]))*reg_e[0]*nfop_Sxwnw_passcar[0])))) feature['private_operational_m£']=float((wprivate_operational+nwprivate_operational)*(10**-6)*standard_disc) #Calculate NPV feature['NPV_m£']=float(feature['S0x_NOx_m£']+feature['S0x_PM25_m£']+feature['S0x_CO2_m£'] +feature['S0x_cyc_m£']+feature['S0x_acnt_netout_m£']+feature['S0x_acnt_medical_m£'] +feature['S0x_acnt_WTP_m£']+feature['S0x_acnt_police_m£']+feature['S0x_acnt_property_m£'] +feature['S0x_acnt_insurance_m£']+feature['user_charges_m£']+feature['PTprovider_revenue_m£'] +feature['travel_time_m£']+feature['PT_operational_m£']+feature['private_operational_m£']) economic_impacts.updateFeature(feature) #update layer changes economic_impacts.commitChanges() #record layer changes and stop editing #Clean layer economic_impacts.startEditing() #start editing economic_impacts.dataProvider().deleteAttributes([4,5,6,7,8,9,10,11,12,13,14,15]) economic_impacts.commitChanges() #record layer changes and stop editing #Save all layers QgsVectorFileWriter.writeAsVectorFormat(Sx_AADF_NOx_grid,path_to_emissionsresults,"utf-8",Sx_AADF_NOx_grid.crs(),"GPKG") QgsVectorFileWriter.writeAsVectorFormat(hhlayer,path_to_hhresults,"utf-8",hhlayer.crs(),"GPKG") QgsVectorFileWriter.writeAsVectorFormat(economic_impacts,path_to_eiresults,"utf-8",economic_impacts.crs(),"GPKG")