Emissions_Calculations_SSB.m 8.43 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 input
Henrik Grythe's avatar
Henrik Grythe committed
24
25
26
fprintf('---------------------------------------------------------------\n')
fprintf('in Emissions_Calculations_SSB    *\n')
fprintf('---------------------------------------------------------------\n')
27
28

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

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

40
TM = readtable(input.files.SSB_Vehicle_dist,'Sheet','MODEL');
41
42
43
44
45
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)
46
47


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

% Congestion weights.  Traffic part delayed (TPD)
% Combine the different congestion levels by volume to determine the
% actual EF:
56
%-----------------------------------
Henrik Grythe's avatar
Henrik Grythe committed
57
TPD  =[1.0, 0.0, 0.0, 0.0, 0.0;...
Henrik Grythe's avatar
Henrik Grythe committed
58
59
60
       0.6, 0.3, 0.1, 0.0, 0.0;...
       0.6, 0.2, 0.1, 0.1, 0.0;...
       0.2, 0.2, 0.2, 0.2, 0.2];
Henrik Grythe's avatar
Henrik Grythe committed
61
62
63
TPD_L = TPD;
TPD_H = TPD;
TPD_B = TPD;
64

Henrik Grythe's avatar
Henrik Grythe committed
65
% EXTRACT 
66
67
68
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
69
70
LEN     = extractfield(RLinks,'DISTANCE');                  % Length of Road      (in kilometres)

71
72
73
74
75
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');
76

Henrik Grythe's avatar
Henrik Grythe committed
77
78
79
80
% Calculated Properties
LW  = 1e-6*sum(L.*LEN)*dayInYear;
HW  = 1e-6*sum(H.*LEN)*dayInYear;
BW  = 1e-6*sum(B.*LEN)*dayInYear;
Henrik Grythe's avatar
Henrik Grythe committed
81
82
83
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)));
Henrik Grythe's avatar
Henrik Grythe committed
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
for com = 1:length(comps)
Henrik Grythe's avatar
Henrik Grythe committed
96
97
98
99
100
    % 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));
101
    fprintf('<--- %s ---\n',char(comps(com)))
102
    
103
    fprintf('Loading large file\n...')
Henrik Grythe's avatar
Henrik Grythe committed
104
105
    tmpfile = sprintf('%sEF_On_AllRoadCond_Municipality_%i_%s.mat',tfold,Tyear,char(comps(com)));
    load(tmpfile)
106
107
    TEF = Trout;

108
109
110
111
112
113
114
115
116
    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)));
117
        
118
119
120
121
122
123
124
125
126
        % 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));
127
        
128
129
130
131
132
        % 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);
133
        
134
135
136
137
        % 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),:)';
138
        
139
140
141
142
        % 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))');
143
        
Henrik Grythe's avatar
Henrik Grythe committed
144
        % Calculate roadLink Emission Factor: (Not used as of now)
145
146
147
148
149
        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
150
    end
151
152
153
154
155
    TLinks = struct2table(RLinks);
    TLinks.LetEM =  EMISS_L;
    TLinks.HeaEM =  EMISS_H;
    TLinks.BusEM =  EMISS_B;
    
Henrik Grythe's avatar
Henrik Grythe committed
156
157
158
    TLinks.Properties.VariableNames(find(ismember(TLinks.Properties.VariableNames,'LetEM'))) = {sprintf('EM_L_%s',char(comps(com)))};
    TLinks.Properties.VariableNames(find(ismember(TLinks.Properties.VariableNames,'HeaEM'))) = {sprintf('EM_H_%s',char(comps(com)))};
    TLinks.Properties.VariableNames(find(ismember(TLinks.Properties.VariableNames,'BusEM'))) = {sprintf('EM_B_%s',char(comps(com)))};
159
160
    RLinks = table2struct(TLinks);
    
Henrik Grythe's avatar
Henrik Grythe committed
161
    % Print some statistical output
162
    fprintf('\n---- NORGE --- %i %s \n',Tyear,char(comps(com)))
Henrik Grythe's avatar
Henrik Grythe committed
163
164
165
    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))
166
    fprintf('---- Totalt  %11.1f   (1000)Ton %s \n\n',1e-9*nansum(EMISS_B+EMISS_H+EMISS_L),char(comps(com)))
167
    
168
    
169
170
171
    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
172
173
174
175
    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)

Henrik Grythe's avatar
Henrik Grythe committed
176
    fprintf('--- %s --->\n',char(comps(com)))
177
end
Henrik Grythe's avatar
Henrik Grythe committed
178
179

save(ofiles.MatlabOutput,'RLinks','-append')
180
fprintf('Added RLinks to output\n')
181
Sn     = RLinks;
182
end