Commit 1dd8a530 authored by Henrik Grythe's avatar Henrik Grythe
Browse files

bugfix HBEFA

parent bdeabd53
......@@ -16,7 +16,7 @@
function Emission_Factor_Model_group_HBEFA()
%--------------------------------------------------------------------------
% The model classes have a static Emission Factor for each driving
% condition. This means that
% condition. This means that
% 22.10.2020 -Henrik Grythe
% Kjeller NILU
%--------------------------------------------------------------------------
......@@ -37,6 +37,7 @@ T = readtable(input.files.SSB_Vehicle_dist,'Sheet','HBEFA778toMODEL');
Tm = readtable(input.files.SSB_Vehicle_dist,'Sheet','MODEL');
fprintf('read Sheet : %s and Sheet %s\n','HBEFA778toMODEL','MODEL')
fprintf('Using weight %s \n',Vehicle_weight)
Tout = readtable(input.files.HBEFA_roads,'Sheet','Roads');
% ModelNumber HBEFA_Weight Uniform_Weight
......@@ -66,7 +67,7 @@ for com =1:length(comps)
% Tables are slow to work with and requires a lot of memory, but given
% the state of HBEFA, they offer good control of emission factors and
% what/where they are missing.
TFout = table;
TFout = Tout;
% Loop all model vehicles
for mod = 1:length(modelN)
% Make the subset of emissions/weights needed to calculate a model
......@@ -82,113 +83,55 @@ for com =1:length(comps)
Nef(i) = length(find(~isnan(EFsub(i,:,:,:,:,:,:))));
end
uef = unique(Nef);
for mr =1:length(uef)
fprintf('MODEL Vehicle: %3i %-42s Found #%i EF\n',modelN(mod),char(Tm.Name(modelN(mod))),uef(mr))
end
% If there are only one amount (1440) of Emission factors
if length(uef) == 1
fprintf('MODEL Vehicle: %3i %-42s Found #%i EF\n',modelN(mod),char(Tm.Name(modelN(mod))),uef)
Tout = 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)
Trout = table;
if ~isnan(EFsub(1,roa,spd,dec,cog,urb))
Trout.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)))};
%fprintf('%s/%s-%i/%i%%/%s\n',char(roads.RoadEnv(urb)),char(roads.RoadType(roa)),roads.RoadSpeeds(spd),roads.RoadGradient(dec),char(roads.Congestion(cog)))
Trout.Nr = k;
Trout.RoadNum = roa;
Trout.SpeedNum = spd;
Trout.DeclNum = dec;
Trout.CongNum = cog;
Trout.EnviNum = urb;
k = k+1;
Tout = [Tout;Trout];
end
end
end
end
end
% If there are only one amount (1440 or No) of Emission factors
% if length(uef) == 1
EFac = zeros(height(Tout),1);
Wght = zeros(height(Tout),1);
for veh = 1:height(Tsub)
switch Vehicle_weight
case 'Uniform'
VW = Tsub.Uniform_Weight(veh);
SW = sum(Tsub.Uniform_Weight);
case 'NERVE'
VW = Tsub.NERVE_Weight(veh);
SW = sum(Tsub.NERVE_Weight);
case 'HBEFA'
VW = Tsub.HBEFA_Weight(veh);
SW = sum(Tsub.HBEFA_Weight);
end
EFac = nan(height(Tout),1);
Wght = nan(height(Tout),1);
for veh = 1:height(Tsub)
switch Vehicle_weight
case 'Uniform'
fprintf('%i of %i %s Weight %3.2f %-30s of %3.2f ',veh,height(Tsub),Vehicle_weight,Tsub.NERVE_Weight(veh),char(Tsub.Name(veh)),sum(Tsub.Uniform_Weight))
case 'NERVE'
fprintf('%i of %i %s Weight %3.2f %-30s of %3.2f ',veh,height(Tsub),Vehicle_weight,Tsub.NERVE_Weight(veh),char(Tsub.Name(veh)),sum(Tsub.NERVE_Weight))
case 'HBEFA'
fprintf('%i of %i %s Weight %3.2f %-30s of %3.2f ',veh,height(Tsub),Vehicle_weight,Tsub.NERVE_Weight(veh),char(Tsub.Name(veh)),sum(Tsub.HBEFA_Weight))
end
for cond = 1:height(Tout)
if veh == 1
switch Vehicle_weight
case 'Uniform'
EFac(cond) = EFsub(veh,Tout.RoadNum(cond),Tout.SpeedNum(cond),Tout.DeclNum(cond),Tout.CongNum(cond),Tout.EnviNum(cond))*Tsub.Uniform_Weight(veh);
Wght(cond) = Tsub.Uniform_Weight(veh);
case 'NERVE'
EFac(cond) = EFsub(veh,Tout.RoadNum(cond),Tout.SpeedNum(cond),Tout.DeclNum(cond),Tout.CongNum(cond),Tout.EnviNum(cond))*Tsub.NERVE_Weight(veh);
Wght(cond) = Tsub.NERVE_Weight(veh);
case 'HBEFA'
EFac(cond) = EFsub(veh,Tout.RoadNum(cond),Tout.SpeedNum(cond),Tout.DeclNum(cond),Tout.CongNum(cond),Tout.EnviNum(cond))*Tsub.HBEFA_Weight(veh);
Wght(cond) = Tsub.HBEFA_Weight(veh);
end
else
switch Vehicle_weight
case 'Uniform'
EFac(cond) = EFac(cond) + EFsub(veh,Tout.RoadNum(cond),Tout.SpeedNum(cond),Tout.DeclNum(cond),Tout.CongNum(cond),Tout.EnviNum(cond))*Tsub.Uniform_Weight(veh);
Wght(cond) = Wght(cond)+ Tsub.Uniform_Weight(veh);
case 'NERVE'
EFac(cond) = EFac(cond) + EFsub(veh,Tout.RoadNum(cond),Tout.SpeedNum(cond),Tout.DeclNum(cond),Tout.CongNum(cond),Tout.EnviNum(cond))*Tsub.NERVE_Weight(veh);
Wght(cond) = Wght(cond)+ Tsub.NERVE_Weight(veh);
case 'HBEFA'
EFac(cond) = EFac(cond) + EFsub(veh,Tout.RoadNum(cond),Tout.SpeedNum(cond),Tout.DeclNum(cond),Tout.CongNum(cond),Tout.EnviNum(cond))*Tsub.HBEFA_Weight(veh);
Wght(cond) = Wght(cond)+ Tsub.HBEFA_Weight(veh);
end
end
end
fprintf('\n')
fprintf('%i of %i %s Weight %3.2f %-30s of %3.2f \n',veh,height(Tsub),Vehicle_weight,VW,char(Tsub.Name(veh)),SW)
for cond = 1:height(Tout)
EFcond(cond,1) = EFsub(veh,Tout.RoadNum(cond),Tout.SpeedNum(cond),Tout.DeclNum(cond),Tout.CongNum(cond),Tout.EnviNum(cond));
end
Tout.EFac = EFac./Wght;
pos = find(ismember(Tout.Properties.VariableNames,{'EFac'}));
Tout.Properties.VariableNames(pos) = Tm.Name(modelN(mod));
if debug_mode
fprintf('%s_%s_weight',char(comps(com)),Vehicle_weight)
fprintf(' %7.1f/%7.1f/%7.1f (mean/max/min) ',mean(EFac./Wght),max(EFac./Wght),min(EFac./Wght))
fprintf(' %7.1f/%7.1f/%7.1f (mean/max/min) \n',mean(Wght),max(Wght),min(Wght))
else
fprintf('\n')
end
else
% if this is the case, there are NaNs in the emission factor
% and we will get wrong results.
for i=1:length(uef)
fprintf('\n\n\n\n\n\n\n\n\n\n#### Found #%i EF\n',uef(i))
if ~isnan(sum(EFcond))
fprintf('Found 1440 emission factors for vehicle\n')
EFac = EFac + EFcond*VW;
Wght = Wght + VW;
elseif ~isnan(nansum(EFcond))
fprintf('Found some emission factors for vehicle\n')
elseif isnan(nansum(EFcond))
fprintf('Found no emission factors for vehicle\n')
end
end
nanmean(EFac)
Tout.EFac = EFac;
pos = find(ismember(Tout.Properties.VariableNames,{'EFac'}));
Tout.Properties.VariableNames(pos) = Tm.Name(modelN(mod));
TFout = [TFout,Tout(:,pos)];
if mod == 1
% Add the first part of the table as well
TFout = Tout;
if debug_mode
fprintf(' %s_%s_weight',char(comps(com)),Vehicle_weight)
fprintf(' %7.1f/%7.1f/%7.1f (mean/max/min) ',mean(EFac./Wght),max(EFac./Wght),min(EFac./Wght))
fprintf(' %7.1f/%7.1f/%7.1f (mean/max/min) \n',mean(Wght),max(Wght),min(Wght))
else
% Concatonate only the column of the new vehicle data if the
% align.
[a,b,c] = intersect(TFout.Name,Tout.Name);
addCol = find(ismember(Tout.Properties.VariableNames,Tm.Name(modelN(mod))));
% Not 100% robust test!!!!?
if length(b)==length(c)
TFout = [TFout,Tout(:,addCol)];
else
fprintf('\n\n\n\n\n\n\n\n\n\n#### Found #%i & %i EF\n',a,b)
end
fprintf('\n')
end
end
EFrdCond = table2array(TFout(:,8:end));
......
......@@ -178,8 +178,8 @@ idl = find(contains(Tout.Properties.VariableNames,'Light'));
idh = find(contains(Tout.Properties.VariableNames,'Heavy'));
idb = find(contains(Tout.Properties.VariableNames,'Buses'));
KDD.TraffDataL = table2array(Tout(:,idl));
KDD.TraffDataH = table2array(Tout(:,idl));
KDD.TraffDataB = table2array(Tout(:,idl));
KDD.TraffDataH = table2array(Tout(:,idh));
KDD.TraffDataB = table2array(Tout(:,idb));
KDD.kommNamesL = Tout.Properties.VariableNames(idl);
KDD.kommNamesH = Tout.Properties.VariableNames(idh);
KDD.kommNamesB = Tout.Properties.VariableNames(idb);
......@@ -189,4 +189,4 @@ save(tfiles.DD_Municipal,'KDD','Tout')
fprintf('Finished processing driving distances\n\n')
% DrivingDistances_per_RoadType = Tout;
% save(ofiles.MatlabOutput,'DrivingDistances_per_RoadType','-append')
end
\ No newline at end of file
end
......@@ -53,25 +53,14 @@ uKomm = unique(KommunerS);
% Congestion weights. Traffic part delayed (TPD)
% Combine the different congestion levels by volume to determine the
% actual EF:
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];
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];
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];
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.2, 0.2, 0.2, 0.4, 0.0];
TPD_L = TPD;
TPD_H = TPD;
TPD_B = TPD;
% EXTRACT
L = extractfield(RLinks,sprintf('L_ADT%04i',Tyear)); % Traffic Volume (# day-1)
......
......@@ -2,26 +2,26 @@ function NERVE_Emissions_Statistical_Oputput()
% Data needed by the MijlDirektoratet sheet
comps = [{'N2O'}];
comps = [{'CO2'},{'FC'},{'CH4'},{'N2O'}];
comps = [{'CO2'},{'FC'},{'CH4'}];
%comps = [{'CO2'},{'FC'}];
%comps = [{'CH4'}];
%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
% comps = [{'CO2'},{'FC'},{'CH4'}];
% comps = [{'CO2'},{'FC'}];
% comps = [{'CH4'}];
% 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
% 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
if ispc
bp = 'N:\Inby';
......@@ -47,7 +47,8 @@ for com = 1:length(comps)
% files that must be read for each year and component!
EF_File = sprintf('%sTemp/EF_On_AllRoadCond_Municipality_%i_%s.mat',pathn,Tyear,char(comps(com)));
load(EF_File)
fileout = sprintf('%sOutput2/NERVE_output_%s_%04i.mat',pathn,char(comps(com)),Tyear);
try
R_EF_File = sprintf('%sTemp/EFA_Table_MODEL_%s_Bio%i.mat',pathn,char(comps(com)),Tyear);
load(R_EF_File,'TFout'); % TFout
......@@ -60,14 +61,14 @@ for com = 1:length(comps)
% files that must be read per year
munFile = sprintf('%sTemp/Municipal_Traffic_Exchange_%i.mat',pathn,Tyear);
RdDistFile = sprintf('%sOutput/RoadTypeDistanceMunicipal%i.mat',pathn,Tyear);
RdDistFile = sprintf('%sOutput2/RoadTypeDistanceMunicipal%i.mat',pathn,Tyear);
SSBCPfile = sprintf('%sTemp/SSB_CarPark_%i.mat',pathn,Tyear);
load(munFile) % TrafficIN TrafficFROM kmne
load(RdDistFile) % KDD
load(SSBCPfile) % Vehicle_dist
fprintf('Reading large RoadLink file...')
RLinks = shaperead(sprintf('%sOutput/Traffic_Emissions_%i',pathn,Tyear));
RLinks = shaperead(sprintf('%sOutput2/Traffic_Emissions_%i',pathn,Tyear));
fprintf('Done\n')
DaysInYear = datenum([Tyear+1 1 1 0 0 0])-datenum([Tyear 1 1 0 0 0]);
......@@ -111,7 +112,7 @@ for com = 1:length(comps)
% EM_INr(k) = sum(EM(f).*IDO(f))+sum(EM(f2).*(1-IDO(f2)));
end
pdata =1e-9*[sum(L_EM),sum(H_EM),sum(B_EM),sum(L_EM+H_EM+B_EM)];
fprintf('\tLIGHT:%6.2f, HEAV:Y%6.2f, BUS:%6.2f, TOT:%6.2f\n',pdata)
fprintf('\tLIGHT:%6.2f, HEAVY:%6.2f, BUS:%6.2f, TOT: %8.3f\n',pdata)
% Vehicle_dist
%------------------------------------------------------------------
......@@ -168,11 +169,11 @@ for com = 1:length(comps)
end
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)
fprintf('\t\t\t\t\t\tTOT_IN : %8.3f diff %4.2f%%\n',pdata)
% END IN
%--------------------------------------------------------------------------
% This method assigns the total Emissions based on the exchange.
% END IN
%------------------------------------------------------------------
% This method assigns the total Emissions based on the exchange.
for komm = 1:size(TrafficFROM,2)
I = find(TrafficFROM(komm,:)>0);
......@@ -189,8 +190,8 @@ for com = 1:length(comps)
end
end
%--------------------------------------------------------------------------
% create *L_FROM(kom)* based on:::: L_IN and exchange
%------------------------------------------------------------------
% 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);
......@@ -203,12 +204,13 @@ for com = 1:length(comps)
end
end
end
fileout = sprintf('%sOutput/NERVE_output_%s_%04i.mat',pathn,char(comps(com)),Tyear);
pdata =[1e-9*nansum(nansum(EM_FROM)),100*nansum(nansum(EM_FROM))/sum(L_EM+H_EM+B_EM)];
fprintf('\t\t\t\t\t\tTOT_FROM: %8.3f diff %4.2f%%\n',pdata)
fprintf('Processed Emissions for %s year %i\n',char(comps(com)),Tyear)
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
addpath('/storage/nilu/Inby/Emission_Group/Emission_Models/HEDGE/MiljodirektoratetData')
NERVE_Miljodirektoratet_Output_generator()
end
\ No newline at end of file
......@@ -39,32 +39,34 @@ Ks = shaperead(Komm_shape);
% Adds a Roads Midpoint based on the starting point and ending point of
% Road
if ~isfield(RLinks,'MIDPOINT_X') || ~isfield(RLinks,'MIDPOINT_Y')
for i = 1:length(RLinks)
x = RLinks(i).X;
y = RLinks(i).Y;
MIDPOINT_X(i)=(x(1)+x(end-1))/2;
MIDPOINT_Y(i)=(y(1)+y(end-1))/2;
end
else
MIDPOINT_X = extractfield(RLinks,'MIDPOINT_X');
MIDPOINT_Y = extractfield(RLinks,'MIDPOINT_Y');
end
% if ~isfield(RLinks,'MIDPOINT_X') || ~isfield(RLinks,'MIDPOINT_Y')
% for i = 1:length(RLinks)
% x = RLinks(i).X;
% y = RLinks(i).Y;
% MIDPOINT_X(i)=(x(1)+x(end-1))/2;
% MIDPOINT_Y(i)=(y(1)+y(end-1))/2;
% end
% else
% MIDPOINT_X = extractfield(RLinks,'MIDPOINT_X');
% MIDPOINT_Y = extractfield(RLinks,'MIDPOINT_Y');
% end
for i = 1:length(RLinks)
x = RLinks(i).X;
y = RLinks(i).Y;
ix = ~isnan(x);
iy = ~isnan(x);
ix = find(~isnan(x));
iy = find(~isnan(y));
START_X(i) = x(ix(1));
START_Y(i) = y(iy(1));
END_X(i) = x(ix(1));
END_Y(i) = y(iy(1));
END_X(i) = x(ix(end));
END_Y(i) = y(iy(end));
MIDPOINT_X(i) = mean(x(ix));
MIDPOINT_Y(i) = mean(y(iy));
end
KOMMS = nan(size(START_X));
KOMME = nan(size(START_X));
KOMM = nan(size(START_X));
KOMM = nan(size(START_X));
for i = 1:length(Ks)
% FIND STARTING MUNICIPALITY
inS = inpolygon(START_X,START_Y,Ks(i).X,Ks(i).Y);
......
Markdown is supported
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