Emissions_Calculations_SSB.m 12.8 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
23
global tfold Tyear SSB_Vehicle_dist comps RLinks Vehicle_dist Vehicle_weight
global debug_mode
24

25
26
27
%--------------------------------------------------------------
% MUNICIPALITY EMISSION FACTORS FOR ALL ROADS
%--------------------------------------------------------------
28
29
30
31



%--------------------------------------------------------------
32
% ROAD LINK
33
34
35
36
%--------------------------------------------------------------
% extract fields needed for calculations
fprintf('\nExtracts necessary fields   \n')
fprintf('Traffic and Vehicle Year:  %i \n',Tyear)
37
fprintf('RoadLinks               :  %i \n',size(RLinks,1))
38
39
40
41
42

file = sprintf('%s%s',tfold,'roads');
fprintf('Warning,using roads file not produced by NERVE model\n%s\n',file)
load(file)

43
44
45
46
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;
47
48
49
50




51
52
53
54
55
56
% 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');
57

58
59
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)
60
61

% Make a cleaned version for writing to file
62
63
Tout = struct2table(RLinks);

64

65
66
67
68
% roads
KommunerS = extractfield(RLinks,'KOMMS');
KommunerE = extractfield(RLinks,'KOMME');
uKomm    = unique(KommunerS);
69
70
71
72

% Congestion weights.  Traffic part delayed (TPD)
% Combine the different congestion levels by volume to determine the
% actual EF:
73
74
75
76
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];
77

78
79
80
81
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];
82

83
84
85
86
87
88
89
90
91
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];
92
93
94
95



for com = 1:length(comps)
96
97
98
99
    fprintf('<--- %s ---\n',char(comps(com)))

    TEF = readtable('OnRoadEF_RoadClasses.xlsx','Sheet',sprintf('%s_%i',char(comps(com)),Tyear),'PreserveVariableNames',1);    
    % load(oEFfile,'TFout','roads');
100
    
101
    ttef = TEF(:,1:7);
102
    
103
104
105
106
107
108
109
110
111
112
    % All road links emissions Light / Heavy / Bus
    nEMISS_L = zeros(size(RLinks));
    nEMISS_H = zeros(size(RLinks));
    nEMISS_B = zeros(size(RLinks));

    
    % 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));    
113
    
114
115
116
117
118
    for komm = 1:length(uKomm)
        % extract all necessary roadLinks:
        f1 = find(KommunerS == uKomm(komm));
        f2 = find(KommunerE == uKomm(komm));
        f  = unique([f1,f2]);
119
        
120
121
122
        MRLinks = RLinks(f);        
        fprintf('\n-%3i Kommune_%04i --- \n',komm,Vehicle_dist.D1_KommNr(komm))
        fprintf('     Has %i Roads     \n',length(f))
123
        
124
125
126
        % extract all roadEFs:
        ref = find(contains(TEF.Properties.VariableNames,{sprintf('%i',uKomm(komm))}));
        Tef = [ttef,TEF(:,ref)];
127
        
128
129
130
131
132
133
134
135
136
137
138
139
140
        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))));

        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');
141
        
142
143
144
145
146
147
        % 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;
148
        
149
150
151
152
153
154
155
        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)))
156
157
        end
        
158
159
160
161
162
163
164
165
166
167
        % *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);
        
        % 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));
168
        
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
        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)+1,:)';
                Hmy(idx) = Hmy(idx) + dayInYear *IDO(r)* H(r) *LEN(r)*TPD_H(RU(r)+1,:)';
                Bmy(idx) = Bmy(idx) + dayInYear *IDO(r)* B(r) *LEN(r)*TPD_B(RU(r)+1,:)';
            else
                Lmy(idx) = Lmy(idx) + dayInYear *IDO(r)* L(r) *LEN(r)*TPD_L(RU(r)+1,:)';
                Hmy(idx) = Hmy(idx) + dayInYear *IDO(r)* H(r) *LEN(r)*TPD_H(RU(r)+1,:)';
                Bmy(idx) = Bmy(idx) + dayInYear *IDO(r)* B(r) *LEN(r)*TPD_B(RU(r)+1,:)';
187
            end
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
            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));
        
        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))
213
        end
214
215
216
217
218
219
220
221
222
223
224
225
226
227
        
        Tef.Light_DD = Lmy;
        Tef.Heavy_DD = Hmy;
        Tef.Buses_DD = Bmy;
        
        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)))
228
    
229
230
231
232
    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');
233
    
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
    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)

    TLinks = struct2table(RLinks);
    TLinks.LetEM =  nEMISS_L;
    TLinks.HeaEM =  nEMISS_H;
    TLinks.BusEM =  nEMISS_B;
249
    
250
251
252
253
    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);
254
end
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
Sn = RLinks;
end

% Adj      = SPD/70;
% Km       = LEN*1e-3;
% wearL    = Adj.*MP_wear_Light.*(365*L).*Km;
% wearH    = Adj.*MP_wear_Heavy.*(365*(H+B)).*Km;
% Tout.EM_TW10  = MP_airbornefraction*(wearL + wearH)';
%
% % fileout = sprintf('%sOutput_%i.mat',tfold,Tyear);
% % %save(fileout,'Class_Weight','TPD','aEM_L','aEM_H','aEM_B','aEM_2W','aEM_Tot','nEM_L','nEM_H','nEM_B','nEM_2W','nEM_Tot')
% % save(fileout,'Class_Weight','TPD','aEM_L','aEM_H','aEM_B','aEM_2W','aEM_Tot')
%
%
% % Finally write the output to a shape file.
% fprintf('Re-Structuring Roads...\n')
% Sn = table2struct(Tout);
%end