Alternative way of coding

1 次查看（过去 30 天）
M.J. Glasbergen
M.J. Glasbergen 2012-3-27
I'm trying to create an application that, based on a simple assumption of water flow, results in flooded area's on a DTM. The principle works, but the way it is done now the tool is extremely slow. I think there should be an alternative way to speed things up using the matrix functionalities in Matlab a bit better, but I don't know how. Any idea anyone?
dbstop if error
[FileName,PathName] = uigetfile('*.asc;*.txt','Selecteer ASCII gridfile');
fid = fopen([PathName FileName], 'r');
Heading = cell(6,1);
for ii = 1:6
AsciiLine = fgetl(fid);
Heading{ii,1} = AsciiLine;
end
ii = 1;
AsciiLine = fgetl(fid);
while ischar(AsciiLine)
CellLine = textscan(AsciiLine,'%f');
Grid(ii,:) = CellLine{1,1}';
ii = ii + 1;
AsciiLine = fgetl(fid);
end
fclose(fid);
TestGrid = Grid;
NoData = -9999;
TestGrid(logical(TestGrid==NoData)) = NaN;
[y x] = ndgrid(1:size(TestGrid,1),1:size(TestGrid,2));
CellSize = char(Heading(5,1));
CellSize = str2double(CellSize(max(regexp(Heading{5,1}, '\s'))+1:end));
Rainfall = 19.8;
InitialDepth = (ones(size(TestGrid))*(Rainfall/1000)).*logical(~isnan(TestGrid));
TotalVolume = sum(sum(InitialDepth))*(CellSize^2);
Waterlevel = TestGrid+InitialDepth;
WaterlevelTemp = TestGrid+InitialDepth;
Stabiel =0;
LoopTeller = 0;
[WetCellsRow WetCellsColumn WetCellsVector] = find(InitialDepth);
while Stabiel == 0
LoopTeller = LoopTeller + 1;
WaterlevelChange = zeros(size(Waterlevel));
tic
for ii = 1:size(WetCellsVector,1)
if ~isnan(TestGrid(WetCellsRow(ii),WetCellsColumn(ii))) && (Waterlevel(WetCellsRow(ii),WetCellsColumn(ii))-TestGrid(WetCellsRow(ii),WetCellsColumn(ii))) > 0
out = min(max(abs(bsxfun(@minus,x,reshape(WetCellsColumn(ii),1,1,[]))),abs(bsxfun(@minus,y,reshape(WetCellsRow(ii),1,1,[])))),[],3);
dist = min(hypot(bsxfun(@minus,y,permute(WetCellsRow(ii),[3 2 1])),bsxfun(@minus,x,permute(WetCellsColumn(ii),[3 2 1]))),[],3);
MinSurroundLevel = Waterlevel(logical(out==1));
MinSurroundLevel = MinSurroundLevel(logical(((Waterlevel(out==1)-Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)))./dist(out==1))==min(((Waterlevel(out==1)-Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)))./dist(out==1)))));
if Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)) >= MinSurroundLevel(1,1)
MeanLevel = (sum(Waterlevel((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)) + Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)))/(size(Waterlevel((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1),1)+1);
if ((Waterlevel(WetCellsRow(ii),WetCellsColumn(ii))-TestGrid(WetCellsRow(ii),WetCellsColumn(ii)))/sum(sum(sum(logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1)))))+(Waterlevel((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)) > MeanLevel;
WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii)) = WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii))+(MeanLevel-Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)));
WaterlevelChange((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1) = WaterlevelChange(((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1))+(MeanLevel-Waterlevel(((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)));
else
WaterlevelChange((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1) = WaterlevelChange((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)+((Waterlevel(WetCellsRow(ii),WetCellsColumn(ii))-TestGrid(WetCellsRow(ii),WetCellsColumn(ii)))/sum(sum(sum(logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1)))));
WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii)) = WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii))+(TestGrid(WetCellsRow(ii),WetCellsColumn(ii)) - Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)));
end
end
end
end
Waterlevel = Waterlevel + WaterlevelChange;
Waterdepth = Waterlevel-TestGrid;
Waterdepth(isnan(Waterdepth)) = 0;
[WetCellsRow WetCellsColumn WetCellsVector] = find(Waterdepth);
TimerResult(LoopTeller) = toc;
if round(WaterlevelTemp(~isnan(Waterlevel))*1000) == round(Waterlevel(~isnan(Waterlevel))*1000)
Stabiel = 1;
CheckVolume = sum(sum(Waterdepth(~isnan(Waterdepth))))*(CellSize^2);
Waterlevel(isnan(Waterlevel)) = NoData;
Waterdepth(isnan(Waterdepth)) = NoData;
Files = ['Waterlevel' 'Waterdepth'];
for ii = 1:2
fid = fopen([PathName Files(ii) datestr(now,'yyyymmdd')], 'wt');
for jj = 1:6
fprintf(fid,'%s\r',char(Heading{jj,1}));
end
GridSize = size(Grid);
fprintf(fid,[repmat('%g\t',1,GridSize(2)-1) '%g\n'],Grid.');
fclose(fid);
end
fid = fopen([PathName 'TimerResult'], 'wt');
for ii = 1:size(TimerResult,1)
fprintf(fid,'%s\r',TimerResult(ii,1));
end
close(fid);
else
WaterlevelTemp = Waterlevel;
end
end

请先登录，再回答此问题。

回答（3 个）

Jan
Jan 2012-10-15
I'd start with cleaning the code. Examples:
  • You can omit the "logical" in logical(out==1), because "out==1" replies a logical already. But instead of creating this value again and again, it would be more efficient and much nicer to do this once only and store the rsult in a new variable.
  • The repeated calculation of "~isnan(Waterlevel)" wastes time also.
  • Files = ['Waterlevel' 'Waterdepth']; fid = fopen([PathName Files(ii)... will most likely not do what you expect, because "Files" is one string and "Files(ii)" is the ii.th character. You need Files = {'Waterlevel' 'Waterdepth'} with curly braces.
  • It is not a good idea to ask a forum to improve lines like these, when there is not even a single comment about the meaning:
dist = min(hypot(bsxfun(@minus,y,permute(WetCellsRow(ii),[3 2 1])), ...
bsxfun(@minus,x,permute(WetCellsColumn(ii),[3 2 1]))),[],3);
  • Finally I strongly suggest to use the PROFILEr to find out, which lines consumes the most time.
M.J. Glasbergen
M.J. Glasbergen 2012-10-15
Thnx Jan,
Unfortunately, most of the time is lost in the amount of repetitions and not in a single line. I know it helps trying to speed up the code by speeding up the individual lines, but in this case I was more or less lokking for a more efficient apporoach. It works it's way through a grid one cell at a time until it reaches a steady state. I think there should be a way to do that part more efficient. I'll see if I can find my m-file with comments for each line (I should have one, might just have to translate the comments). I've posted this question in March already (before I created a commented version) without recieving any replies and was now just trying to see if it would help if I changed the title of the question :-)
M.J. Glasbergen
M.J. Glasbergen 2012-11-22
This time with comments
if true
% Code stops on error while debugging
dbstop if error
%Select file [FileName,PathName] = uigetfile('*.asc;*.txt','Select ASCII gridfile'); fid = fopen([PathName FileName], 'r');
%Blanc array for headings Heading = cell(6,1);
%Get headings from file and fill heading array for ii = 1:6
AsciiLine = fgetl(fid);
Heading{ii,1} = AsciiLine;
end
ii = 1; AsciiLine = fgetl(fid);
%Converting lines with data to a 'grid' while ischar(AsciiLine)
CellLine = textscan(AsciiLine,'%f');
Grid(ii,:) = CellLine{1,1}';
ii = ii + 1;
AsciiLine = fgetl(fid);
end
%Close file fclose(fid);
%Define value for NoData in 'grid' and convert to NaN TestGrid = Grid; NoData = -9999; TestGrid(logical(TestGrid==NoData)) = NaN;
%Code for the replacement of the bwdist function from the image processing toolbox. %Thnx to Andrei Bobrov via MatlabCentral: %[i1,j1] = find(yourarray) %[y x] = ndgrid(1:size(yourarray,1),1:size(yourarray,2)) %out = min(max(abs(bsxfun(@minus,x,reshape(j1,1,1,[]))),abs(bsxfun(@minus,y,reshape(i1,1,1,[])))),[],3) [y x] = ndgrid(1:size(TestGrid,1),1:size(TestGrid,2));
%Determining the grid cell size from the headings CellSize = char(Heading(5,1)); CellSize = str2double(CellSize(max(regexp(Heading{5,1}, '\s'))+1:end));
%Rainfall Rainfall = 19.8; %mm
%Waterdepth on the grid when defined rain falls on the grid InitialDepth = (ones(size(TestGrid))*(Rainfall/1000)).*logical(~isnan(TestGrid));
%Watervolume (for water balance check) TotalVolume = sum(sum(InitialDepth))*(CellSize^2);
%Waterlevel in m + AD (grid level + waterdepth) Waterlevel = TestGrid+InitialDepth; WaterlevelTemp = TestGrid+InitialDepth;
Stabiel =0; LoopTeller = 0;
%Only look at the 'wet' cells [WetCellsRow WetCellsColumn WetCellsVector] = find(InitialDepth);
%As long as there is no stable answer, run the loop while Stabiel == 0
%Increase counter
LoopTeller = LoopTeller + 1;
%Blanc array for water level change
WaterlevelChange = zeros(size(Waterlevel));
%Start timer (for debugging only)
tic
%For all wet cells (each wet cell is evaluated individually):
for ii = 1:size(WetCellsVector,1)
%If there is a value in the current cell and there is water on the cell (waterlevel > gridlevel)
if ~isnan(TestGrid(WetCellsRow(ii),WetCellsColumn(ii))) && (Waterlevel(WetCellsRow(ii),WetCellsColumn(ii))-TestGrid(WetCellsRow(ii),WetCellsColumn(ii))) > 0
%Determine the distance from the current cell to the surrounding cells. (max) 8 surrounding cells het the value 1, (max) 16 around them the value 2 etc. Bepalen van de afstand van de huidige cel naar de omliggende
out = min(max(abs(bsxfun(@minus,x,reshape(WetCellsColumn(ii),1,1,[]))),abs(bsxfun(@minus,y,reshape(WetCellsRow(ii),1,1,[])))),[],3);
%Determine the distance from the middle of the current cell to the middle of the surrounding cells. (max) 4 directly attached cells (left, right, top and bottom) get the value 1, (max) 4 diagonally attached cells get the value SQRT(2) etc
dist = min(hypot(bsxfun(@minus,y,permute(WetCellsRow(ii),[3 2 1])),bsxfun(@minus,x,permute(WetCellsColumn(ii),[3 2 1]))),[],3);
%Determine the surrounding water levels and the largest hydraulic head between the current cell and the surrounding cells. Water will flow from the current cell to the one with the largest hydraulic head
MinSurroundLevel = Waterlevel(logical(out==1));
MinSurroundLevel = MinSurroundLevel(logical(((Waterlevel(out==1)-Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)))./dist(out==1))==min(((Waterlevel(out==1)-Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)))./dist(out==1)))));
%Water will only flow if one of the surrounding cells has a lower value
if Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)) >= MinSurroundLevel(1,1)
%Determine the average water level from the current cell and the lowest surrounding cell(s)
MeanLevel = (sum(Waterlevel((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)) + Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)))/(size(Waterlevel((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1),1)+1);
%If the resulting water level in the surrounding cells reaches a value that is higher than the average level if the water from the current cell is added...
if ((Waterlevel(WetCellsRow(ii),WetCellsColumn(ii))-TestGrid(WetCellsRow(ii),WetCellsColumn(ii)))/sum(sum(sum(logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1)))))+(Waterlevel((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)) > MeanLevel;
...then only the amount of water resulting in the average water level will be transported from the current cell. The movement of water is put in a temporary array
WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii)) = WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii))+(MeanLevel-Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)));
...then only the amount of water resulting in the average water level will be transported to the surrounding cell(s). The movement of water is put in a temporary array
WaterlevelChange((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1) = WaterlevelChange(((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1))+(MeanLevel-Waterlevel(((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)));
%Otherwise
else
...all the water from the current cell will go to the surrounding cell. The movement of water is put in a temporary array
WaterlevelChange((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1) = WaterlevelChange((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)+((Waterlevel(WetCellsRow(ii),WetCellsColumn(ii))-TestGrid(WetCellsRow(ii),WetCellsColumn(ii)))/sum(sum(sum(logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1)))));
WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii)) = WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii))+(TestGrid(WetCellsRow(ii),WetCellsColumn(ii)) - Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)));
end
end
end
end
%The new levels will be calculated when all water movements have been calculated
Waterlevel = Waterlevel + WaterlevelChange;
%The new waterdepths are calculated. NaN values will become 0
Waterdepth = Waterlevel-TestGrid;
Waterdepth(isnan(Waterdepth)) = 0;
%Determine the wet cells again
[WetCellsRow WetCellsColumn WetCellsVector] = find(Waterdepth);
%Stop the Timer and save the result to monitor the progress.
TimerResult(LoopTeller) = toc;
%If the resulting water level array is the same as the one from the previous timestep, the result is stable and new Ascii grids will be created
if round(WaterlevelTemp(~isnan(Waterlevel))*1000) == round(Waterlevel(~isnan(Waterlevel))*1000)
Stabiel = 1;
CheckVolume = sum(sum(Waterdepth(~isnan(Waterdepth))))*(CellSize^2);
Waterlevel(isnan(Waterlevel)) = NoData;
Waterdepth(isnan(Waterdepth)) = NoData;
Files = ['Waterlevel' 'Waterdepth'];
for ii = 1:2
fid = fopen([PathName Files(ii) datestr(now,'yyyymmdd')], 'wt');
for jj = 1:6
fprintf(fid,'%s\r',char(Heading{jj,1}));
end
GridSize = size(Grid);
fprintf(fid,[repmat('%g\t',1,GridSize(2)-1) '%g\n'],Grid.');
fclose(fid);
end
fid = fopen([PathName 'TimerResult'], 'wt');
for ii = 1:size(TimerResult,1)
fprintf(fid,'%s\r',TimerResult(ii,1));
end
close(fid);
else
WaterlevelTemp = Waterlevel;
% imagesc (Waterlevel); figure(gcf)
end
end end

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by