Commit 825dfbcd authored by Henrik Grythe's avatar Henrik Grythe
Browse files

Have a working version of road link emission

parent edd24c91
......@@ -32,13 +32,12 @@ if ~ismember(Vehicle_source,{'SSB'})
end
fprintf('\tin Emission_Factor_Model_group_HBEFA\n')
fprintf('in Emission_Factor_Model_group_HBEFA *\n')
fprintf('--- Grouping Vehicles according to ---\n--- Model SSB HBEFA merger specifications ---\n%s \n',SSB_Vehicle_dist)
T = readtable(SSB_Vehicle_dist,'Sheet','HBEFA778toMODEL');
T = readtable(SSB_Vehicle_dist,'Sheet','HBEFA778toMODEL');
Tm = readtable(SSB_Vehicle_dist,'Sheet','MODEL');
fprintf('read Sheet : %s and Sheet %s\n','HBEFA778toMODEL','MODEL')
% ModelNumber HBEFA_Weight Uniform_Weight
modelN = unique(T.ModelNumber);
......@@ -49,7 +48,6 @@ fprintf('Found %i Model Vehicles \n',length(modelN))
% Loop all components (comps) to do conversion for
for com =1:length(comps)
fprintf('<--- \nSpec| %s\n',char(comps(com)))
iEFfile = sprintf('%s/EFA_matrix41_RAW_%s.mat',tfold,char(comps(com)));
oEFfile = sprintf('%s/EFA_Table_MODEL_%s',tfold,char(comps(com)));
......@@ -85,11 +83,11 @@ for com =1:length(comps)
fprintf('MODEL Vehicle: %3i %-42s Found #%i EF',modelN(mod),char(Tm.Name(modelN(mod))),uef)
Tout = table;
k = 1;
for roa = 1:size(EFA,2)
for spd = 1:size(EFA,3)
for dec = 1:size(EFA,4)
for cog = 1:size(EFA,5)
for urb = 1:size(EFA,6)
for roa = 1:size(EF_AVG,2)
for spd = 1:size(EF_AVG,3)
for dec = 1:size(EF_AVG,4)
for cog = 1:size(EF_AVG,5)
for urb = 1:size(EF_AVG,6)
Trout = table;
if ~isnan(EFsub(1,roa,spd,dec,cog,urb))
Trout.Name = {sprintf('%s/%s-%i/%i%%/%s',char(roads.RoadEnv(urb)),char(roads.RoadType(roa)),roads.RoadSpeeds(spd),roads.RoadGradient(dec),char(roads.Congestion(cog)))};
......@@ -165,7 +163,7 @@ for com =1:length(comps)
% Concatonate only the column of the new vehicle data if the
% align.
[a,b,c] = intersect(TFout.Name,Tout.Name);
addCol = find(ismember(Tout.Properties.VariableNames,Tm.Name(modelN(mod))));
addCol = find(ismember(Tout.Properties.VariableNames,Tm.Name(modelN(mod))));
% Not 100% robust test!!!!?
if length(b)==length(c)
......
......@@ -15,7 +15,10 @@
%--------------------------------------------------------------------------
function Emission_Factor_Test_HBEFA_to_MODEL_conversion(com,ConVfile)
global EFA EF_AVG comps roads
fprintf('\tin Emission_Factor_Test_HBEFA_to_MODEL_conversion\n')
fprintf('in Emission_Factor_Test_HBEFA_to_MODEL_conversion * \n')
Traff = table;
k = 1;
for roa = 1:size(EF_AVG,2)
......
......@@ -13,7 +13,6 @@ VD = Vehicle_dist.Vdist;
% Loop all components to be calculated:
for com = 1:length(comps)
% Load the component file with
fprintf('<--- %s\n',char(comps(com)))
TFout = readtable('modelHBEFA.xlsx','Sheet',sprintf('%s_%s_Weight',char(comps(com)),Vehicle_weight),'PreserveVariableNames',1);
......@@ -24,13 +23,18 @@ for com = 1:length(comps)
Trout = TFout(:,1:7);
for komm = 1: length(Vehicle_dist.D1_KommNr)
if debug_mode
fprintf('Sum Weights: Komm: %04i@\n \tT:%f\n \tL:%f\n \tH:%f\n \tB:%f \n',...
Vehicle_dist.D1_KommNr(komm),sum(VD(komm,:)),...
sum(VD(komm,LightVehiclesIdx)),...
sum(VD(komm,HeavyVehiclesIdx)),...
sum(VD(komm,BusesVehiclesIdx)))
else
fprintf('Komm:%04i;',Vehicle_dist.D1_KommNr(komm))
if rem(komm,15)==0;fprintf('\n'); end
end
if abs(sum(VD(komm,:))-3)>1e-5
VD(komm,LightVehiclesIdx) = VD(komm,LightVehiclesIdx)/sum(VD(komm,LightVehiclesIdx));
VD(komm,HeavyVehiclesIdx) = VD(komm,HeavyVehiclesIdx)/sum(VD(komm,HeavyVehiclesIdx));
......@@ -46,21 +50,31 @@ for com = 1:length(comps)
for i= 1:height(TFout)
EFs = VD(komm,:).*table2array(TFout(i,idef));
% rows2vars(TFout(1,idef))
EFLight(i,1) = nansum(EFs(LightVehiclesIdx));
EFHeavy(i,1) = nansum(EFs(HeavyVehiclesIdx));
EFBuses(i,1) = nansum(EFs(BusesVehiclesIdx));
end
Trout.Light = EFLight;
Trout.Heavy = EFHeavy;
Trout.Buses = EFBuses;
Trout.Buses = EFBuses;
Trout.Properties.VariableNames(find(ismember(Trout.Properties.VariableNames,'Light'))) = {sprintf('EF_Light_%04i',Vehicle_dist.D1_KommNr(komm))};
Trout.Properties.VariableNames(find(ismember(Trout.Properties.VariableNames,'Heavy'))) = {sprintf('EF_Heavy_%04i',Vehicle_dist.D1_KommNr(komm))};
Trout.Properties.VariableNames(find(ismember(Trout.Properties.VariableNames,'Buses'))) = {sprintf('EF_Buses_%04i',Vehicle_dist.D1_KommNr(komm))};
end
writetable(Trout,'OnRoadEF_RoadClasses.xlsx','Sheet',sprintf('%s_%i',char(comps(com)),Tyear))
switch char(comps(com))
case 'CO2'
fprintf('### ADD CALL TO: Emission_Factor_mix_in_biofuels()\n')
end
ofile = 'OnRoadEF_RoadClasses.xlsx';
OnRoadEF_RoadClasses = Trout;
save('OnRoadEF_RoadClasses.mat','OnRoadEF_RoadClasses')
writetable(Trout,ofile,'Sheet',sprintf('%s_%i',char(comps(com)),Tyear))
fprintf('%s--- >\n',char(comps(com)))
end
end
\ No newline at end of file
function Emission_Factors_Road_DrivingDistance_IN_Municipalities()
global RLinks tfold Tyear Vehicle_dist debug_mode SSB_Vehicle_dist
TM = readtable(SSB_Vehicle_dist,'Sheet','MODEL');
LightVehiclesIdx = TM.ClassNum==1|TM.ClassNum==2;
BusesVehiclesIdx = TM.ClassNum==3|TM.ClassNum==4;
HeavyVehiclesIdx = TM.ClassNum==5|TM.ClassNum==6|TM.ClassNum==7;
T = readtable('modelHBEFA.xlsx','Sheet','Roads');
file = sprintf('%s%s',tfold,'roads');
fprintf('### Warning,using roads file not produced by NERVE model\n%s\n',file)
load(file)
% All road links kilometer per type of road Light / Heavy / Bus
nDD_L = zeros(height(T),1);
nDD_H = zeros(height(T),1);
nDD_B = zeros(height(T),1);
%--------------------------------------------------------------
% ROAD LINK Calculations
%--------------------------------------------------------------
% extract fields needed for calculations
fprintf('\nExtracts necessary fields \n')
fprintf('Traffic and Vehicle Year: %i \n',Tyear)
fprintf('RoadLinks : %i \n',size(RLinks,1))
dayInYear = round(datenum([Tyear+1, 1, 1, 12 0 0 ])-datenum([Tyear, 1, 1, 12, 0 ,0]));
fprintf('Found %3i days in year : %i \n',dayInYear,Tyear)
KommunerS = extractfield(RLinks,'KOMMS');
KommunerE = extractfield(RLinks,'KOMME');
Kommuner = extractfield(RLinks,'KOMM');
uKomm = unique(KommunerS);
% Congestion weights. Traffic part delayed (TPD)
% Combine the different congestion levels by volume to determine the
% actual EF:
fprintf('### Warning, using parameterized CONGESTION\n%s\n',file)
TPD_L=[1.0, 0.0, 0.0, 0.0, 0.0;...
0.6, 0.4, 0.0, 0.0, 0.0;...
0.6, 0.2, 0.2, 0.0, 0.0;...
0.5, 0.2, 0.15, 0.15, 0.0];
TPD_H=[1.0, 0.0, 0.0, 0.0, 0.0;...
0.6, 0.4, 0.0, 0.0, 0.0;...
0.6, 0.2, 0.2, 0.0, 0.0;...
0.5, 0.2, 0.15, 0.15, 0.0];
TPD_B=[1.0, 0.0, 0.0, 0.0, 0.0;...
0.6, 0.4, 0.0, 0.0, 0.0;...
0.6, 0.2, 0.2, 0.0, 0.0;...
0.5, 0.2, 0.15, 0.15, 0.0];
%-----------------------------------
TPD =[1.0, 0.0, 0.0, 0.0, 0.0;...
0.6, 0.4, 0.0, 0.0, 0.0;...
0.6, 0.2, 0.2, 0.0, 0.0;...
0.5, 0.2, 0.15, 0.15, 0.0];
Tout = T;
for komm = 1:length(uKomm)
% extract all necessary roadLinks:
f1 = find(KommunerS == uKomm(komm));
f2 = find(KommunerE == uKomm(komm));
f3 = find(Kommuner == uKomm(komm));
f = unique([f1,f2,f3]);
MRLinks = RLinks(f);
fprintf('\n-%3i Kommune_%04i --- \n',komm,Vehicle_dist.D1_KommNr(komm))
fprintf(' Has %i (%i) Roads \n',length(f1),length(f))
L = extractfield(MRLinks,sprintf('L_ADT%04i',Tyear)); % Traffic Volume (# day-1)
H = extractfield(MRLinks,sprintf('H_ADT%04i',Tyear)); % Traffic Volume (# day-1)
B = extractfield(MRLinks,sprintf('B_ADT%04i',Tyear)); % Traffic Volume (# day-1)
HBEFA = extractfield(MRLinks,'HBEFA_EQIV'); % Type of Road (Access->MW)
SPD = extractfield(MRLinks,'SPEED'); % Speed on Road (km hr-1)
DEC = extractfield(MRLinks,'SLOPE'); % Verticality of Road (%)
ENV = extractfield(MRLinks,'URBAN'); % Urbanity of Road (URB/RUR)
LEN = extractfield(MRLinks,'DISTANCE'); % Length of Road (in kilometres)
RU = extractfield(MRLinks,'RUSH_DELAY');
% Classification of rush hour traffic. 0-4 (4 is congested)
ru = zeros(size(RU));
ru(RU<=0.1) = 1;
ru(RU>0.1 & RU<=4) = 2;
ru(RU>4 & RU<=15) = 3;
ru(RU>15) = 4;
IDO = extractfield(MRLinks,'IDO');
KOMS = extractfield(MRLinks,'KOMMS');
if debug_mode
fprintf(' Light Traffic L=%7.1f (1 000 000) Km TDL=%9.3f \n',1e-6*sum(L.*IDO.*LEN)*dayInYear,1e-6*sum(Vehicle_dist.modelTD(komm,LightVehiclesIdx)))
fprintf(' Heavy Traffic H=%7.1f (1 000 000) Km TDH=%9.3f \n',1e-6*sum(H.*IDO.*LEN)*dayInYear,1e-6*sum(Vehicle_dist.modelTD(komm,HeavyVehiclesIdx)))
fprintf(' Buses Traffic B=%7.1f (1 000 000) Km TDB=%9.3f \n',1e-6*sum(B.*IDO.*LEN)*dayInYear,1e-6*sum(Vehicle_dist.modelTD(komm,BusesVehiclesIdx)))
end
% *my is the total distance driven on each ROAD link Type in Tyear
Lmy = zeros(size(T,1),1);
Hmy = zeros(size(T,1),1);
Bmy = zeros(size(T,1),1);
for r =1:length(L)
dec = find(roads.RoadGradient == DEC(r));
env = find(roads.RoadEnvID == ENV(r)+1);
SPD(r) = round(SPD(r)/10)*10;
if SPD(r) < 30; SPD(r)=30; end
if SPD(r) > 110; SPD(r)=110; end
spd = find(roads.RoadSpeeds == SPD(r));
name = {sprintf('%s/%s-%i/%i%%/%s',char(roads.RoadEnv(env)),char(roads.RoadType(HBEFA(r))),roads.RoadSpeeds(spd),roads.RoadGradient(dec))};
idx = find(contains(T.Name,name));
if KOMS(r) == uKomm(komm)
Lmy(idx) = Lmy(idx) + dayInYear *IDO(r)* L(r) *LEN(r)*TPD_L(ru(r),:)';
Hmy(idx) = Hmy(idx) + dayInYear *IDO(r)* H(r) *LEN(r)*TPD_H(ru(r),:)';
Bmy(idx) = Bmy(idx) + dayInYear *IDO(r)* B(r) *LEN(r)*TPD_B(ru(r),:)';
else
Lmy(idx) = Lmy(idx) + dayInYear *IDO(r)* L(r) *LEN(r)*TPD_L(ru(r),:)';
Hmy(idx) = Hmy(idx) + dayInYear *IDO(r)* H(r) *LEN(r)*TPD_H(ru(r),:)';
Bmy(idx) = Bmy(idx) + dayInYear *IDO(r)* B(r) *LEN(r)*TPD_B(ru(r),:)';
end
end
if abs((sum(Lmy)/(sum(L.*IDO.*LEN)*dayInYear))-1)>1e-3
fprintf('### more than 1%% shift %f \n',sum(Lmy)/(sum(L.*IDO.*LEN)*dayInYear))
end
Tout.Light_DD = Lmy;
Tout.Heavy_DD = Hmy;
Tout.Buses_DD = Bmy;
Tout.Properties.VariableNames(find(ismember(Tout.Properties.VariableNames,'Light_DD'))) = {sprintf('DD_Light_%04i',Vehicle_dist.D1_KommNr(komm))};
Tout.Properties.VariableNames(find(ismember(Tout.Properties.VariableNames,'Heavy_DD'))) = {sprintf('DD_Heavy_%04i',Vehicle_dist.D1_KommNr(komm))};
Tout.Properties.VariableNames(find(ismember(Tout.Properties.VariableNames,'Buses_DD'))) = {sprintf('DD_Buses_%04i',Vehicle_dist.D1_KommNr(komm))};
end
save(sprintf('Municipal_DrivingDistances_per_RoadType_%i.mat',Tyear),'Tout')
writetable(Tout,'Municipal_DrivingDistances_per_RoadType.xlsx','Sheet',sprintf('DD_%i',Tyear))
fprintf('\n---- NORGE --- \n')
L = extractfield(RLinks,sprintf('L_ADT%04i',Tyear)); % Traffic Volume (# day-1)
H = extractfield(RLinks,sprintf('H_ADT%04i',Tyear)); % Traffic Volume (# day-1)
B = extractfield(RLinks,sprintf('B_ADT%04i',Tyear)); % Traffic Volume (# day-1)
LEN = extractfield(RLinks,'DISTANCE');
LW = 1e-6*sum(L.*LEN)*dayInYear;
HW = 1e-6*sum(H.*LEN)*dayInYear;
BW = 1e-6*sum(B.*LEN)*dayInYear;
LTD = 1e-6*sum(sum(Vehicle_dist.modelTD(:,LightVehiclesIdx)));
HTD = 1e-6*sum(sum(Vehicle_dist.modelTD(:,HeavyVehiclesIdx)));
BTD = 1e-6*sum(sum(Vehicle_dist.modelTD(:,BusesVehiclesIdx)));
fprintf(' Light Traffic L=%7.1f (1 000 000) Km TDL=%7.1f (%5.1f%%) \n',LW,LTD,100*LW/LTD)
fprintf(' Heavy Traffic H=%7.1f (1 000 000) Km TDH=%7.1f (%5.1f%%) \n',HW,HTD,100*HW/HTD)
fprintf(' Buses Traffic B=%7.1f (1 000 000) Km TDB=%7.1f (%5.1f%%) \n',BW,BTD,100*BW/BTD)
end
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
% % EMISS_* are the road link emissions for each link in the
% % municipality
% EMISS_L = zeros(size(MRLinks));
% EMISS_H = zeros(size(MRLinks));
% EMISS_B = zeros(size(MRLinks));
%
% for r =1:length(L)
% dec = find(roads.RoadGradient == DEC(r));
% env = find(roads.RoadEnvID == ENV(r)+1);
% SPD(r) = round(SPD(r)/10)*10;
% if SPD(r) < 30; SPD(r)=30; end
% if SPD(r) > 110; SPD(r)=110; end
% spd = find(roads.RoadSpeeds == SPD(r));
% name = {sprintf('%s/%s-%i/%i%%/%s',char(roads.RoadEnv(env)),char(roads.RoadType(HBEFA(r))),roads.RoadSpeeds(spd),roads.RoadGradient(dec))};
%
% idx = find(contains(Tef.Name,name));
% if KOMS(r) == uKomm(komm)
% Lmy(idx) = Lmy(idx) + dayInYear *IDO(r)* L(r) *LEN(r)*TPD_L(ru(r),:)';
% Hmy(idx) = Hmy(idx) + dayInYear *IDO(r)* H(r) *LEN(r)*TPD_H(ru(r),:)';
% Bmy(idx) = Bmy(idx) + dayInYear *IDO(r)* B(r) *LEN(r)*TPD_B(ru(r),:)';
% else
% Lmy(idx) = Lmy(idx) + dayInYear *IDO(r)* L(r) *LEN(r)*TPD_L(ru(r),:)';
% Hmy(idx) = Hmy(idx) + dayInYear *IDO(r)* H(r) *LEN(r)*TPD_H(ru(r),:)';
% Bmy(idx) = Bmy(idx) + dayInYear *IDO(r)* B(r) *LEN(r)*TPD_B(ru(r),:)';
% end
% % EMISS_L(r) = sum(Lmy(idx).*table2array(Tef(idx,Lef)))/sum(Lmy(idx));
% % EMISS_H(r) = sum(Hmy(idx).*table2array(Tef(idx,Hef)))/sum(Lmy(idx));
% % EMISS_B(r) = sum(Bmy(idx).*table2array(Tef(idx,Bef)))/sum(Lmy(idx));
% EMISS_L(r) = sum(Lmy(idx).*table2array(Tef(idx,Lef)));
% EMISS_H(r) = sum(Hmy(idx).*table2array(Tef(idx,Hef)));
% EMISS_B(r) = sum(Bmy(idx).*table2array(Tef(idx,Bef)));
% end
%
%
% tEML = Lmy.*table2array(Tef(:,Lef));
% tEMH = Hmy.*table2array(Tef(:,Hef));
% tEMB = Bmy.*table2array(Tef(:,Bef));
%
%
% % nDD_L = nDD_L + Lmy
% % nDD_H = nDD_H + Hmy
% % nDD_B = nDD_B + Bmy
%
% nEMISS_L(f) = nEMISS_L(f) + EMISS_L;
% nEMISS_H(f) = nEMISS_H(f) + EMISS_H;
% nEMISS_B(f) = nEMISS_B(f) + EMISS_B;
%
%
% if debug_mode
% fprintf('---- Lette %11.1f Kilometer (%3.0f%%) \n' ,sum(L),100*sum(Lmy)/sum(Lmy+Hmy+Bmy))
% fprintf('---- Tunge %11.1f Kilometer (%3.0f%%) \n' ,sum(H),100*sum(Hmy)/sum(Lmy+Hmy+Bmy))
% fprintf('---- Busser %11.1f Kilometer (%3.0f%%) \n\n',sum(B),100*sum(Bmy)/sum(Lmy+Hmy+Bmy))
%
% fprintf('---- Lette %11.1f Gram %s (%3.0f%%)\n' ,sum(tEML),char(comps(com)),100*sum(tEML)/sum(tEML+tEMH+tEMB))
% fprintf('---- Tunge %11.1f Gram %s (%3.0f%%)\n' ,sum(tEMH),char(comps(com)),100*sum(tEMH)/sum(tEML+tEMH+tEMB))
% fprintf('---- Busser %11.1f Gram %s (%3.0f%%)\n\n',sum(tEMB),char(comps(com)),100*sum(tEMB)/sum(tEML+tEMH+tEMB))
%
% fprintf('---- Lette %11.1f Gram/Kilometer \n' ,sum(tEML)/sum(Lmy))
% fprintf('---- Tunge %11.1f Gram/Kilometer \n' ,sum(tEMH)/sum(Hmy))
% fprintf('---- Busser %11.1f Gram/Kilometer \n\n',sum(tEMB)/sum(Bmy))
% end
%
%
%
% KEF(komm,1) = sum(Tef(:,Lef).*Tef.Light_DD)/sum(Tef.Light_DD)
%
%
% Tef.Properties.VariableNames(find(ismember(Tef.Properties.VariableNames,'Light_DD'))) = {sprintf('DD_Light_%04i',Vehicle_dist.D1_KommNr(komm))};
% Tef.Properties.VariableNames(find(ismember(Tef.Properties.VariableNames,'Heavy_DD'))) = {sprintf('DD_Heavy_%04i',Vehicle_dist.D1_KommNr(komm))};
% Tef.Properties.VariableNames(find(ismember(Tef.Properties.VariableNames,'Buses_DD'))) = {sprintf('DD_Buses_%04i',Vehicle_dist.D1_KommNr(komm))};
%
%
%
% end
%
%
% fprintf('\n---- NORGE --- \n')
% fprintf('---- Lette %11.1f (1000)Ton %s (%3.0f%%)\n' ,1e-9*sum(nEMISS_L),char(comps(com)),100*nansum(nEMISS_L)/nansum(nEMISS_L+nEMISS_H+nEMISS_B))
% fprintf('---- Tunge %11.1f (1000)Ton %s (%3.0f%%)\n' ,1e-9*nansum(nEMISS_H),char(comps(com)),100*nansum(nEMISS_H)/nansum(nEMISS_L+nEMISS_H+nEMISS_B))
% fprintf('---- Busser %11.1f (1000)Ton %s (%3.0f%%)\n' ,1e-9*nansum(nEMISS_B),char(comps(com)),100*nansum(nEMISS_B)/nansum(nEMISS_L+nEMISS_H+nEMISS_B))
% fprintf('---- Totalt %11.1f (1000)Ton %s \n\n',1e-9*nansum(nEMISS_B+nEMISS_H+nEMISS_L),char(comps(com)))
%
% L = extractfield(RLinks,sprintf('L_ADT%04i',Tyear)); % Traffic Volume (# day-1)
% H = extractfield(RLinks,sprintf('H_ADT%04i',Tyear)); % Traffic Volume (# day-1)
% B = extractfield(RLinks,sprintf('B_ADT%04i',Tyear)); % Traffic Volume (# day-1)
% LEN = extractfield(RLinks,'DISTANCE');
% LW = 1e-6*sum(L.*LEN)*dayInYear;
% HW = 1e-6*sum(H.*LEN)*dayInYear;
% BW = 1e-6*sum(B.*LEN)*dayInYear;
% LTD = 1e-6*sum(sum(Vehicle_dist.modelTD(:,LightVehiclesIdx)));
% HTD = 1e-6*sum(sum(Vehicle_dist.modelTD(:,HeavyVehiclesIdx)));
% BTD = 1e-6*sum(sum(Vehicle_dist.modelTD(:,BusesVehiclesIdx)));
%
% fprintf(' Light Traffic L=%7.1f (1 000 000) Km TDL=%7.1f (%5.1f%%) \n',LW,LTD,100*LW/LTD)
% fprintf(' Heavy Traffic H=%7.1f (1 000 000) Km TDH=%7.1f (%5.1f%%) \n',HW,HTD,100*HW/HTD)
% fprintf(' Buses Traffic B=%7.1f (1 000 000) Km TDB=%7.1f (%5.1f%%) \n',BW,BTD,100*BW/BTD)
% fprintf(' Light Traffic L=%7.1f g/Km\n',sum(nEMISS_L)/LW)
% fprintf(' Heavy Traffic H=%7.1f g/Km\n',sum(nEMISS_H)/HW)
% fprintf(' Buses Traffic B=%7.1f g/Km\n',sum(nEMISS_B)/BW)
%
%
%
%
%
% TLinks = struct2table(RLinks);
% TLinks.LetEM = nEMISS_L;
% TLinks.HeaEM = nEMISS_H;
% TLinks.BusEM = nEMISS_B;
%
% TLinks.Properties.VariableNames(find(ismember(TLinks.Properties.VariableNames,'LetEM'))) = {sprintf('EM_Light_%s',char(comps(com)))};
% TLinks.Properties.VariableNames(find(ismember(TLinks.Properties.VariableNames,'HeaEM'))) = {sprintf('EM_Heavy_%s',char(comps(com)))};
% TLinks.Properties.VariableNames(find(ismember(TLinks.Properties.VariableNames,'BusEM'))) = {sprintf('EM_Buses_%s',char(comps(com)))};
% RLinks = table2struct(TLinks);
% end
......@@ -30,6 +30,7 @@ end
if ismember(Vehicle_source,{'SSB'})
Emission_Factors_OnRoadAllCond()
Emission_Factors_Road_DrivingDistance_IN_Municipalities()
[Sn] = Emissions_Calculations_SSB();
return
end
......
......@@ -22,14 +22,10 @@ function [Sn] = Emissions_Calculations_SSB()
global tfold Tyear SSB_Vehicle_dist comps RLinks Vehicle_dist Vehicle_weight
global debug_mode
%--------------------------------------------------------------
% MUNICIPALITY EMISSION FACTORS FOR ALL ROADS
%--------------------------------------------------------------
%--------------------------------------------------------------
% ROAD LINK
% ROAD LINK Calculations
%--------------------------------------------------------------
% extract fields needed for calculations
fprintf('\nExtracts necessary fields \n')
......@@ -44,22 +40,11 @@ TM = readtable(SSB_Vehicle_dist,'Sheet','MODEL');
LightVehiclesIdx = TM.ClassNum==1|TM.ClassNum==2;
BusesVehiclesIdx = TM.ClassNum==3|TM.ClassNum==4;
HeavyVehiclesIdx = TM.ClassNum==5|TM.ClassNum==6|TM.ClassNum==7;
% HBEFA = extractfield(RLinks,'HBEFA_EQIV'); % Type of Road (Access->MW)
% SPD = extractfield(RLinks,'SPEED'); % Speed on Road (km hr-1)
% DEC = extractfield(RLinks,'SLOPE'); % Verticality of Road (%)
% ENV = extractfield(RLinks,'URBAN'); % Urbanity of Road (URB/RUR)
% LEN = extractfield(RLinks,'DISTANCE')*1000; % Length of Road (in metres)
% RU = extractfield(RLinks,'RUSH_DELAY');
dayInYear = round(datenum([Tyear+1, 1, 1, 12 0 0 ])-datenum([Tyear, 1, 1, 12, 0 ,0]));
fprintf('Found %i days in year: %i\n',dayInYear,Tyear)
% Make a cleaned version for writing to file
Tout = struct2table(RLinks);
% Tout = struct2table(RLinks);
% roads
......@@ -92,145 +77,102 @@ TPD =[1.0, 0.0, 0.0, 0.0, 0.0;...
for com = 1:length(comps)
fprintf('<--- %s ---\n',char(comps(com)))
L = extractfield(RLinks,sprintf('L_ADT%04i',Tyear)); % Traffic Volume (# day-1)
H = extractfield(RLinks,sprintf('H_ADT%04i',Tyear)); % Traffic Volume (# day-1)
B = extractfield(RLinks,sprintf('B_ADT%04i',Tyear)); % Traffic Volume (# day-1)
HBEFA = extractfield(RLinks,'HBEFA_EQIV'); % Type of Road (Access->MW)
SPD = extractfield(RLinks,'SPEED'); % Speed on Road (km hr-1)
DEC = extractfield(RLinks,'SLOPE'); % Verticality of Road (%)
ENV = extractfield(RLinks,'URBAN'); % Urbanity of Road (URB/RUR)
LEN = extractfield(RLinks,'DISTANCE'); % Length of Road (in kilometres)
RU = extractfield(RLinks,'RUSH_DELAY');
TEF = readtable('OnRoadEF_RoadClasses.xlsx','Sheet',sprintf('%s_%i',char(comps(com)),Tyear),'PreserveVariableNames',1);
% load(oEFfile,'TFout','roads');
ttef = TEF(:,1:7);
% All road links emissions Light / Heavy / Bus
nEMISS_L = zeros(size(RLinks));
nEMISS_H = zeros(size(RLinks));
nEMISS_B = zeros(size(RLinks));
% Classification of rush hour traffic. 0-4 (4 is congested)
ru = zeros(size(RU));
ru(RU<=0.1) = 1;
ru(RU>0.1 & RU<=4) = 2;
ru(RU>4 & RU<=15) = 3;
ru(RU>15) = 4;
% All road links kilometer per type of road Light / Heavy / Bus
nDD_L = zeros(height(ttef));
nDD_H = zeros(height(ttef));
nDD_B = zeros(height(ttef));
for komm = 1:length(uKomm)
% extract all necessary roadLinks:
f1 = find(KommunerS == uKomm(komm));
f2 = find(KommunerE == uKomm(komm));
f = unique([f1,f2]);
MRLinks = RLinks(f);
fprintf('\n-%3i Kommune_%04i --- \n',komm,Vehicle_dist.D1_KommNr(komm))
fprintf(' Has %i Roads \n',length(f))
% extract all roadEFs:
ref = find(contains(TEF.Properties.VariableNames,{sprintf('%i',uKomm(komm))}));
Tef = [ttef,TEF(:,ref)];
Lef = find(ismember(Tef.Properties.VariableNames,sprintf('EF_Light_%04i',Vehicle_dist.D1_KommNr(komm))));
Hef = find(ismember(Tef.Properties.VariableNames,sprintf('EF_Heavy_%04i',Vehicle_dist.D1_KommNr(komm))));
Bef = find(ismember(Tef.Properties.VariableNames,sprintf('EF_Buses_%04i',Vehicle_dist.D1_KommNr(komm))));
IDO = extractfield(RLinks,'IDO');
KOMS = extractfield(RLinks,'KOMMS');
L = extractfield(MRLinks,sprintf('L_ADT%04i',Tyear)); % Traffic Volume (# day-1)
H = extractfield(MRLinks,sprintf('H_ADT%04i',Tyear)); % Traffic Volume (# day-1)
B = extractfield(MRLinks,sprintf('B_ADT%04i',Tyear)); % Traffic Volume (# day-1)
HBEFA = extractfield(MRLinks,'HBEFA_EQIV'); % Type of Road (Access->MW)
SPD = extractfield(MRLinks,'SPEED'); % Speed on Road (km hr-1)
DEC = extractfield(MRLinks,'SLOPE'); % Verticality of Road (%)
ENV = extractfield(MRLinks,'URBAN'); % Urbanity of Road (URB/RUR)
LEN = extractfield(MRLinks,'DISTANCE'); % Length of Road (in kilometres)
RU = extractfield(MRLinks,'RUSH_DELAY');
% Classification of rush hour traffic. 0-4 (4 is congested)
ru = zeros(size(RU));
ru(RU<=0.1) = 1;
ru(RU>0.1 & RU<=4) = 2;
ru(RU>4 & RU<=15) = 3;
ru(RU>15) = 4;
IDO = extractfield(MRLinks,'IDO');
KOMS = extractfield(MRLinks,'KOMMS');
if debug_mode
fprintf(' Light Traffic L=%7.1f (1 000 000) Km TDL=%7.1f \n',1e-6*sum(L.*IDO.*LEN)*dayInYear,1e-6*sum(Vehicle_dist.modelTD(komm,LightVehiclesIdx)))
fprintf(' Heavy Traffic H=%7.1f (1 000 000) Km TDH=%7.1f \n',1e-6*sum(H.*IDO.*LEN)*dayInYear,1e-6*sum(Vehicle_dist.modelTD(komm,HeavyVehiclesIdx)))
fprintf(' Buses Traffic B=%7.1f (1 000 000) Km TDB=%7.1f \n',1e-6*sum(B.*IDO.*LEN)*dayInYear,1e-6*sum(Vehicle_dist.modelTD(komm,BusesVehiclesIdx)))
end
% All road links emissions Light / Heavy / Bus
EMISS_L = zeros(size(RLinks));
EMISS_H = zeros(size(RLinks));
EMISS_B = zeros(size(RLinks));
Link_emission_factor= zeros(size(RLinks));
for com = 1:length(comps)
fprintf('<--- %s ---\n',char(comps(com)))
fprintf('Loading large file\n...')
TEF = readtable('OnRoadEF_RoadClasses.xlsx','Sheet',sprintf('%s_%i',char(comps(com)),Tyear),'PreserveVariableNames',1);
fprintf('Loaded.\n')
tef = TEF.Name;
for r =1:length(L)
% Find the municipality to use Emission factor from:
% Find the correct emission factor column in TEF:
komm = extractfield(RLinks(r),'KOMMS');
Lef = find(ismember(TEF.Properties.VariableNames,sprintf('EF_Light_%04i',komm)));
Hef = find(ismember(TEF.Properties.VariableNames,sprintf('EF_Heavy_%04i',komm)));
Bef = find(ismember(TEF.Properties.VariableNames,sprintf('EF_Buses_%04i',komm)));
% *my is the total distance driven on each link in Tyear
Lmy = zeros(size(Tef,1),1);
Hmy = zeros(size(Tef,1),1);
Bmy = zeros(size(Tef,1),1);
% Compose the name type of the road:
dec = find(roads.RoadGradient == DEC(r));
env = find(roads.RoadEnvID == ENV(r)+1);
SPD(r) = round(SPD(r)/10)*10;
if SPD(r) < 30; SPD(r)=30; end
if SPD(r) > 110; SPD(r)=110; end
spd = find(roads.RoadSpeeds == SPD(r));