Difference between revisions of "LISFLOOD-FP and MATLAB"
Line 71: | Line 71: | ||
end | end | ||
fclose('all'); | 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== | ==Making movies from .wd files== |
Revision as of 15:50, 18 March 2008
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);