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

Bugfix simplification of complex road placement

parent 1dd8a530
function NERVE_Emissions_Statistical_Oputput()
% Data needed by the MijlDirektoratet sheet
comps = [{'N2O'}];
% comps = [{'N2O'}];
comps = [{'CO2'},{'FC'},{'CH4'},{'N2O'}];
% comps = [{'CO2'},{'FC'},{'CH4'}];
% comps = [{'CO2'},{'FC'}];
......@@ -29,40 +29,55 @@ else
bp = '/storage/nilu/Inby';
end
pathn = sprintf('%s/Emission_Group/Emission_Models/HEDGE/',bp);
tfold = sprintf('%sTemp/' ,pathn);
ofold = sprintf('%sOutput2/',pathn);
ifold = sprintf('%sInput/' ,pathn);
pfold = sprintf('%sPlots/' ,pathn);
use_ido = 0;
makeplots = 1;
% Files that must be read once only
% 'Model_Vehicles_Merge_SSB_and_HBEFA_Vehicles.xlsx'
TM = readtable(sprintf('%sInput/Model_Vehicles_Merge_SSB_and_HBEFA_Vehicles.xlsx',pathn),'Sheet','MODEL');
model_class_file = sprintf('%sModel_Vehicles_Merge_SSB_and_HBEFA_Vehicles.xlsx',ifold);
TM = readtable(model_class_file,'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;
PCVehiclesIdx = TM.ClassNum==1;
LCVVehiclesIdx = TM.ClassNum==2;
BusesVehiclesIdx = TM.ClassNum==3|TM.ClassNum==4;
ATVehiclesIdx = TM.ClassNum==6|TM.ClassNum==7;
RTVehiclesIdx = TM.ClassNum==5;
for com = 1:length(comps)
% files that must be read per species
fprintf('%s\n',char(comps(com)))
fprintf('\n\n%s\n',char(comps(com)))
for Tyear = 2009:2019
% 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)));
EF_File = sprintf('%sEF_On_AllRoadCond_Municipality_%i_%s.mat',tfold,Tyear,char(comps(com)));
load(EF_File)
fileout = sprintf('%sOutput2/NERVE_output_%s_%04i.mat',pathn,char(comps(com)),Tyear);
fileout = sprintf('%sNERVE_output_%s_%04i.mat',ofold,char(comps(com)),Tyear);
try
R_EF_File = sprintf('%sTemp/EFA_Table_MODEL_%s_Bio%i.mat',pathn,char(comps(com)),Tyear);
try % this should only work for CO2
R_EF_File = sprintf('%sEFA_Table_MODEL_%s_Bio%i.mat',tfold,char(comps(com)),Tyear);
load(R_EF_File,'TFout'); % TFout
fprintf('Found Bio File! ')
catch
R_EF_File = sprintf('%sTemp/EFA_Table_MODEL_%s.mat',pathn,char(comps(com)));
fprintf('Found Bio File! %s',R_EF_File)
catch % all non CO2 species
R_EF_File = sprintf('%sEFA_Table_MODEL_%s.mat',tfold,char(comps(com)));
load(R_EF_File,'TFout'); % TFout
fprintf('No Bio File! %s',R_EF_File)
end
fprintf('%i\n',Tyear)
fprintf('\n%i\n',Tyear)
% files that must be read per year
munFile = sprintf('%sTemp/Municipal_Traffic_Exchange_%i.mat',pathn,Tyear);
RdDistFile = sprintf('%sOutput2/RoadTypeDistanceMunicipal%i.mat',pathn,Tyear);
SSBCPfile = sprintf('%sTemp/SSB_CarPark_%i.mat',pathn,Tyear);
munFile = sprintf('%sMunicipal_Traffic_Exchange_%i.mat',tfold,Tyear);
SSBCPfile = sprintf('%sSSB_CarPark_%i.mat',tfold,Tyear);
RdDistFile = sprintf('%sRoadTypeDistanceMunicipal%i.mat',ofold,Tyear);
load(munFile) % TrafficIN TrafficFROM kmne
load(RdDistFile) % KDD
load(SSBCPfile) % Vehicle_dist
......@@ -75,6 +90,7 @@ for com = 1:length(comps)
KOMM = extractfield(RLinks,'KOMMS');
KOMMe = extractfield(RLinks,'KOMME');
KOMMm = extractfield(RLinks,'KOMM');
IDO = extractfield(RLinks,'IDO');
L_ADT = extractfield(RLinks,sprintf('L_ADT%i',Tyear));
H_ADT = extractfield(RLinks,sprintf('H_ADT%i',Tyear));
......@@ -92,12 +108,23 @@ for com = 1:length(comps)
% find *L_IN(kom)*
ukomm = unique(KOMM);
for k=1:length(ukomm)
if use_ido
f = find(KOMM == ukomm(k));
f2 = find((KOMMe == ukomm(k))&(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)));
else
f = find(KOMMm == ukomm(k));
L_IN(k) = sum(LW(f));
H_IN(k) = sum(HW(f));
B_IN(k) = sum(BW(f));
end
end
fprintf('Trafikkarbeid Road : Oppsummert \n')
fprintf('Light %12.0f : %12.0f %4.2f%%\n',sum(LW),sum(L_IN),100*sum(LW)/sum(L_IN))
fprintf('Heavy %12.0f : %12.0f %4.2f%%\n',sum(HW),sum(H_IN),100*sum(HW)/sum(H_IN))
fprintf('Buses %12.0f : %12.0f %4.2f%%\n',sum(BW),sum(B_IN),100*sum(BW)/sum(B_IN))
L_EM = extractfield(RLinks,sprintf('EM_L_%s',char(comps(com))));
H_EM = extractfield(RLinks,sprintf('EM_H_%s',char(comps(com))));
......@@ -105,15 +132,18 @@ for com = 1:length(comps)
EM = L_EM+H_EM+B_EM;
for k=1:length(ukomm)
if use_ido
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)));
% 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)));
else
f = find(KOMMm == ukomm(k));
EM_INr(k) = sum(EM(f));
end
end
fprintf('Utslipp Road : Oppsummert \n----------------------\n')
pdata =1e-9*[sum(L_EM),sum(H_EM),sum(B_EM),sum(L_EM+H_EM+B_EM)];
fprintf('\tLIGHT:%6.2f, HEAVY:%6.2f, BUS:%6.2f, TOT: %8.3f\n',pdata)
fprintf('Roads LIGHT:%6.2f, HEAVY:%6.2f, BUS:%6.2f, TOT RDs : %8.3f\n',pdata)
% Vehicle_dist
%------------------------------------------------------------------
% find *NV(kom,veh)* from:: SSB data
......@@ -156,25 +186,27 @@ for com = 1:length(comps)
% Estimate *EM_IN(kom,Veh) * based on:::: EF_IN and exchange
for k=1:length(ukomm)
Ltcomp = sum(ORdComp(k,LightVehiclesIdx));
Htcomp = sum(ORdComp(k,HeavyVehiclesIdx));
Btcomp = sum(ORdComp(k,BusesVehiclesIdx));
for veh = 1:length(vehicles)
type = TM.Model_Class(veh);
if type <= 2
EM_IN(k,veh) = L_IN(k)*ORdComp(k,veh)*EF_IN(k,veh);
EM_IN(k,veh) = L_IN(k)*ORdComp(k,veh)/Ltcomp*EF_IN(k,veh);
elseif (type == 5 || type == 6 || type == 7)
EM_IN(k,veh) = H_IN(k)*ORdComp(k,veh)*EF_IN(k,veh);
EM_IN(k,veh) = H_IN(k)*ORdComp(k,veh)/Htcomp*EF_IN(k,veh);
elseif (type == 3 || type == 4)
EM_IN(k,veh) = B_IN(k)*ORdComp(k,veh)*EF_IN(k,veh);
EM_IN(k,veh) = B_IN(k)*ORdComp(k,veh)/Btcomp*EF_IN(k,veh);
end
end
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\t\tTOT_IN : %8.3f diff %4.2f%%\n',pdata)
pdata =[1e-9*nansum(nansum(EM_IN(:,LightVehiclesIdx ))),1e-9*nansum(nansum(EM_IN(:,HeavyVehiclesIdx ))),1e-9*nansum(nansum(EM_IN(:,BusesVehiclesIdx ))),1e-9*nansum(nansum(EM_IN)),100*nansum(nansum(EM_IN))/sum(L_EM+H_EM+B_EM)];
fprintf('NERVE LIGHT:%6.2f, HEAVY:%6.2f, BUS:%6.2f, TOT_IN : %8.3f TOT:diff %4.2f%%\n',pdata)
% 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);
......@@ -189,7 +221,6 @@ for com = 1:length(comps)
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)
......@@ -205,10 +236,196 @@ for com = 1:length(comps)
end
end
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)
pdata =[1e-9*nansum(nansum(EM_FROM(:,LightVehiclesIdx ))),1e-9*nansum(nansum(EM_FROM(:,HeavyVehiclesIdx ))),1e-9*nansum(nansum(EM_FROM(:,BusesVehiclesIdx ))),1e-9*nansum(nansum(EM_FROM)),100*nansum(nansum(EM_FROM))/sum(L_EM+H_EM+B_EM)];
fprintf('NERVE LIGHT:%6.2f, HEAVY:%6.2f, BUS:%6.2f, TOT_FROM: %8.3f TOT:diff %4.2f%%\n',pdata)
if makeplots
f = figure('visible','off');
scatter(L_IN,L_FROM,sum(NV,2).^.7,ukomm,'filled')
set(gca,'XScale','log','YScale','log');
cbh = colorbar;
cbh.Label.String = 'Kommunenummer';
hold on
plot([min([min(L_IN),min(L_FROM)]),max([max(L_IN),max(L_FROM)])],[min([min(L_IN),min(L_FROM)]),max([max(L_IN),max(L_FROM)])],'k--','LineWidth',3)
grid on
ylabel('Distanse kjørt med biler med opphave i kommune x')
xlabel('Distanse kjørt på veiene i kommune x')
tn = sprintf('Trafikkarbeid Alle Biler og kommuner %i',Tyear );
title(tn)
x= gcf;
x.Position=[1986 1 2495 1361];
tn = regexprep(tn,' ','_');
saveas(gcf,strcat(pfold,tn,'IN_FROM.jpg'))
f = figure('visible','off');
scatter(L_FROM,sum(TD,2),sum(NV,2).^.7,ukomm,'filled')
set(gca,'XScale','log','YScale','log');
cbh = colorbar;
cbh.Label.String = 'Kommunenummer';
hold on
plot([min([min(L_IN),min(L_FROM)]),max([max(L_IN),max(L_FROM)])],[min([min(L_IN),min(L_FROM)]),max([max(L_IN),max(L_FROM)])],'k--','LineWidth',3)
grid on
xlabel('Distanse kjørt med biler med opphave i kommune x')
ylabel('Distanse kjørt med biler (SSB)registrert i kommune x')
tn = sprintf('Trafikkarbeid Alle Biler og kommuner %i',Tyear );
title(tn)
x= gcf;
x.Position=[1986 1 2495 1361];
tn = regexprep(tn,' ','_');
saveas(gcf,strcat(pfold,tn,'FROM_SSBDD.jpg'))
f = figure('visible','off');
scatter((L_IN+H_IN+B_IN)',nansum(EM_IN,2)./(L_IN+H_IN+B_IN)',sum(NV,2).^.7,ukomm,'filled')
set(gca,'XScale','log');
cbh = colorbar;
cbh.Label.String = 'Kommunenummer';
grid on
xlabel('Distanse kjørt med biler med opphave i kommune x')
ylabel('Utslippsfaktor (g/km)')
tn = sprintf('Trafikkarbeid og Utslippsfaktor Alle Biler og kommuner %i %s',Tyear,,char(comps(com)));
hold on
scatter((L_IN+H_IN+B_IN)',nansum(EM_IN,2)./(L_IN+H_IN+B_IN)',5,'w','filled')
title(tn)
x= gcf;
x.Position=[1986 1 2495 1361];
tn = regexprep(tn,' ','_');
saveas(gcf,strcat(pfold,tn,'.jpg'))
f = figure('visible','off');
scatter((L_IN)',nansum(EM_IN(:,LightVehiclesIdx),2)./(L_IN)',sum(NV,2).^.7,ukomm,'filled')
set(gca,'XScale','log');
cbh = colorbar;
cbh.Label.String = 'Kommunenummer';
grid on
xlabel('Distanse kjørt med biler med opphave i kommune x')
ylabel('Utslippsfaktor (g/km)')
tn = sprintf('Trafikkarbeid og Utslippsfaktor Lette Biler %i %s',Tyear,char(comps(com)));
hold on
scatter((L_IN)',nansum(EM_IN(:,LightVehiclesIdx),2)./(L_IN)',5,'w','filled')
title(tn)
x= gcf;
x.Position=[1986 1 2495 1361];
tn = regexprep(tn,' ','_');
saveas(gcf,strcat(pfold,tn,'.jpg'))
f = figure('visible','off');
scatter((H_IN)',nansum(EM_IN(:,HeavyVehiclesIdx),2)./(H_IN)',sum(NV,2).^.7,ukomm,'filled')
set(gca,'XScale','log');
cbh = colorbar;
cbh.Label.String = 'Kommunenummer';
grid on
xlabel('Distanse kjørt med biler med opphave i kommune x')
ylabel('Utslippsfaktor (g/km)')
tn = sprintf('Trafikkarbeid og Utslippsfaktor Tunge Kjøretøy %i %s',Tyear,char(comps(com)));
hold on
scatter((H_IN)',nansum(EM_IN(:,HeavyVehiclesIdx),2)./(H_IN)',5,'w','filled')
title(tn)
x= gcf;
x.Position=[1986 1 2495 1361];
tn = regexprep(tn,' ','_');
saveas(gcf,strcat(pfold,tn,'.jpg'))
f = figure('visible','off');
scatter((B_IN)',nansum(EM_IN(:,BusesVehiclesIdx),2)./(B_IN)',sum(NV,2).^.7,ukomm,'filled')
set(gca,'XScale','log');
cbh = colorbar;
cbh.Label.String = 'Kommunenummer';
grid on
xlabel('Distanse kjørt med biler med opphave i kommune x')
ylabel('Utslippsfaktor (g/km)')
tn = sprintf('Trafikkarbeid og Utslippsfaktor Busser %i %s',Tyear,char(comps(com)));
hold on
scatter((B_IN)',nansum(EM_IN(:,BusesVehiclesIdx),2)./(B_IN)',5,'w','filled')
title(tn)
x= gcf;
x.Position=[1986 1 2495 1361];
tn = regexprep(tn,' ','_');
saveas(gcf,strcat(pfold,tn,'.jpg'))
close all
f = figure('visible','off');
bh1 = boxplot(EF_FROM(:,PCVehiclesIdx),'positions',(1:sum(PCVehiclesIdx))-.2,'symbol','r.','Color',[.4 .4 .4]);
hold on
bh2 = boxplot(EF_IN (:,PCVehiclesIdx),'positions',(1:sum(PCVehiclesIdx))+.2,'symbol','b.');
set(gca,'XTick',1:sum(PCVehiclesIdx),'XTickLabel',TM.Name(PCVehiclesIdx),'XTickLabelRotation',45)
grid on
tn = sprintf('Utslippsfaktor (EF IN) Personbil %i %s',Tyear,char(comps(com)) );
title(tn)
x= gcf;
x.Position=[1986 1 2495 1361];
text('Position',[10.0 0],'Interpreter','latex','String','EM IN','Color','b','FontSize',26);
text('Position',[20.0 0],'Interpreter','latex','String','EM FROM','Color',[.4 .4 .4],'FontSize',26);
tn = regexprep(tn,' ','_');
saveas(gcf,strcat(pfold,tn,'.jpg'))
f = figure('visible','off');
bh1 = boxplot(EF_FROM(:,LCVVehiclesIdx),'positions',(1:sum(LCVVehiclesIdx))-.2,'symbol','r.','Color',[.4 .4 .4]);
hold on
bh2 = boxplot(EF_IN (:,LCVVehiclesIdx),'positions',(1:sum(LCVVehiclesIdx))+.2,'symbol','b.');
set(gca,'XTick',1:sum(LCVVehiclesIdx),'XTickLabel',TM.Name(LCVVehiclesIdx),'XTickLabelRotation',45)
grid on
tn = sprintf('Utslippsfaktor (EF IN) Varebil %i %s',Tyear,char(comps(com)));
title(tn)
x= gcf;
x.Position=[1986 1 2495 1361];
text('Position',[10.0 0],'Interpreter','latex','String','EM IN','Color','b','FontSize',26);
text('Position',[20.0 0],'Interpreter','latex','String','EM FROM','Color',[.4 .4 .4],'FontSize',26);
tn = regexprep(tn,' ','_');
saveas(gcf,strcat(pfold,tn,'.jpg'))
f = figure('visible','off');
bh1 = boxplot(EF_FROM(:,BusesVehiclesIdx),'positions',(1:sum(BusesVehiclesIdx))-.2,'symbol','r.','Color',[.4 .4 .4]);
hold on
bh2 = boxplot(EF_IN (:,BusesVehiclesIdx),'positions',(1:sum(BusesVehiclesIdx))+.2,'symbol','b.');
set(gca,'XTick',1:sum(BusesVehiclesIdx),'XTickLabel',TM.Name(BusesVehiclesIdx),'XTickLabelRotation',45)
grid on
tn = sprintf('Utslippsfaktor (EF IN) Busser %s %i',char(comps(com)),Tyear );
title(tn)
x= gcf;
x.Position=[1986 1 2495 1361];
text('Position',[10.0 0],'Interpreter','latex','String','EM IN','Color','b','FontSize',26);
text('Position',[20.0 0],'Interpreter','latex','String','EM FROM','Color',[.4 .4 .4],'FontSize',26);
tn = regexprep(tn,' ','_');
saveas(gcf,strcat(pfold,tn,'.jpg'))
f = figure('visible','off');
bh1 = boxplot(EF_FROM(:,RTVehiclesIdx),'positions',(1:sum(RTVehiclesIdx))-.2,'symbol','r.','Color',[.4 .4 .4]);
hold on
bh2 = boxplot(EF_IN (:,RTVehiclesIdx),'positions',(1:sum(RTVehiclesIdx))+.2,'symbol','b.');
set(gca,'XTick',1:sum(RTVehiclesIdx),'XTickLabel',TM.Name(RTVehiclesIdx),'XTickLabelRotation',45)
tn = sprintf('Utslippsfaktor (EF IN) Lastebiler %s %i',char(comps(com)),Tyear );
grid on
title(tn)
x= gcf;
x.Position=[1986 1 2495 1361];
text('Position',[10.0 0],'Interpreter','latex','String','EM IN','Color','b','FontSize',26);
text('Position',[20.0 0],'Interpreter','latex','String','EM FROM','Color',[.4 .4 .4],'FontSize',26);
tn = regexprep(tn,' ','_');
saveas(gcf,strcat(pfold,tn,'.jpg'))
f = figure('visible','off');
bh1 = boxplot(EF_FROM(:,ATVehiclesIdx),'positions',(1:sum(ATVehiclesIdx))-.2,'symbol','r.','Color',[.4 .4 .4]);
hold on
bh2 = boxplot(EF_IN (:,ATVehiclesIdx),'positions',(1:sum(ATVehiclesIdx))+.2,'symbol','b.');
set(gca,'XTick',1:sum(ATVehiclesIdx),'XTickLabel',TM.Name(ATVehiclesIdx),'XTickLabelRotation',45)
grid on
tn = sprintf('Utslippsfaktor (EF IN) Trekkbiler %s %i',char(comps(com)),Tyear );
title(tn)
x= gcf;
x.Position=[1986 1 2495 1361];
text('Position',[10.0 0],'Interpreter','latex','String','EM IN','Color','b','FontSize',26);
text('Position',[20.0 0],'Interpreter','latex','String','EM FROM','Color',[.4 .4 .4],'FontSize',26);
tn = regexprep(tn,' ','_');
saveas(gcf,strcat(pfold,tn,'.jpg'))
close all
end
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')
close all
end
end
addpath('/storage/nilu/Inby/Emission_Group/Emission_Models/HEDGE/MiljodirektoratetData')
......
......@@ -25,6 +25,7 @@ startkommune = extractfield(RLinks,'KOMMS');
endekommune = extractfield(RLinks,'KOMME');
midtkommune = extractfield(RLinks,'KOMM');
komVeg = nan(size(startkommune));
komVeg(startkommune == endekommune&(startkommune == midtkommune))=1;
komVeg(startkommune ~= endekommune) = 0;
......@@ -36,18 +37,84 @@ fprintf('Roads shared between municipalities %7i \n',length(find(komVeg==0)))
fprintf('Roads without municipality %7i \n',length(find(isnan(komVeg))))
fprintf('----------------------------------------\n')
Keep_RLinks = RLinks(find(komVeg));
t_RLinks = RLinks(komVeg==0);
Komms = extractfield(Ks,'KOMMUNENUM');
Fraction = zeros(size(t_RLinks));
% First deal with roads that are continouus.
for i=1:length(t_RLinks)
startkommune = extractfield(t_RLinks(i),'KOMMS');
endekommune = extractfield(t_RLinks(i),'KOMME');
midtkommune = extractfield(t_RLinks(i),'KOMM');
xs = extractfield(t_RLinks(i),'X');
individual(i) = length(find(isnan(xs)));
subsegments(i) = length(find(~isnan(xs)))-1;
end
idx = individual==1 & subsegments==1;
Straight_single_roads = find(idx);
Curvy_single_roads = find(individual==1 & ~idx);
Complex_roads = find(individual>1 & subsegments~=1);
fprintf('----------------------------------------\n')
fprintf('Straight_single_roads municipalities %7i \n',length(Straight_single_roads))
fprintf('Curvy_single_roads municipalities %7i \n',length(Curvy_single_roads))
fprintf('Complex_roads municipality %7i \n',length(Complex_roads))
fprintf('Unexplained roads municipality %7i \n',length(find(komVeg==0))-length(Complex_roads)-length(Curvy_single_roads)-length(Straight_single_roads) )
fprintf('----------------------------------------\n')
RDcomplexity(Straight_single_roads) = 1;
RDcomplexity(Curvy_single_roads) = 2;
RDcomplexity(Complex_roads) = 3;
% keep a backup ht = t_RLinks;
% trim some of the subnodes off the roads a bit so they are not too long,
% as the test can run out of memory.
for i = 1:length(t_RLinks)
if(RDcomplexity(i)==2 && subsegments(i)>500)
xs = extractfield(t_RLinks(i),'X');
ys = extractfield(t_RLinks(i),'Y');
lxy = length(xs);
startx = xs(1);
starty = ys(1);
stopx = xs(end-1);
stopy = ys(end-1);
while lxy > 450
segD(1) = sqrt((xs(1)-xs(2))^2 +(ys(1)-ys(2))^2);
for r = 2:length(xs)
segD(r) = sqrt((xs(r-1)-xs(r))^2 +(ys(r-1)-ys(r))^2);
end
idx = find(segD <= min(segD)+1e-0);
segD(idx) = NaN;
idx2 = ~isnan(segD);
xs = xs(idx2);
ys = ys(idx2);
lxy = length(xs);
clear segD
end
newx = [startx,xs,stopx,NaN];
newy = [starty,ys,stopy,NaN];
t_RLinks(i).X = newx;
t_RLinks(i).Y = newy;
end
end
fail = 1;
Fraction = zeros(size(t_RLinks));
for i=1:length(t_RLinks)
startkommune = extractfield(t_RLinks(i),'KOMMS');
endekommune = extractfield(t_RLinks(i),'KOMME');
midtkommune = extractfield(t_RLinks(i),'KOMM');
fprintf('%-4i',i)
fprintf('%-5i',i)
if ~isnan(startkommune)
pla = find(Komms == startkommune);
for p=1:length(pla)
......@@ -71,6 +138,7 @@ for i=1:length(t_RLinks)
fail = fail+1;
end
end
if rem(i,50)==0; fprintf('\n'); end
end
KeepT = struct2table(Keep_RLinks);
t_T = struct2table(t_RLinks);
......
......@@ -63,6 +63,8 @@ fprintf('Light=%4.1f%%\n',100*(mean(Scale.ToYearL)-1))
fprintf('Heavy=%4.1f%%\n',100*(mean(Scale.ToYearH)-1))
fprintf('Buses=%4.1f%%\n',100*(mean(Scale.ToYearB)-1))
% Tabelize Roads
Tlink = struct2table(RLinks);
......
......@@ -23,20 +23,15 @@ for s=1:length(S)
in = inpolygon(S(s).X,S(s).Y,E.X,E.Y);
nn = sum(isnan(S(s).X));
DST(s) = S(s).DISTANCE;
if sum(in) == 0
% No vertices inside the domain: Road is outside and will be cut
if sum(in) == 0 % No vertices inside the domain: Road is outside and will be cut
fra = 0;
elseif sum(in) == length(S(s).X)-nn
% All vertices inside the domain Road is inside and will be used
elseif sum(in) == length(S(s).X)-nn % All vertices inside the domain Road is inside and will be used
fra = 1;
else % There are intersections, partly in domain. Calculate the part inside the domain.
xt = S(s).X(in);
yt = S(s).Y(in);
% if i==23817; stop; end
% NEED A TEST FOR SIZE AS IT CAN EASILY RUN OUT OF MEMORY DOING
% THIS TEST
% This does not necessarily find all the intersection
% This test sometimes fail to find all the intersections
if length(S(s).X) < 500
[xi,yi] = polyxpoly(S(s).X,S(s).Y,E.X,E.Y);
else
......@@ -46,28 +41,30 @@ for s=1:length(S)
[xi,yi] = polyxpoly(S(s).X([1,end-1]),S(s).Y([1,end-1]),E.X,E.Y);
end
if ~isempty(xi)
if length(xi)==1
if length(xt)==1
if ~isempty(xi) % At least one intersection with the kommune border is found.
if length(xi)==1 % Exactly one intersection is found
if length(xt)==1 % Its only one point inside.
tx = [xt,xi,NaN];
ty = [yt,yi,NaN];
nd = sqrt((xt-xi)^2+(yt-yi)^2);
else
else % There are several points inside domain.
% find the nearesty point to the intersection
clear d1
for nx = 1:length(xt)
d1(nx) = sqrt((xi-xt(nx))^2 +(yi-yt(nx))^2);
end
itp = find(d1==min(d1));
% Replace the closest point inside municipality by the
% intersection point. For distance calculations.
xt(itp) = xi;
yt(itp) = yi;
for nx = 1:length(xt)-1
d(nx) = sqrt((xt(nx)-xt(nx+1))^2 +(yt(nx)-yt(nx+1))^2);
end
tx = [xt,NaN];
ty = [yt,NaN];
for nx = 1:length(tx)-2
d(nx) = sqrt((tx(nx)-tx(nx+1))^2 +(ty(nx)-ty(nx+1))^2);
end
nd = sum(d); clear d
nd = nansum(d);
clear d
end
else
% multiple intersections:
......@@ -78,13 +75,13 @@ for s=1:length(S)
nd = sqrt((min(xst)-max(xst)).^2 +(min(yst)-max(yst)).^2);
end
else
warning('No intersection found for RL:\n')
nd = S(s).DISTANCE;
else % Did not find intersection cant deal with this right now
warning(sprintf('\nNo intersection found for RL: Border?\n'))
nd = (S(s).DISTANCE*1000);
end
fra = nd/S(s).DISTANCE;
if fra >1.01
fra = nd/(S(s).DISTANCE*1000);
if fra >1.01
fra = nd/(S(s).DISTANCE)
end
if debug_mode
fprintf('%04.4f / %04.4f (%04.1f%%) in StartGrid \n',DST(s)*fra,DST(s),fra*100)
......
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