MATLAB Examples

Contents

The Geomorphic Flood Index: GFI=ln(hr/H)

Authors: Caterina Samela (caterina.samela@unibas.it); Salvatore Manfreda

Date: December 2016

Description: This script returns the Geomorphic Flood Index (Samela et al., 2016; Samela, C., 2015; Manfreda et al. 2015)

Required Input data:

   demvoid   m x n matrix with elevation values (raster grid Digital Elevation Model)
   demcon    Hydrologically conditioned or filled Digital Elevation Model
   FlowDir   Flow Direction
   FlowAcc   Flow accumulation
   G         Slope (matrix with tangent of the 8-connected neighborhood gradient of a Digital Elevation Model)
clc
clear all
close all

Loading Input data

To load the input data, we suggest the use of the function 'rasterread' from the 'TopoToolbox' by Schwanghart & Kuhn (2010) (https://topotoolbox.wordpress.com/).

% EXAMPLE
[demcon,X,Y,header] = rasterread('demcon_huc507_bigsandy_guyandotte.txt');
[demvoid,~,~,~] = rasterread('demvoid_huc507_bigsandy_guyandotte.txt');
[FlowDir,~,~,~] = rasterread('flowdir_huc507_bigsandy_guyandotte.txt');
[FlowAcc,~,~,~] = rasterread('flowacc_huc507_bigsandy_guyandotte.txt');
[G,~,~,~] = rasterread('Slope_huc507_bigsandy_guyandotte.txt');

% cellsize: horizontal spacing between data points
cellsize = abs(Y(2)-Y(1));

figure(1)
subplot(2,2,1)
imagesc(X(1,:),Y(:,1),demvoid)
axis image; axis xy; axis off; colorbar
title('DEM')
subplot(2,2,2)
imagesc(X(1,:),Y(:,1),FlowDir); axis image; axis xy; axis off; colorbar
title('Flow Direction')
subplot(2,2,3)
imagesc(X(1,:),Y(:,1),FlowAcc); axis image; axis xy; axis off; colorbar
title('Flow Accumulation')
subplot(2,2,4)
imagesc(X(1,:),Y(:,1),G); axis image; axis xy; axis off; colorbar
title('Slope')
% Create textbox
annotation(figure(1),'textbox',...
    [0.410974563953488 0.917492836553016 1 0.083756345177665],...
    'String',{'INPUT DATA'},...
    'FontWeight','bold',...
    'FontSize',14,...
    'LineStyle','none');

Identification of the Drainage Network

To define a channel matrix, a hydrologically conditioned DEM (or a filled DEM) must be used. The drainage network is identified on the basis of an area-slope criterion proposed by Giannoni et al. (2005). Specifically, the channel is assumed to start from locations where the quantity AS^k (where A is the upslope contributing area, S is the local slope) exceeds a given threshold, while its path to the outlet is identified by following the maximum slope direction. The exponent k takes on values within the range 1.2–2.7, with an average of 1.7 (see Giannoni et al., 2005). Here is assumed k = 1.7 , and the threshold is set equal to 10^5 m^2.

dem=demcon;

[a, b]= size(dem);
dem(isnan(dem))=100000;
CHANNEL=(zeros(size(dem)));
CHANNEL=(FlowAcc.*cellsize^2.*(G+0.0001).^1.7)>10^5;
CHANNEL=double(CHANNEL);


for i=1+1:a-1
    for j=1+1:b-1
        if CHANNEL(i,j)>0
            x=i; y=j;
            while dem(x,y) < 10000 && x<a-1 && x>1 && y<b-1 && y>1

                   switch FlowDir(x,y)
                        case 1
                            x=x; y=y+1;
                        case 128
                            x=x-1; y=y+1;
                        case 64
                            x=x-1; y=y;
                        case 32
                            x=x-1; y=y-1;
                        case 16
                            x=x; y=y-1;
                        case 8
                            x=x+1; y=y-1;
                        case 4
                            x=x+1; y=y;
                        case 2
                            x=x+1; y=y+1;
                   end
                   if(CHANNEL(x,y)==1)
                       break;
                   else
                        CHANNEL(x,y)=1;
                   end
            end

        end
    end
%     disp(i)
end

Calculation of:

1) the the contributing area in the nearest section of the drainage network hydrologically connected to the point under exam, here named "Ariver"[m^2]; 2) the Elevation difference to the nearest channel, H[m].

dem=demcon;
[a, b]= size(dem);
H=single(zeros(size(dem)));
Ariver = zeros(size(dem));

