NERVE_Emissions_Statistical_Oputput.m 8.29 KB
Newer Older
1
2
3
4
5
function NERVE_Emissions_Statistical_Oputput()
% Data needed by the MijlDirektoratet sheet
comps  = [{'N2O'}];
comps  = [{'CO2'},{'FC'},{'CH4'},{'N2O'}];
comps  = [{'CO2'},{'FC'},{'CH4'}];
6
7
8
%comps  = [{'CO2'},{'FC'}];
%comps  = [{'CH4'}];

9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
%x NV(kom,veh) %#
%x TD(kom,veh) %km
%x ORdComp(kom,veh) % frac
%x FRdComp(kom,veh) % frac
%x L_IN(kom)        % km/yr
%x H_IN(kom)        % km/yr
%x B_IN(kom)        % km/yr

%Check the below values against Roads: 
%s L_FROM(kom)       % km/yr
%s H_FROM(kom)       % km/yr
%s B_FROM(kom)       % km/yr
%s EF_IN(kom,veh)    % g/km
%s EM_IN(kom,veh)    % g/yr
%s EF_FROM(kom,veh)  % g/km
%s EM_FROM(kom,veh)  % g/yr

% Files that must be read once only
% 'Model_Vehicles_Merge_SSB_and_HBEFA_Vehicles.xlsx'
TM = readtable('Input/Model_Vehicles_Merge_SSB_and_HBEFA_Vehicles.xlsx','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;

for com = 1:length(comps)
    % files that must be read per species
35

36
37
    fprintf('%s\n',char(comps(com)))

38
    for Tyear = 2009:2019
39
40
41
42
43
44
45
46
47
48
49
50
51
        % files that must be read for each year and component!
        EF_File = sprintf('Temp/EF_On_AllRoadCond_Municipality_%i_%s.mat',Tyear,char(comps(com)));
        load(EF_File)

        try
            R_EF_File = sprintf('Temp/EFA_Table_MODEL_%s_Bio%i.mat',char(comps(com)),Tyear);
            load(R_EF_File,'TFout'); % TFout
            fprintf('Found Bio File!   ')
        catch
            R_EF_File = sprintf('Temp/EFA_Table_MODEL_%s.mat',char(comps(com)));
            load(R_EF_File,'TFout'); % TFout
        end
        fprintf('%i',Tyear)
52
53

        % files that must be read per year
54
        munFile    = sprintf('Temp/Municipal_Traffic_Exchange_%i.mat',Tyear);
55
        RdDistFile = sprintf('Output/RoadTypeDistanceMunicipal%i.mat',Tyear);
56
57
        SSBCPfile  = sprintf('Temp/SSB_CarPark_%i.mat',Tyear);
        load(munFile)    % TrafficIN TrafficFROM kmne
58
        load(RdDistFile) % KDD
59
        load(SSBCPfile)  % Vehicle_dist
60
61
62
        
                
        RLinks = shaperead(sprintf('Output/Traffic_Emissions_%i',Tyear));
63
64

        DaysInYear = datenum([Tyear+1 1 1 0 0 0])-datenum([Tyear 1 1 0 0 0]);
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
        
        KOMM       = extractfield(RLinks,'KOMMS');
        KOMMe      = extractfield(RLinks,'KOMME');
        IDO        = extractfield(RLinks,'IDO');
        L_ADT      = extractfield(RLinks,sprintf('L_ADT%i',Tyear));
        H_ADT      = extractfield(RLinks,sprintf('H_ADT%i',Tyear));
        B_ADT      = extractfield(RLinks,sprintf('B_ADT%i',Tyear));
        DISTANCE   = extractfield(RLinks,'DISTANCE');
        LW         = DaysInYear*L_ADT.*DISTANCE;
        HW         = DaysInYear*H_ADT.*DISTANCE;
        BW         = DaysInYear*B_ADT.*DISTANCE;
        ND_NERVE_L = sum(LW);
        ND_NERVE_H = sum(HW);
        ND_NERVE_B = sum(BW);
        ND_NERVE   = ND_NERVE_L + ND_NERVE_H + ND_NERVE_B;

        % Ø----------------------------------------------------------------
        % find *L_IN(kom)*
        ukomm = unique(KOMM);
        for k=1:length(ukomm)
            f  = find(KOMM  == ukomm(k));
            f2 = find((KOMMe == ukomm(k))&(KOMM ~= ukomm(k)) );
            L_IN(k) = sum(LW(f).*IDO(f))+sum(LW(f2).*(1-IDO(f2)));
            H_IN(k) = sum(HW(f).*IDO(f))+sum(HW(f2).*(1-IDO(f2)));
            B_IN(k) = sum(BW(f).*IDO(f))+sum(BW(f2).*(1-IDO(f2)));
        end
        
        L_EM      = extractfield(RLinks,sprintf('EM_L_%s',char(comps(com))));
        H_EM      = extractfield(RLinks,sprintf('EM_H_%s',char(comps(com))));
        B_EM      = extractfield(RLinks,sprintf('EM_B_%s',char(comps(com))));
        
        EM = L_EM+H_EM+B_EM; 
        for k=1:length(ukomm)
            f  = find(KOMM  == ukomm(k));
            f2 = find((KOMMe == ukomm(k))&(KOMM ~= ukomm(k)) );
            EM_INr(k) = sum(EM(f).*IDO(f))+sum(EM(f2).*(1-IDO(f2)));
101
102
%             EM_INr(k) = sum(EM(f).*IDO(f))+sum(EM(f2).*(1-IDO(f2)));
%             EM_INr(k) = sum(EM(f).*IDO(f))+sum(EM(f2).*(1-IDO(f2)));
103
        end
104
105
106
        pdata =1e-9*[sum(L_EM),sum(L_EM),sum(L_EM),sum(L_EM+H_EM+B_EM)];
        fprintf('\tLIGHT:%6.2f,  HEAV:Y%6.2f,  BUS:%6.2f,  TOT:%6.2f\n',pdata)

107
108
109
110
111
112
        % Vehicle_dist
        %------------------------------------------------------------------
        % find *NV(kom,veh)* from:: SSB data
        % find *TD(kom,veh)* from:: SSB data
        % find *OrdComp(kom,veh)* from:: SSB data
        % find *FrdComp(kom,veh)* from:: SSB data
113
114
        NV      = Vehicle_dist.modelNV;
        TD      = Vehicle_dist.modelTD;
115
116
        ORdComp = Vehicle_dist.Vdist;
        FRdComp = Vehicle_dist.VdistFROM;        
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
        
        %------------------------------------------------------------------
        % Estimate *EF_IN(kom,Veh)  * based on:::: EF_IN and exchange
        
        KDD_L = KDD.TraffDataL;
        KDD_H = KDD.TraffDataH;
        KDD_B = KDD.TraffDataB;
        
        % DOES NOT INCLUDE BIO!
        RdEFs = table2array(TFout(:,8:end));       
        vehicles = TFout.Properties.VariableNames(8:end);
        for k=1:length(ukomm)
            LtrafficSitW_IN = KDD_L(:,k)/sum(KDD_L(:,k));
            HtrafficSitW_IN = KDD_H(:,k)/sum(KDD_H(:,k));
            BtrafficSitW_IN = KDD_B(:,k)/sum(KDD_B(:,k));
            for veh = 1:length(vehicles)
                if (TM.Model_Class(veh) ==1 || TM.Model_Class(veh) ==2)
                    EF_IN(k,veh) = sum(RdEFs(:,veh).*LtrafficSitW_IN);
                elseif(TM.Model_Class(veh) ==5 ||TM.Model_Class(veh) ==6 ||TM.Model_Class(veh) ==7)
                    EF_IN(k,veh) = sum(RdEFs(:,veh).*HtrafficSitW_IN);
                elseif(TM.Model_Class(veh) ==3 ||TM.Model_Class(veh) ==4)
                    EF_IN(k,veh) = sum(RdEFs(:,veh).*BtrafficSitW_IN);
                else
                    fprintf('fail\n')
                end
            end
        end
        
        %------------------------------------------------------------------
        % Estimate *EM_IN(kom,Veh)  * based on:::: EF_IN and exchange
        
        for k=1:length(ukomm)
            for veh = 1:length(vehicles)
                type = TM.Model_Class(veh);
                if type <= 2
152
                    EM_IN(k,veh) = L_IN(k)*ORdComp(k,veh)*EF_IN(k,veh);
153
                elseif (type == 5 || type == 6 || type == 7)
154
                    EM_IN(k,veh) = H_IN(k)*ORdComp(k,veh)*EF_IN(k,veh);
155
                elseif (type == 3 || type == 4)
156
                    EM_IN(k,veh) = B_IN(k)*ORdComp(k,veh)*EF_IN(k,veh);
157
158
159
160
                end
            end
        end

161
162
163
        pdata =[1e-9*nansum(nansum(EM_IN)),100*nansum(nansum(EM_IN))/sum(L_EM+H_EM+B_EM)];
        fprintf('\t\t\t\t\tTOT_IN %8.1f diff %4.2f%%\n',pdata)
 
164
165
166
167
168
169
170
171
172
173
174
175
% END IN
%--------------------------------------------------------------------------
% This method assigns the total Emissions based on the exchange.

        for komm = 1:size(TrafficFROM,2)
            I            = find(TrafficFROM(komm,:)>0);
            Tshare       = (TrafficFROM(komm,I)/100);
            L_FROM(komm) = nansum(L_IN(I).*Tshare);
            H_FROM(komm) = nansum(H_IN(I).*Tshare);
            B_FROM(komm) = nansum(B_IN(I).*Tshare);
        end
        
176
        EF_FROM = zeros(size(EF_IN));
177
178
179
180
181
182
183
184
185
186
187
188
        for komm=1:size(EF_IN,1)
            for Veh=1:size(EF_IN,2)
                 EF_FROM(komm,Veh) = EF_FROM(komm,Veh)+nansum(EF_IN(:,Veh).*(TrafficFROM(:,komm)/100));
             end
        end
 
        %--------------------------------------------------------------------------
        % create *L_FROM(kom)* based on:::: L_IN and exchange
        for k=1:length(ukomm)
            for veh = 1:length(vehicles)
                type = TM.Model_Class(veh);
                if type <= 2
189
                    EM_FROM(k,veh) = L_FROM(k)*FRdComp(k,veh)*EF_FROM(k,veh);
190
                elseif (type == 5 || type == 6 || type == 7)
191
                    EM_FROM(k,veh) = H_FROM(k)*FRdComp(k,veh)*EF_FROM(k,veh);
192
                elseif (type == 3 || type == 4)
193
                    EM_FROM(k,veh) = B_FROM(k)*FRdComp(k,veh)*EF_FROM(k,veh);
194
195
196
197
198
                end
            end
        end

        fileout = sprintf('Output/NERVE_output_%s_%04i.mat',char(comps(com)),Tyear);
199
        fprintf('Processed Emissions for %s year %i\n',char(comps(com)),Tyear)
200
201
202
203
204
205
        save(fileout,'NV','TD','L_IN','H_IN','B_IN','L_FROM','H_FROM','B_FROM','ORdComp','FRdComp','EF_IN','EF_FROM','EM_IN','EM_FROM')
    end
end


end