LISFLOOD-FP and MATLAB
Revision as of 17:48, 18 February 2009 by Tim-fewtrell (talk | contribs)
This page contains various Matlab functions that may be useful for analysing results from the LISFLOOD-FP model.
File import and export
Importing ascii raster files
This function will import an ascii raster file using the string filename. To use it, copy and past text into a textfile called ascii_reader.m
function [dem, ncols, nrows, xllcorner, yllcorner, cellsize] = ascii_reader (filename)
% [DEM, ncols, nrows, xllcorner, yllcorner, cellsize] = ascii_reader (filename)
% This function reads arc .asc files
% requires filename string as input
% j neal
% 18/6/07
%% read an ascii file
fin = fopen(filename,'r');
A = fscanf(fin,'%s',1); ncols = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); nrows = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); xllcorner = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); yllcorner = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); cellsize = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); nodata = fscanf(fin,'%f',1); %#ok<NASGU>
dem = fscanf(fin,'%f',[ncols, nrows]);
dem = dem';
fclose('all');
You can view the imported file using the command:
imagesc(dem);
Exporting ascii raster files
This function will create an ascii raster file called filename from the 2D matrix DEM. To use it, copy and past text into a textfile called ascii_write.m
function ascii_write (filename, DEM, xllcorner, yllcorner, cellsize, nodata)
% ascii_write (filename, DEM, xllcorner, yllcorner, cellsize, nodata)
% only filename (string) and DEM (2D matrix) are essential
% this function writes ascii raster files for arc
% j neal
% 6/3/2008
if nargin < 2,
error('Requires at least two input arguments');
end
if nargin < 4,
xllcorner = 0;
yllcorner = 0;
end
if nargin < 5,
cellsize = 1;
end
if nargin < 6
nodata = -9999;
end
A = size(DEM);
fout = fopen (filename,'w');
fprintf(fout, 'ncols %.0f\n', A(2));
fprintf(fout, 'nrows %.0f\n', A(1));
fprintf(fout, 'xllcorner %f\n', xllcorner);
fprintf(fout, 'yllcorner %f\n', yllcorner);
fprintf(fout, 'cellsize %f\n', cellsize);
fprintf(fout, 'NODATA_value %f\n', nodata);
for ii = 1:A(1)
B = DEM(ii,:);
fprintf(fout, '%1.3f ', B);
fprintf(fout, '\n');
end
fclose('all');
Importing .profile files
This function imports .profile files
function [data, M] = profileread (resroot, numfiles, t, string1, string2)
% This function reads LISFLOOD-FP .profile files
%
% DATA = PROFILEREAD (RESROOT)
% where RESROOT is a string containing the resroot from the LISFLOOD-FP
% .par file. Absolute path names can also be used e.g. 'C:\LISFLOOD\resroot'
% DATA is a structure containing .header {cell aray}
% .segmentN (numerical array for segment N)
%
% DATA = PROFILEREAD (RESROOT, NUMFILES)
% NUMFILES = number of files to read. Reads m files assuming that the files
% are numbered such that resroot-000m.profile.
% DATA(m) is a structure containing m files
%
% DATA = PROFILEREAD (RESROOT, NUMFILES, 1)
% specifies the profiles were exported a user specified time such that
% resroot-000x-t.profile
%
% [DATA, M] = PROFILEREAD (RESROOT, NUMFILES, 1, string1, string2)
% instructs the function to plot string1 or sting1 against string2... these
% strings muct match the headernames in the .profile file e.g. 'Flow'
% returns movieframes in the structure M. Also plots the data.
%
% Note: LISFLOOD-FP will name multiple files 0 to m-1 but DATA(m) will count from 1 to m
%
% j.neal
% 12/5/2008
%%
if nargin < 1,
error('Requires filename');
end
if nargin < 2,
% only one file
numfiles = 0;
end
if nargin < 2,
t = 0;
end
if nargin < 4,
% no plotting required
plotting = 0;
M = 0;
elseif nargin < 5,
% use plotter
plotting = 1;
elseif nargin < 6,
plotting = 2;
end
%% Arrange data read
if t == 1
fileex = '-t.profile';
else
fileex = '.profile';
end
if numfiles == 0
% read in a single .profile file
filename = [resroot,fileex];
data = read_profile (filename);
elseif numfiles < 9999
% read a series of .profile files
for i = 1:numfiles
ii = i-1;
if ii < 10
resroot2 = [resroot,'-000',num2str(ii),fileex];
elseif ii < 100
resroot2 = [resroot,'-00',num2str(ii),fileex];
elseif ii < 1000
resroot2 = [resroot,'-0',num2str(ii),fileex];
else
resroot2 = [resroot,'-',num2str(ii),fileex];
end
%returnes structure with format data(i).segmentX... also includes .header
[data(i), seg_num] = read_profile (resroot2); %#ok<AGROW>
end
else
error ('Too many files or incorrect input to numfiles');
end
%% Plotting under development
% account for 0 seg num
seg_num = seg_num +1;
% if plotting == 0 % do nothing... no plots
if plotting == 1
% find the column you want to plot
I = strmatch(string1 , data(1).header);
% plot a single variable
A = size(data);
for i = 1:A(2)
figure(i);
for j = 1:seg_num
k = j-1;
subplot(1,seg_num,j);
eval(['plot(data(i).segment',num2str(k),'(:,',num2str(I),'));']);
tstring = ['Plot of ',resroot,num2str(i),' segment ',num2str(k)];
title(tstring);
% get movie frame
M(i,j) = getframe(gcf); %#ok<AGROW>
end
end
elseif plotting == 2
% find the columns you want to plot
I = strmatch(string1 , data(1).header);
I2 = strmatch(string2 , data(1).header);
A = size(data);
for i = 1:A(2)
figure(i);
for j = 1:seg_num
k = j-1;
subplot(1,seg_num,j);
eval(['plot(data(i).segment',num2str(k),'(:,',num2str(I),'),data(i).segment',num2str(k),'(:,',num2str(I2),'));']);
% get movie frame
tstring = ['Plot of ',resroot,num2str(i),' segment ',num2str(k)];
title(tstring);
M(i,j) = getframe(gcf); %#ok<AGROW>
end
end
end
%% read .profile file
function [data, seg_num] = read_profile (filename)
fin = fopen(filename,'r');
more_segs = 1;
% read 'Channel_segment:' from file
temp = fscanf(fin,'%s',1); %#ok<NASGU>
while more_segs == 1
% read filename
seg_num = fscanf(fin,'%f',1);
for i = 1:11
data.header{i} = fscanf(fin,'%s',1);
end
testline = 1;
j = 0;
while testline == 1
j = j+1;
% reads the first bit of the line as a string
A = fscanf(fin, '%s',1);
% this string is either a numer of a new channel sgment of the end
% of the file
% if end of file
ST = feof (fin);
if ST == 1
% reached the end of the file exit both while loops
testline = 0;
more_segs = 0;
else
TF = strcmp('Channel_segment:', A);
if TF == 1
% new segment has been found
testline = 0;
else
% read rest for line into data
B(j,1) = str2double(A); %#ok<AGROW>
for k = 2:11
B(j,k) = fscanf(fin, '%f',1); %#ok<AGROW>
end
end
end
end
eval(['data.segment',num2str(seg_num),' = B;']);
clear B
end
fclose('all');
Writing parameter files
This function can be used to write many LISFLOOD-FP parameter files with N different channel and floodplain roughness, using your own N-by-2 matrix called roughness. You will need to edit the code for your application by commenting unwanted lines with a percentage symbol.
function prm_writer (roughness)
% This function writes LISFLOOD-FP parameter files numbered 1 to n using
% the roughness values in the n-by-2 matrix roughness. yuo will need to
% edit the code for your application.
[num_sims] = size(roughness);
for i = 1:num_sims(1)
a = ['mymodel', num2str(i),'.par']; b = num2str(i);
fout = fopen(a,'w');
fprintf(fout,'DEMfile mydem.dem.ascii\n');
fprintf(fout,'resroot myres%s\n',b);
fprintf(fout,'dirroot mydir\n');
fprintf(fout,'sim_time 245000.0\n');
fprintf(fout,'initial_tstep 10.0\n');
fprintf(fout,'massint 300.0\n');
fprintf(fout,'saveint 10000.0\n');
%fprintf(fout,'overpass 135000.0\n');
fprintf(fout,'fpfric %1.3f\n',roughness(i,2));
fprintf(fout,'nch %1.3f\n',roughness(i,1));
%fprintf(fout,'manningfile mymanfile%s.n.ascii\n',b);
fprintf(fout,'riverfile myriver.river\n');
fprintf(fout,'bdyfile mybdy.bdy\n');
fprintf(fout,'stagefile mystage.stage\n');
fprintf(fout,'bcifile mybci.bci\n');
%fprintf(fout,'startfile mystart.old\n');
%fprintf(fout,'checkpoint 2\n');
%fprintf(fout,'checkfile mycheck.chkpnt\n');
fprintf(fout,'elevoff\n');
fprintf(fout,'qoutput\n');
fprintf(fout,'diffusive\n');
fclose('all');
end
Plotting flow vectors from .Qx and .Qy files
This function will plot flow vectors from a .Qx and .Qy over a dem using the matlab functions image and quiver
function lisflood_flow_plotter (qxfile, qyfile, demfile, north, east, south, west, scaling)
% This function plots flow vectors from LISFLOOD_FP Qx and QY files
%
% lisflood_flow_plotter (qxfile, qyfile, demfile);
%
% these input variables are strings containing the relative of full
% pathnames of the .Qx .Qy and ascii dem files to be plotted.
%
% lisflood_flow_plotter (qxfile, qyfile, demfile, N, E, S, W);
%
% N, E, W and S are optional they specify the northernmost, southeernmost
% ect cell to be displayed. this can be used to crop the domain.
%
% lisflood_flow_plotter (qxfile, qyfile, demfile, N, E, S, W, scaling);
%
% scaling is an optional term used to manually scale the vector lengthes
% drawn by quiver (the arrow plotting function).
%
% J Neal
% 20/02/2008
%% decide if crop is required and give imput error message
if nargin < 3,
error('Requires at least three input arguments');
end
if nargin < 7,
crop = 0;
else
crop = 1;
end
if nargin < 7,
scaling = 1;
end
% Note don't do this for large areas ?
%% Task 1: read ascii files
% read Qx file
[QX, ncolsqx, nrowsqx, xllcornerqx, yllcornerqx, cellsizeqx] = read_file (qxfile); %#ok<NASGU>
% read Qy file
[QY, ncolsqy, nrowsqy, xllcornerqy, yllcornerqy, cellsizeqy] = read_file (qyfile); %#ok<NASGU>
% read DEM
[DEM, ncolsdem, nrowsdem, xllcornerdem, yllcornerdem, cellsizedem] = read_file (demfile);
%% crop data if required
if crop == 1;
DEM = DEM(north:south,west:east);
QX = QX(north:south,west:east+1);
QY = QY(north:south+1,west:east);
xllcornerdem = xllcornerdem + west*cellsizedem - cellsizedem;
yllcornerdem = yllcornerdem + (nrowsdem - south) * cellsizedem;
[nrowsqx, ncolsqx] = size(QX); [nrowsqy, ncolsqy] = size(QY); [nrowsdem, ncolsdem] = size(DEM);
end
%% Taks 2: Plot dem and generate figure
figure1 = figure('PaperSize',[20.98 29.68],'Color',[1 1 1]);
colormap('gray');
% Create axes
axes1 = axes('Parent',figure1,'YTickLabel',{yllcornerdem + nrowsdem*cellsizedem ,yllcornerdem},...
'YTick',[0.5,nrowsdem+0.5],...
'YDir','reverse',...
'XTickLabel',{xllcornerdem,xllcornerdem+ncolsdem*cellsizedem},...
'XTick',[0.5, ncolsdem + 0.5],...
'Layer','top',...
'DataAspectRatio',[1 1 1]...
%,'Clim',[13 30]... % uncomment to set colourmap limits manually
);
box('on');
hold('all');
image(DEM,'Parent',axes1,'CDataMapping','scaled');
axis('image');
xlabel('Easting (m)');
ylabel('Northing (m)');
%% Task 3: calculate flow vectors
% pre allocate memory for loop... it will run faster this way
xlocs = zeros(ncolsqx*nrowsqx+ncolsqy*nrowsqy,1);
ylocs = xlocs;
U = xlocs;
V = xlocs;
%work out X and Y locations
for i = 1:nrowsqx
for j = 1:ncolsqx
xlocs((i-1)*ncolsqx+j,1) = j - 0.5;
ylocs((i-1)*ncolsqx+j,1) = i;
U((i-1)*ncolsqx+j,1) = QX(i,j);
end
end
for i = 1:nrowsqy
for j = 1:ncolsqy
xlocs(ncolsqx*nrowsqx + (i-1)*ncolsqy +j,1) = j;
ylocs(ncolsqx*nrowsqx + (i-1)*ncolsqy +j,1) = i - 0.5;
V(ncolsqx*nrowsqx+(i-1)*ncolsqy+j,1) = QY(i,j);
end
end
% work out number of nonzeros elements
num_flows = nnz(U+V);
xlocs2 = zeros(num_flows,1);
ylocs2 = xlocs2;
U2 = xlocs2;
V2 = xlocs2;
% remover zero flow locations
k = 1;
for i = 1:ncolsqx*nrowsqx+ncolsqy*nrowsqy
if (U(i,1) == 0) && (V(i,1) == 0)
% do nothing
else
xlocs2(k,1) = xlocs(i,1);
ylocs2(k,1) = ylocs(i,1);
U2(k,1) = U(i,1);
V2(k,1) = V(i,1);
k = k+1;
end
end
% check numbers are correct
if k == num_flows + 1
% all is ok
else
disp('Problem with flows');
end
%% Task 4: overlay flow vectors
disp('Plotting data... if this is taking a long time you may need fewer locations');
quiver (xlocs2,ylocs2,U2,V2, scaling);
% quiver (xloc2,yloc2,U2,V2,'Parent',figure1);
disp('Done');
%% function to read LISFLOOD ascii files using filename
function [OUT, ncols, nrows, xllcorner, yllcorner, cellsize] = read_file (filename)
% read an ascii file
fin = fopen(filename,'r');
A = fscanf(fin,'%s',1); ncols = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); nrows = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); xllcorner = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); yllcorner = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); cellsize = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); nodata = fscanf(fin,'%f',1); %#ok<NASGU>
OUT = fscanf(fin,'%f',[ncols, nrows]);
OUT = OUT';
fclose('all');
Making movies from .wd files
This function will convert a series of .wd .Qx or .Qy file into a .avi movie with the same filename
function LISFLOOD_mov (resroot,fileex,vartype,num_snaps,snapint,dem,demCS,depthCS,framesPS,movQ)
% LISFLOOD_mov generates movies of LISFLOOD-FP output
%
% LISFLOOD_mov (resroot,fileex,vartype,num_snaps,snapint,dem,demCS,depthCS,framesPS,movQ)
%
% where
% resroot is the LISFLOOD resroot (relative to m file or absolute)(string)
% fileex is the file extension '.wd' (can be chaned to plot WSE and flows)
% vartype is the name of the variable being plotted e.g. 'Depth'
% num_snaps is the number of snapshot times
% snapint is the time interval between each LISFLOOD-FP snapshot (seconds)
% dem is the name of the demfile (relative or absolute)
% demCS is the range of dem hights to be plotted e.g. demCS = [0,20];
% depthCS is the range of depth values to be plotted e.g. depthCS = [0,10];
% framesPS number of frames per second (each plot is a simgle frame)
% movQ is the movie quality between 1 and 100 (100 is best)
%
% 8/11/2007. J Neal.
%
% EXAMPLE:
% LISFLOOD_mov ('C:\Experiment1\res22','.wd','Water depth',24,10000,...
% 'C:\Experiment1\dem10m.dem.ascii',[0,60], [0,10],1,100);
%% Movie maker
% read dem filenam is filename and 1 specifies the 0.00 should not be NaN's
[DEM, ncols, nrows, xllcorner, yllcorner, cellsize] = read_file (dem); %#ok<NASGU>
% generate empty .avi file
% A = [resroot,'.avi'];
% mov = avifile(A);
% loop through snapshots
for i = 1:num_snaps
if i < 10
resfile = [resroot,'-000',num2str(i),fileex];
elseif i < 100
resfile = [resroot,'-00',num2str(i),fileex];
elseif i < 1000
resfile = [resroot,'-0',num2str(i),fileex];
else
resfile = [resroot,'-',num2str(i),fileex];
end
% read in deoth file for resfile i (2 specifies the 0.00 should be NaN
[DEPTH, ncols, nrows, xllcorner, yllcorner, cellsize] = read_file (resfile);
% plot data as figure 1
plotter (DEPTH, DEM, nrows, ncols, cellsize, xllcorner, yllcorner,demCS,depthCS,snapint,i,vartype);
% get movie frame
M(i) = getframe(gcf); %#ok<AGROW>
hold off;
% mov = addframe(mov,M(i));
% close the figure
close all
end
% % close .avi file
% mov = close(mov);
% generate empty .avi file
disp('Generating .avi movie file');
A = [resroot,fileex,'.avi'];
movie2avi(M,A,'fps',framesPS,'quality',movQ);
disp(A);
%% Read LISFLOOD ascii files using filename
function [OUT, ncols, nrows, xllcorner, yllcorner, cellsize] = read_file (filename)
% read an ascii file
fin = fopen(filename,'r');
A = fscanf(fin,'%s',1); ncols = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); nrows = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); xllcorner = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); yllcorner = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); cellsize = fscanf(fin,'%f',1); %#ok<NASGU>
A = fscanf(fin,'%s',1); nodata = fscanf(fin,'%f',1); %#ok<NASGU>
OUT = fscanf(fin,'%f',[ncols, nrows]);
OUT = OUT';
fclose('all');
% convert nodata to NaN;
for ii = 1:nrows
for jj = 1:ncols
if OUT(ii,jj) == nodata
OUT(ii,jj) = NaN;
elseif OUT(ii,jj) == -99.000
% lisflood uses -99.000 as nodat so will also remove these
OUT(ii,jj) = NaN;
else
% do nothing
end
end
end
%% Plotter data for movie slide
function plotter (DEPTH, DEM, nrows, ncols, cellsize, xllcorner, yllcorner,demCS,depthCS,snapint,i,vartype)
% number of colors in colormap
num_colors = 64;
% DEM and variable color map
cmap = [gray(num_colors); jet(num_colors)];
% scale DEM and DEPTH data to fit with new colormap
demCS2 = demCS(2)-demCS(1);
dem_scalar = num_colors/demCS2;
DEM2 = DEM*dem_scalar;
depthCS2 = depthCS(2)-depthCS(1);
depth_scalar = num_colors/depthCS2;
DEPTH2 = (DEPTH*depth_scalar)+num_colors+1;
%Plot the DEM
figure1 = figure('PaperSize',[20.98 29.68],'Color',[1 1 1]);
colormap(cmap);
axes1 = axes('Parent',figure1,'YTickLabel',{yllcorner+(nrows*cellsize),yllcorner+cellsize},...
'YTick', [1,nrows],...
'YDir','reverse',...
'XTickLabel',{xllcorner,xllcorner+(ncols*cellsize)-cellsize},...
'XTick',[1,ncols],...
'Layer','top');
box('on');
hold('all');
image(DEM2,'Parent',axes1);
axis('equal');axis('tight');
% Overlay the water depth (or other variable)
H = image(DEPTH2);
% code to make depth plot layer transparent
OP = ones(nrows,ncols);
for ii = 1:nrows
for jj = 1:ncols
if DEPTH(ii,jj) > 0
else
OP(ii,jj) = 0;
end
end
end
% make transparent. depth layer must be H
set(H,'AlphaData',OP);
% label colourbar
color_breaks = 5;
colorbreak = num_colors/color_breaks;
SCALE = zeros(color_breaks*2,1);
SCALE(1) = 0.0001;
for j = 1:color_breaks*2
SCALE(j+1) = colorbreak*j;
end
SCALE(color_breaks+1) = num_colors+1;
VIS_SCALE = zeros(color_breaks*2,1);
VIS_SCALE(1) = demCS(1);
for j = 1:color_breaks-1
VIS_SCALE(j+1) = ((demCS2/color_breaks)*j)+demCS(1);
end
VIS_SCALE(color_breaks+1) = depthCS(1);
for j = 1:color_breaks
VIS_SCALE(j+1+color_breaks) = ((depthCS2/color_breaks)*j)+depthCS(1);
end
colorbar('YTick',SCALE,'YTickLabel',VIS_SCALE);
% Label axis
xlabel('BNG easting (m)');
ylabel('BNG northing (m)');
A = [vartype,' over DEM at time ',num2str((i*snapint)),'s'];
title(A);