for i=1+1:a-1
    for j=1+1:b-1
        if dem(i,j)<10000
            x=i; y=j;
            Ld=0; dm=0;
            while CHANNEL(x,y) == 0 && x<a-1 && x>1 && y<b-1 && y>1 && FlowDir(x,y)>-9999
                switch FlowDir(x,y)
                        case 1
                            x=x; y=y+1;
                            Ld=Ld + cellsize;
                        case 128
                            x=x-1; y=y+1;
                            Ld=Ld + cellsize*2^0.5;
                        case 64
                            x=x-1; y=y;
                            Ld=Ld + cellsize;
                        case 32
                            x=x-1; y=y-1;
                            Ld=Ld + cellsize*2^0.5;
                        case 16
                            x=x; y=y-1;
                            Ld=Ld + cellsize;
                        case 8
                            x=x+1; y=y-1;
                            Ld=Ld + cellsize*2^0.5;
                        case 4
                            x=x+1; y=y;
                            Ld=Ld + cellsize;
                        case 2
                            x=x+1; y=y+1;
                            Ld=Ld + cellsize*2^0.5;
                end
            end
            Ariver(i,j)=FlowAcc(x,y);
            H(i,j)=dem(i,j)-dem(x,y);
         end
    end
%      disp(i)
end

H(isnan(demvoid))= NaN;
Ariver(i+1,:)=0;
Ariver(:,j+1)=0;
Ariver=single(Ariver);
H(H==0) = 0.00001; % Cells inside the channel

clear dem a b i j x y Ld dm

Estimation of the water level 'hr' in the nearest element of the drainage network.

hr is computed as a function of the contributing area (Ariver) in the nearest section of the drainage network hydrologically connected to the location under exam (‘r’ stands for ‘river’) using an hydraulic scaling relation (Leopold and Maddock, 1953): h=aA^n Values of the parameters a and n can be taken from literature or calibrated. Here, we provide the parameters calibrated for the Ohio River basin (Samela et al., 2015):

a = 0.1035;
n = 0.4057;

% hydraulic scaling relation (A[km^2])
hr = a.*(((Ariver+1).*cellsize^2)./10^6).^n;

Calculation of the Geomorphic Flood Index, GFI = ln(hr/H)

hronH=hr./H;
ln_hronH=real(log(hronH));

figure(2)
subplot(2,2,1);
colormap('default')
imagesc(X(1,:),Y(:,1),CHANNEL);
axis image; axis xy; axis off; colorbar
title('Channel')
subplot(2,2,2);
imagesc(X(1,:),Y(:,1),H);
axis image; axis xy; axis off; colorbar
title('Elevation difference, H')
subplot(2,2,3);
imagesc(X(1,:),Y(:,1),Ariver);
axis image; axis xy; axis off; colorbar
title('Contributing area, A_{river}')
subplot(2,2,4);
imagesc(X(1,:),Y(:,1),ln_hronH);
axis image;axis xy; axis off; colorbar
title('Geomorphic Flood Index, ln(h_r/H)')

Writing GFI matrix to an ESRI ArcGIS ASCII file

To save the Geomorphic Flood Index calculated, we suggest the use of the function 'rasterwrite', from the 'TopoToolbox' by Schwanghart & Kuhn (2010) (https://topotoolbox.wordpress.com/). Rasterwrite writes a matrix to an ESRI ArcGIS ASCII file.

rasterwrite('GFI.txt',X,Y,ln_hronH);

Reference

1. Giannoni, F., Roth, G., Rudari, R., 2005. A procedure for drainage network identification from geomorphology and its application to the prediction of the hydrologic response. Adv. Water Resour. 28, 567–581. 2. Manfreda, S., Samela, C., Gioia, A., Consoli, G., Iacobellis, V., Giuzio, L., Cantisani, A., & Sole, A. (2015). Flood-Prone Areas Assessment Using Linear Binary Classifiers based on flood maps obtained from 1D and 2D hydraulic models, Natural Hazards, 79(2), 735-754. 3. Manfreda, S., Samela, C., Sole, A., & Fiorentino, M. (2014b). Flood-prone Areas Assessment Using Linear Binary Classifiers based on Morphological Indices, ASCE-Vulnerability, Uncertainty, and Risk Quantification, Mitigation, and Management (pp. 2002-2011). 4. Samela, C. (2015). Flood hazard mapping over large regions using geomorphic procedures (Doctoral dissertation, Università degli Studi della Basilicata, Italy). 5. Samela, C., Manfreda, S., De Paola, F., Giugni, M., Sole, A., & Fiorentino, M. (2016). DEM-Based Approaches for the Delineation of Flood-Prone Areas in an Ungauged Basin in Africa, Journal of Hydrologic Engineering, 06015010. 6. Schwanghart, W., Kuhn, N.J. (2010): TopoToolbox: a set of Matlab functions for topographic analysis. Environmental Modelling & Software, 25, 770-781. [DOI: 10.1016/j.envsoft.2009.12.002].