LISFLOOD-FP and MATLAB

This page contains various Matlab functions that may be useful for analysing results from the LISFLOOD-FP model.

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 A = fscanf(fin,'%s',1); nrows = fscanf(fin,'%f',1);          %#ok A = fscanf(fin,'%s',1); xllcorner = fscanf(fin,'%f',1);      %#ok A = fscanf(fin,'%s',1); yllcorner = fscanf(fin,'%f',1);      %#ok A = fscanf(fin,'%s',1); cellsize = fscanf(fin,'%f',1);       %#ok A = fscanf(fin,'%s',1); nodata = fscanf(fin,'%f',1);         %#ok 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 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 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 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 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 for k = 2:11 B(j,k) = fscanf(fin, '%f',1); %#ok end end end end eval(['data.segment',num2str(seg_num),' = B;']); clear B end fclose('all');

Importing .stage files
This function imports .stage files and allows for checkpointing restarts

function [metadata, stagedata] = stage_read (filename, num_locations, checkpoint) % This function reads lisflood .stage files % [metadata, stagedata] = stage_read (filename, num_locations) % % if checkpointing was used an additional argument is required to look for % checkpointed data % [metadata, stagedata] = stage_read (filename, num_locations, % checkpoint=1); % % The code will also check for repetition and remove duplicate lines % % J Neal % 26/6/09 if nargin < 2, error('Requires at least two input arguments'); end if nargin < 3, checkpoint = 0; end fin = fopen(filename, 'r'); header1 = fgets(fin); %#ok header2 = fgets(fin); %#ok header3 = fgets(fin); %#ok % reads the metadata metadata = fscanf (fin, '%f', [4, num_locations]); metadata = metadata'; header4 = fgets(fin); %#ok header5 = fgets(fin); %#ok<NASGU> header6 = fgets(fin); %#ok<NASGU> header7 = fgets(fin); %#ok<NASGU> if checkpoint == 0 % reads stage data stagedata = fscanf (fin, '%f', [num_locations+1,inf]); stagedata = stagedata'; else % scan for checkpointing lines line = 1; checks = 0; checkpoint = 0; while 1 tline = fgetl(fin); % this will return -1 at end of file if tline == -1 checks = checks+1; checkpoint(checks+1) = line; %#ok<AGROW> line = line-1; break; end TF = strcmp (tline, '####################################################### Checkpoint restart ########################################################'); if TF == 1 % this line is a checkpoint line checks = checks+1; checkpoint(checks+1) = line; %#ok<AGROW> end % the data could be read in at this point but this could be very % slow for large stage files line = line+1; end %now we know how may times it has checkpointed and where disp (['Number of checkpoints = ',num2str(checks)]); fclose('all'); % re open file and read header again! fin = fopen(filename, 'r'); header1 = fgets(fin); %#ok<NASGU> header2 = fgets(fin); %#ok<NASGU> header3 = fgets(fin); %#ok<NASGU> % reads the metadata metadata = fscanf (fin, '%f', [4, num_locations]); metadata = metadata'; header4 = fgets(fin); %#ok<NASGU> header5 = fgets(fin); %#ok<NASGU> header6 = fgets(fin); %#ok<NASGU> header7 = fgets(fin); %#ok<NASGU> % read stage data if checks == 2 % reads stage data as there are no checkpoints stagedata = fscanf (fin, '%f', [num_locations+1,line]); stagedata = stagedata'; else stagedata = []; for i = 2:checks tempstagedata = fscanf (fin, '%f', [num_locations+1,checkpoint(i-1)-checkpoint(i)-1]); tempstagedata = tempstagedata'; % this is not efficient but shjould only happen a few time stagedata = [stagdata;tempstagedata]; % read the checkpoint line tline = fgetl(fin); %#ok<NASGU> end % we now have lots of data that may overlap so we only need unique % rows stagedata = unique(stagedata,'rows'); end 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);