NERVE_Emissions_Statistical_Oputput.m 8.56 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
%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

26
27
28
29
30
31
32
33
if ispc
    bp = 'N:\Inby';
else
   bp = '/storage/nilu/Inby';
end
pathn = sprintf('%s/Emission_Group/Emission_Models/HEDGE/',bp);


34
35
% Files that must be read once only
% 'Model_Vehicles_Merge_SSB_and_HBEFA_Vehicles.xlsx'
36
TM = readtable(sprintf('%sInput/Model_Vehicles_Merge_SSB_and_HBEFA_Vehicles.xlsx',pathn),'Sheet','MODEL');
37
38
39
40
41
42
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
43

44
45
    fprintf('%s\n',char(comps(com)))

46
    for Tyear = 2009:2019
47
        % files that must be read for each year and component!
48
        EF_File = sprintf('%sTemp/EF_On_AllRoadCond_Municipality_%i_%s.mat',pathn,Tyear,char(comps(com)));
49
50
51
        load(EF_File)

        try
52
            R_EF_File = sprintf('%sTemp/EFA_Table_MODEL_%s_Bio%i.mat',pathn,char(comps(com)),Tyear);
53
54
55
            load(R_EF_File,'TFout'); % TFout
            fprintf('Found Bio File!   ')
        catch
56
            R_EF_File = sprintf('%sTemp/EFA_Table_MODEL_%s.mat',pathn,char(comps(com)));
57
58
            load(R_EF_File,'TFout'); % TFout
        end
59
        fprintf('%i\n',Tyear)
60
61

        % files that must be read per year
62
63
64
        munFile    = sprintf('%sTemp/Municipal_Traffic_Exchange_%i.mat',pathn,Tyear);
        RdDistFile = sprintf('%sOutput/RoadTypeDistanceMunicipal%i.mat',pathn,Tyear);
        SSBCPfile  = sprintf('%sTemp/SSB_CarPark_%i.mat',pathn,Tyear);
65
        load(munFile)    % TrafficIN TrafficFROM kmne
66
        load(RdDistFile) % KDD
67
        load(SSBCPfile)  % Vehicle_dist
68
        
69
70
71
        fprintf('Reading large RoadLink file...')
        RLinks = shaperead(sprintf('%sOutput/Traffic_Emissions_%i',pathn,Tyear));
        fprintf('Done\n')
72
73

        DaysInYear = datenum([Tyear+1 1 1 0 0 0])-datenum([Tyear 1 1 0 0 0]);
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
101
102
103
104
105
106
107
108
109
        
        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)));
110
111
%             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)));
112
        end
113
        pdata =1e-9*[sum(L_EM),sum(H_EM),sum(B_EM),sum(L_EM+H_EM+B_EM)];
114
115
        fprintf('\tLIGHT:%6.2f,  HEAV:Y%6.2f,  BUS:%6.2f,  TOT:%6.2f\n',pdata)

116
117
118
119
120
121
        % 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
122
123
        NV      = Vehicle_dist.modelNV;
        TD      = Vehicle_dist.modelTD;
124
125
        ORdComp = Vehicle_dist.Vdist;
        FRdComp = Vehicle_dist.VdistFROM;        
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
152
153
154
155
156
157
158
159
160
        
        %------------------------------------------------------------------
        % 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
161
                    EM_IN(k,veh) = L_IN(k)*ORdComp(k,veh)*EF_IN(k,veh);
162
                elseif (type == 5 || type == 6 || type == 7)
163
                    EM_IN(k,veh) = H_IN(k)*ORdComp(k,veh)*EF_IN(k,veh);
164
                elseif (type == 3 || type == 4)
165
                    EM_IN(k,veh) = B_IN(k)*ORdComp(k,veh)*EF_IN(k,veh);
166
167
168
169
                end
            end
        end

170
171
172
        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)
 
173
174
175
176
177
178
179
180
181
182
183
184
% 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
        
185
        EF_FROM = zeros(size(EF_IN));
186
187
188
189
190
191
192
193
194
195
196
197
        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
198
                    EM_FROM(k,veh) = L_FROM(k)*FRdComp(k,veh)*EF_FROM(k,veh);
199
                elseif (type == 5 || type == 6 || type == 7)
200
                    EM_FROM(k,veh) = H_FROM(k)*FRdComp(k,veh)*EF_FROM(k,veh);
201
                elseif (type == 3 || type == 4)
202
                    EM_FROM(k,veh) = B_FROM(k)*FRdComp(k,veh)*EF_FROM(k,veh);
203
204
205
206
                end
            end
        end

207
        fileout = sprintf('%sOutput/NERVE_output_%s_%04i.mat',pathn,char(comps(com)),Tyear);
208
        fprintf('Processed Emissions for %s year %i\n',char(comps(com)),Tyear)
209
210
211
212
213
214
        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