Emissions_Calculations_SSB.m 8.54 KB
Newer Older
1
2
%--------------------------------------------------------------------------
% This file is part of NERVE
3
%
4
5
6
% NERVE is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation version 3.
7
%
8
9
10
11
% NERVE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
12
%
13
14
15
% You should have received a copy of the GNU General Public License
% along with NERVE.  If not, see <https://www.gnu.org/licenses/>.
%--------------------------------------------------------------------------
16
function [Sn] = Emissions_Calculations_SSB()
17
18
19
20
21
%--------------------------------------------------------------------------
%
% 22.10.2020 -Henrik Grythe
% Kjeller NILU
%--------------------------------------------------------------------------
22
global tfold Tyear SSB_Vehicle_dist comps RLinks Vehicle_dist Vehicle_weight
Henrik Grythe's avatar
Henrik Grythe committed
23
global debug_mode ofiles
24
25

%--------------------------------------------------------------
Henrik Grythe's avatar
Henrik Grythe committed
26
% ROAD LINK EMISSIONS Calculations
27
28
29
30
%--------------------------------------------------------------
% extract fields needed for calculations
fprintf('\nExtracts necessary fields   \n')
fprintf('Traffic and Vehicle Year:  %i \n',Tyear)
31
fprintf('RoadLinks               :  %i \n',size(RLinks,1))
32
33

file = sprintf('%s%s',tfold,'roads');
Henrik Grythe's avatar
Henrik Grythe committed
34
fprintf('### Warning, using roads file not produced by NERVE model\n%s\n',file)
35
36
load(file)

37
38
39
40
41
42
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;
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)
43
44


45
46
47
% roads
KommunerS = extractfield(RLinks,'KOMMS');
KommunerE = extractfield(RLinks,'KOMME');
Henrik Grythe's avatar
Henrik Grythe committed
48
uKomm     = unique(KommunerS);
49
50
51
52

% Congestion weights.  Traffic part delayed (TPD)
% Combine the different congestion levels by volume to determine the
% actual EF:
53
54
55
56
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];
57

58
59
60
61
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];
62

63
64
65
66
67
68
69
70
71
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];
72

Henrik Grythe's avatar
Henrik Grythe committed
73
% EXTRACT 
74
75
76
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)
Henrik Grythe's avatar
Henrik Grythe committed
77
78
LEN     = extractfield(RLinks,'DISTANCE');                  % Length of Road      (in kilometres)

79
80
81
82
83
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)
RU      = extractfield(RLinks,'RUSH_DELAY');
84

85
86
87
88
89
90
% 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;
91

92
93
IDO     = extractfield(RLinks,'IDO');
KOMS    = extractfield(RLinks,'KOMMS');
94

95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
% 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)));
114
        
115
116
117
118
119
120
121
122
123
        % 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));
        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));
124
        
125
126
127
128
129
        % Calculate the annual driving distance on the roadlink in
        % Kilometers:
        DL = L(r) *dayInYear* LEN(r) *IDO(r);
        HL = H(r) *dayInYear* LEN(r) *IDO(r);
        BL = B(r) *dayInYear* LEN(r) *IDO(r);
130
        
131
132
133
134
        % Divide into given congestion levels
        Lmy(idx) =   DL*TPD_L(ru(r),:)';
        Hmy(idx) =   HL*TPD_H(ru(r),:)';
        Bmy(idx) =   BL*TPD_B(ru(r),:)';
135
        
136
137
138
139
        % Calculate the emissions:
        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))');
140
141
        
        
Henrik Grythe's avatar
Henrik Grythe committed
142
        % Calculate roadLink Emission Factor: (Not used as of now)
143
144
145
146
147
        Link_Light_emission_factor(r) = EMISS_L(r)/sum(Lmy(idx));
        Link_Heavy_emission_factor(r) = EMISS_H(r)/sum(Hmy(idx));
        Link_Buses_emission_factor(r) = EMISS_B(r)/sum(Bmy(idx));
        Link_emission_factor(r)       = (EMISS_L(r)+EMISS_H(r)+EMISS_B(r))/(sum(Lmy(idx))+sum(Hmy(idx))+sum(Bmy(idx)));
        if rem(r,10000)==0; fprintf('Roads Calculated %i/%i\n  L:%f\n  H:%f\n  B:%f\nEF: %f\n',r,length(RLinks),sum(EMISS_L)*1e-9,sum(EMISS_H)*1e-9,sum(EMISS_B)*1e-9,Link_emission_factor(r)); end
148
    end
149
150
151
152
153
154
155
156
157
158
    TLinks = struct2table(RLinks);
    TLinks.LetEM =  EMISS_L;
    TLinks.HeaEM =  EMISS_H;
    TLinks.BusEM =  EMISS_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);
    
Henrik Grythe's avatar
Henrik Grythe committed
159
    % Print some statistical output
160
    fprintf('\n---- NORGE --- \n')
Henrik Grythe's avatar
Henrik Grythe committed
161
162
163
    fprintf('---- Lette   %11.1f   (1000)Ton %s (%3.0f%%)\n'  ,1e-9*sum(EMISS_L),char(comps(com)),100*nansum(EMISS_L)/nansum(EMISS_L + EMISS_H + EMISS_B))
    fprintf('---- Tunge   %11.1f   (1000)Ton %s (%3.0f%%)\n'  ,1e-9*nansum(EMISS_H),char(comps(com)),100*nansum(EMISS_H)/nansum(EMISS_L + EMISS_H + EMISS_B))
    fprintf('---- Busser  %11.1f   (1000)Ton %s (%3.0f%%)\n'  ,1e-9*nansum(EMISS_B),char(comps(com)),100*nansum(EMISS_B)/nansum(EMISS_L + EMISS_H + EMISS_B))
164
    fprintf('---- Totalt  %11.1f   (1000)Ton %s \n\n',1e-9*nansum(EMISS_B+EMISS_H+EMISS_L),char(comps(com)))
165
    
166
167
168
169
170
171
172
173
174
175
    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)
176
177
178
    fprintf('     Light Traffic L=%7.1f g/Km\n',1e-6*sum(EMISS_L)/LW)
    fprintf('     Heavy Traffic H=%7.1f g/Km\n',1e-6*sum(EMISS_H)/HW)
    fprintf('     Buses Traffic B=%7.1f g/Km\n',1e-6*sum(EMISS_B)/BW)
Henrik Grythe's avatar
Henrik Grythe committed
179
    fprintf('--- %s --->\n',char(comps(com)))
180
end
Henrik Grythe's avatar
Henrik Grythe committed
181
182

save(ofiles.MatlabOutput,'RLinks','-append')
183
fprintf('Added RLinks to output\n')
184
Sn     = RLinks;
185
end