Commit baf26d18 authored by Henrik Grythe's avatar Henrik Grythe
Browse files

Fixed (1?) bug inducing NaN emissions on certain roads

parent 8bd6f12e
......@@ -4,7 +4,8 @@ function Emission_Factor_Model_group_HBEFA()
% 22.10.2020 -Henrik Grythe
% Kjeller NILU
%--------------------------------------------------------------------------
global tfold comps SSB_Vehicle_dist Vehicle_source
global tfold comps SSB_Vehicle_dist Vehicle_source debug_mode
global EFA EF_AVG roads
% global rfid
......@@ -19,7 +20,7 @@ Vehicle_weight = 'Uniform';
fprintf('\tin Emission_Factor_Model_group_HBEFA\n')
fprintf('--- Grouping Vehicles according to \n--- Model SSB HBEFA merger specifications \n\n')
fprintf('--- Grouping Vehicles according to ---\n--- Model SSB HBEFA merger specifications ---\n%s \n',SSB_Vehicle_dist)
T = readtable(SSB_Vehicle_dist,'Sheet','HBEFA778toMODEL');
......@@ -27,8 +28,10 @@ T = readtable(SSB_Vehicle_dist,'Sheet','HBEFA778toMODEL');
modelN = unique(T.ModelNumber);
fprintf('Found %i Model Vehicles \n',length(modelN))
ConVfile = sprintf('%s/EFA_conversion.xlsx',tfold);
for com =1:length(comps)
fprintf('Spec| %s\n',char(comps(com)))
fprintf('<--- \nSpec| %s\n',char(comps(com)))
iEFfile = sprintf('%s/EFA_matrix41_RAW_%s',tfold,char(comps(com)));
oEFfile = sprintf('%s/EFA_matrix41_MODEL_%s',tfold,char(comps(com)));
load(iEFfile);
......@@ -40,7 +43,7 @@ for com =1:length(comps)
% Defines the EFA
EFA = zeros(length(modelN),size(EF_AVG,2),size(EF_AVG,3),size(EF_AVG,4),size(EF_AVG,5),size(EF_AVG,6));
for mod=1:length(modelN)
for mod = 1:length(modelN)
Tsub = T(T.ModelNumber==modelN(mod),:);
%##################
% errrrrorrrr Need to add column "HBEFA_Num" to the sheet!!!!!
......@@ -84,7 +87,10 @@ for com =1:length(comps)
end
end % switch
end % for mod
Emission_Factor_Test_HBEFA_to_MODEL_conversion(com,ConVfile)
save(oEFfile,'EFA','T');
fprintf('Saved a temp-file for Emission Factors Model:\n%s\n',oEFfile)
......
function Emission_Factor_Process_HBEFA_Matrix_Raw()
%--------------------------------------------------------------------------
% 22.10.2020 -Henrik Grythe
% Kjeller NILU
%--------------------------------------------------------------------------
% Script for reading and processing HBEFA .csv file of subsegments.
%
%
% Basically taken from NERVE, but includes some updates for the format of
% HBEFA4.1.
%
% Makes a 6-D matix of the emissions factor, convenient for speedy
% HBEFA4.1. Makes a 6-D matix of the emissions factor, convenient for speedy
% processing.
% 22.10.2020 -Henrik Grythe
% Kjeller NILU
%--------------------------------------------------------------------------
%
% % Path to emissions
% HBEFA_path = '/storage/nilu/Inby/Emission_Group/Emission_Factors/Traffic/HBEFA_Raw_data/HOT_HBEFA_41/';
%
% % List the avalable species SLA have downloaded Emission factors for
% comps = [{'NOx'},{'CO2'},{'BC'},{'Be'},{'CO'},{'HC'},{'NMHC'},{'NO2'},{'PM'},{'PN'},{'NH3'},{'N2O'},{'CH4'},{'FC'},{'FC_MJ'}];
global do_preProcessing_HBEFA debug_mode
global HBEFA_path tfold comps
fprintf('* call PreProcess_HBEFA_Matrix_Raw *\n')
fprintf('\tEmission_Factor_Process_HBEFA_Matrix_Raw *\n')
if ~do_preProcessing_HBEFA
fprintf('We Already Have Necessary HBEFA Emission Factors \n')
fprintf('... Continuing Without new HBEFA calculations \n ... \n')
......@@ -44,9 +35,9 @@ NfldList = [{'Year'},{'TrafficScenario'},{'RoadCat'},{'IDSubsegment'},{'KM'},{'x
{'EFA_weighted_100_'},{'EFA_WTT'},{'EFA_WTT_0_'},{'EFA_WTT_100_'},{'EFA_WTW'},{'EFA_WTW_0_'},{'EFA_WTW_100_'},{'AmbientCondPattern'}];
% Loop over all compounds in comps list.
redo = 1;
for com = 1:length(comps)
for com = 1:length(comps)
% assumed a naming convention
ifile1 = sprintf('%sEFA_HOT_Subsegm_%s.csv',HBEFA_path,char(comps(com)));
ifile2 = sprintf('%sHGV/EFA_HOT_Subsegm_%s_hgv.csv',HBEFA_path,char(comps(com)));
......@@ -73,14 +64,13 @@ for com = 1:length(comps)
end
ShoulBeNumeric = logical(ShoulBeNumeric);
% Thers a problem in some of the files where columns are read as
% different types of variablse. Therefore access each table and check
% if the right columns are numeric.
numericVars1 = varfun(@isnumeric,Tn1,'output','uniform');
numericVars2 = varfun(@isnumeric,Tn2,'output','uniform');
% table 1 first
% table 1 first
NT1 =table;
for i = 1:width(Tn1)
fprintf('%02i :: %s\n',i,char(Tn1.Properties.VariableNames(i)))
......@@ -112,11 +102,14 @@ for com = 1:length(comps)
end
NT2.Properties.VariableNames = Tn2.Properties.VariableNames;
% Now that both fields have the same numeric columns,
% combine the two tables
Tn = [NT1;NT2];
% Test that
Old = [Tn1(1,:);Tn2(1,:);Tn1(end,:);Tn2(end,:)]
New = [NT1(1,:);NT2(1,:);NT1(end,:);NT2(end,:)]
if debug_mode
% Test that
Old = [Tn1(1,:);Tn2(1,:);Tn1(end,:);Tn2(end,:)]
New = [NT1(1,:);NT2(1,:);NT1(end,:);NT2(end,:)]
end
writetable(Tn,ifile3)
clear Tn1 Tn2 NT1 NT2
......@@ -168,6 +161,7 @@ for com = 1:length(comps)
T.ENGINE = Tn.SizeClasse;
T.EUROCAT = Tn.EmConcept;
traffCon = unique(T.ROAD_COND);
% Numerical columns: with silightly simplified names.
T.GRADIENT = Gradient';
......@@ -382,7 +376,7 @@ for com = 1:length(comps)
% save(sprintf('EFA_matrix_RAW_%s',char(comps(com))),'roads','EF_AVG','EF_000','EF_100')
oEFfile = sprintf('%s/EFA_matrix41_RAW_%s',tfold,char(comps(com)));
save(oEFfile,'roads','EF_AVG','EF_000','EF_100');
save(oEFfile,'roads','EF_AVG','EF_000','EF_100','traffCon');
fprintf('Saved temporary file\n%s\n',oEFfile)
clear Gradient Tn T RTID roads EF_*
......
function Emission_Factor_Test_HBEFA_to_MODEL_conversion(com,ConVfile)
global EFA EF_AVG comps roads
fprintf('\tin Emission_Factor_Test_HBEFA_to_MODEL_conversion\n')
Traff = table;
k = 1;
for roa = 1:size(EF_AVG,2)
for spd = 1:size(EF_AVG,3)
for dec = 1:size(EF_AVG,4)
for cog = 1:size(EF_AVG,5)
for urb = 1:size(EF_AVG,6)
Tb = table;
Tb.name = {sprintf('%s/%s-%i/%i%%/%s',char(roads.RoadEnv(urb)),char(roads.RoadType(roa)),roads.RoadSpeeds(spd),roads.RoadGradient(dec),char(roads.Congestion(cog)))};
Tb.roadTypeID = roa;
Tb.roadType = roads.RoadType(roa);
Tb.roadSpeedID = spd;
Tb.roadSpeed = roads.RoadSpeeds(spd);
Tb.roadGradID = dec;
Tb.roadGrad = roads.RoadGradient(dec);
Tb.roadCongID = cog;
Tb.roadCong = roads.Congestion(cog);
Tb.roadCongID = urb;
Tb.roadCong = roads.RoadEnv(urb);
VEH = zeros(1,size(EF_AVG,1));
VEHp = zeros(1,size(EF_AVG,1));
for veh = 1:size(EF_AVG,1)
if ~isnan(EF_AVG(veh,roa,spd,dec,cog,urb))
VEH(veh) = VEH(veh)+1;
end
if EF_AVG(veh,roa,spd,dec,cog,urb)>0
VEHp(veh) = VEHp(veh)+1;
end
end
nums1(k) = sum(VEH);
nums2(k) = sum(VEHp);
k = k+1;
Traff = [Traff;Tb];
end
end
end
end
end
Traff.NumVehEF_HBEFA = nums1';
Traff.NumVehEFz_HBEFA = nums2';
Traff2 = table;
k = 1;
for roa = 1:size(EFA,2)
for spd = 1:size(EFA,3)
for dec = 1:size(EFA,4)
for cog = 1:size(EFA,5)
for urb = 1:size(EFA,6)
Traff2.name(k) = {sprintf('%s/%s-%i/%i%%/%s',char(roads.RoadEnv(urb)),char(roads.RoadType(roa)),roads.RoadSpeeds(spd),roads.RoadGradient(dec),char(roads.Congestion(cog)))};
VEH = zeros(1,size(EFA,1));
VEHp = zeros(1,size(EFA,1));
for veh = 1:size(EFA,1)
if ~isnan(EFA(veh,roa,spd,dec,cog,urb))
VEH(veh) = VEH(veh)+1;
end
if EFA(veh,roa,spd,dec,cog,urb)>0
VEHp(veh) = VEHp(veh)+1;
end
end
nums1(k) = sum(VEH);
nums2(k) = sum(VEHp);
k = k+1;
end; end; end; end; end;
Traff2.NumVehEF_MODEL = nums1';
Traff2.NumVehEFz_MODEL = nums2';
errs = 0;
for i=1:height(Traff)
idx = find(ismember(Traff2.name,Traff.name(i)));
if (Traff.NumVehEF_HBEFA(i))>0 && (Traff2.NumVehEF_MODEL(idx))>0
elseif (Traff.NumVehEF_HBEFA(i))>0
errs = errs +1;
else
end
end
if errs >0
fprintf('#### SHEET CONVERISION ERROR ####\n')
else
fprintf('CONVERISION TESTED OK.\n')
end
Traff = [Traff,Traff2(:,2:end)];
fprintf('HBEFA : %i / %i / %i \n',round(min(Traff.NumVehEF_HBEFA)),round(mean(Traff.NumVehEF_HBEFA)),round(max(Traff.NumVehEF_HBEFA)) )
fprintf('MODEL : %i / %i / %i \n',round(min(Traff.NumVehEF_MODEL)),round(mean(Traff.NumVehEF_MODEL)),round(max(Traff.NumVehEF_MODEL)) )
Traff = Traff(Traff.NumVehEF_MODEL>0|Traff.NumVehEF_HBEFA>0,:);
writetable(Traff,ConVfile,'Sheet',sprintf('%s',char(comps(com))))
end
\ No newline at end of file
......@@ -73,8 +73,6 @@ LightVehiclesIdx = TM.ClassNum==1|TM.ClassNum==2;
HeavyVehiclesIdx = TM.ClassNum==3|TM.ClassNum==4;
BusesVehiclesIdx = TM.ClassNum==5|TM.ClassNum==6|TM.ClassNum==7;
VD = Vehicle_dist.Vdist;
%--------------------------------------------------------------
......
......@@ -25,6 +25,7 @@ function [RLinks] = PreProcess_Roads_Input_Traffic()
% Roads_Add_Width()
% Roads_Scale_Traffic_to_Year()
% Roads_Congestion_Parameters()
% Roads_Fix_Roadfields
%
% 20.09.2020 -Henrik Grythe
% Kjeller NILU
......@@ -61,6 +62,7 @@ if use_temporary_files
for i=1:length(fields)
fprintf('%s\n',char(fields(i)))
end
return
catch
fprintf('### No temporary file found ###\n')
prj = read_projection(traffile);
......@@ -92,7 +94,7 @@ end
% Test if required fields are there:
for i=1:length(required_fields)
if isfield(RLinks,required_fields(i))
if isfield(RLinks,required_fields(i)) || use_temporary_files
else
required_fields(i)
error('Missing required field (possibly others)')
......@@ -126,7 +128,9 @@ end
%--------------------------------------------------------------------------
RLinks = Roads_Remove_NoTrafficRoads(RLinks);
RLinks = Fix_Roadfields(RLinks);
RLinks = Roads_Fix_Roadfields(RLinks);
RLinks = Roads_Recalc_DISTANCE(RLinks);
%--------------------------------------------------------------------------
if isfield(RLinks,'KAPTID_06_') && ~isfield(RLinks,'KOTID_MORGEN')
......@@ -180,66 +184,3 @@ if use_temporary_files
end
end
%
% try
% prj = read_projection(traffile);
% fprintf('Reading large file: \n\t %s ...',fname)
% RLinks = shaperead(traffile);
% fprintf('done\n')
% fields = fieldnames(RLinks);
% fprintf('Fields:\n')
% for i=1:length(fields)
% fprintf('%s,',char(fields(i)))
% if rem(i,10)==0;fprintf('\n');end
% end
% fprintf('\n')
% catch
% RLinks=[];
% warning('Did not find road links at defined path')
% fprintf('continuing...\n')
% return
% end
%
% %
%
% a = strfind(traffile,'/');
% if isempty(a)
% fname = [];
% ofile = strcat(tfold,tfiles.RL);
% else
% fname = traffile(a(end)+1:end);
% ofile = strcat(tfold,tfiles.RL);
% end
% % Try to load a temporary file:
% try
% fprintf('Trying to read temp file\n%s_temp...\n',tfiles.RL)
% RLinks = shaperead(tfiles.RL);
% fprintf('Found Pre-Processed file:\n %s \n',tfiles.RL)
% fields = fieldnames(RLinks);
% fprintf('Fields:\n')
% for i=1:length(fields)
% fprintf('%s\n',char(fields(i)))
% end
% return
% catch
% fprintf('### No temporary file found ###\n')
% % try to read the road shape
% try
% prj = read_projection(traffile);
% fprintf('Reading large file: \n\t %s ...',fname)
% RLinks = shaperead(traffile);
% fprintf('done\n')
% fields = fieldnames(RLinks);
% fprintf('Fields:\n')
% for i=1:length(fields)
% fprintf('%s,',char(fields(i)))
% if rem(i,10)==0;fprintf('\n');end
% end
% fprintf('\n')
% catch
% RLinks=[];
% warning('Did not find road links at defined path')
% fprintf('continuing...\n')
% return
% end
% end
\ No newline at end of file
......@@ -27,7 +27,7 @@ function [Sn] = Roads_Add_Municipality(RLinks,Komm_shape)
fprintf('* call Add_ROAD_Municipality *\n')
fprintf('Adding kommune field based on geo file: \n%s\n',Komm_shape)
Ks = shaperead(Komm_shape);
prj = HEDGE_read_projection(Komm_shape);
prj = read_projection(Komm_shape);
% Adds a Roads Midpoint based on the starting point and ending point of
% Road
......
......@@ -15,6 +15,6 @@ else
tfiles.RL = sprintf('%sHEDGE_RL_%4i_%s_temp',tfold,Tyear,traffile);
end
ofiles.RLShape = sprintf('%sHEDGE_Emissions_%4i',ofold,Tyear);
ofiles.RLShape = sprintf('%sTraffic_Emissions_%4i',ofold,Tyear);
end
\ No newline at end of file
......@@ -5,7 +5,7 @@ function Vehicle_Distribution_Statistical_output(Vehicle_data,T)
% Tb.ModNum = T.ModelNumber;
% Tb.ModName = T.Name;
% Tb.Class = str2num(char(T.ClassNum));
%
%
% Tb.OsloNV = modelNV(1,:)';
% Tb.OsloADD = modelTD(1,:)';
% Tb.OsloPct = modelTDfrac(1,:)'*100;
......@@ -18,19 +18,29 @@ National2D = zeros(size(National,2)*size(National,3),size(National,1));
Ts = table;
for k = 1:size(National,1)
t = 1;
Ttemp = table;
for j = 1:size(National,3)
for i = 1:size(National,2)
Ts.Num(t) = t;
Ts.Euro(t) = Vehicle_data.D4_euro(j);
Ts.Kategori(t) = Vehicle_data.D3_nybiltype(i);
Ts.Var1(t) = National(k,i,j);
t=t+1;
if k==1
Ttemp = table;
Ttemp.Num = t;
Ttemp.Euro = Vehicle_data.D4_euro(j);
Ttemp.Kategori = Vehicle_data.D3_nybiltype(i);
Ttemp.Var1 = National(k,i,j);
Ts = [Ts;Ttemp];
t = t+1;
else
Ts.Var1(t) = National(k,i,j);
t = t+1;
end
end
end
vn = find(ismember(Ts.Properties.VariableNames,'Var1'));
Ts.Properties.VariableNames(vn)= {sprintf('x%i',Vehicle_data.D1_yrs(k))};
end
writetable(T,'National_SSB_Vehicle_Number_Distribution.xlsx','Sheet','NationalNumberVehicles')
file = 'National_SSB_Vehicle_Number_Distribution.xlsx';
writetable(T,file,'Sheet','NationalNumberVehicles')
fprintf('Wrote National Vehicle stats file :\n%s\n',file)
% National Driving Distance
National = squeeze(sum(Vehicle_data.YTD,2));
k=1;
......@@ -38,13 +48,21 @@ National2D = zeros(size(National,2)*size(National,3),size(National,1));
Ts2 = table;
for k = 1:size(National,1)
t = 1;
Ttemp = table;
for j = 1:size(National,3)
for i = 1:size(National,2)
Ts2.Num(t) = t;
Ts2.Euro(t) = Vehicle_data.D4_euro(j);
Ts2.Kategori(t) = Vehicle_data.D3_nybiltype(i);
Ts2.Var1(t) = National(k,i,j);
t=t+1;
if k==1
Ttemp = table;
Ttemp.Num = t;
Ttemp.Euro = Vehicle_data.D4_euro(j);
Ttemp.Kategori = Vehicle_data.D3_nybiltype(i);
Ttemp.Var1 = National(k,i,j);
Ts2 = [Ts2;Ttemp];
t=t+1;
else
Ts2.Var1(t) = National(k,i,j);
t = t+1;
end
end
end
vn = find(ismember(Ts2.Properties.VariableNames,'Var1'));
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment