LISFLOOD-FP and MATLAB

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');

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

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);