Example_MatFlood_Reduction_Factor.m

%%% Apply reduction factor to the simulated flood depth %%%
% This code should be used after EXAMPLE_MATFLOOD_VARFWL.m

clear; close all; clc

% Change the path to where the script/ toolbox are
folder = pwd;

cd(folder)

% add the toolbox
addpath(genpath(folder))

%% Load results

load Flood_depth_Boston.mat

% colormaps (for plotting)
terrain     = load('TerrainColorMap','mycmap');

%% If observations are available, they can be used to determine the value of the reduction factor

% high water marks (flood observations)
hwms =[-71.04839425	42.32466284
-71.03133002	42.33080898
-71.03103399	42.34006521
-71.02385351	42.33505098
-71.04583567	42.30988461
-71.04711496	42.32096035
-71.04165948	42.34250886
-71.03392032	42.32994154
-71.04859513	42.31526909];

%% Plot topography

figure; pcolor(SF_s2.lon,SF_s2.lat,SF_s2.z); shading flat; colorbar
hold all; 

h1 = plot(SF_s2.wb(1),SF_s2.wb(2),'.m','MarkerSize',20);

ax = gca;
colormap(ax,terrain.mycmap)
clim([0 10])

dz = 0:1:30;
[~, h11]= contour(SF_s2.lon,SF_s2.lat,SF_s2.z,dz,'color','k','LineWidth',.1,'showtext','off');

h3 = plot(hwms(:,1),hwms(:,2),...
    'o','MarkerSize',5,'MarkerFaceColor','r','MarkerEdgeColor','k');

title('Elevation data (m), Boston city')
legend([h1,h3],'Point to indicate water','High water marks','Location','Best')

%% Section 1: univariate reduction factor (spatially uniform reduction factor)

% 1.1 Find the position of the high water marks on the grid ***************
% *************************************************************************
pos_hwms = nan(size(hwms,1),1);

for i = 1: size(hwms,1)

    [~,pos_hwms(i)] =...
        min(abs(distance(SF_s2.lat(:),SF_s2.lon(:),hwms(i,2),hwms(i,1))));

end

%% 1.2 Find the optimal values of reduction factor (rf) ********************
% *************************************************************************

% values of reduction factor to test
rf = linspace(0.5e-3,2e-3,10);

% store the number of high water marks within the flooded area for each rf
hwms_in_f = nan(length(rf),1);

for i = 1: length(rf)

    disp([num2str(i) ' out of ' num2str(length(rf))])

    % i reduction factor
    rf_i = rf(i);

    if i == 1
        % the first run takes more time, the next runs will use some of the
        % previous calculations (prevcalc) in order to save computational time
        
        tic
        RF_s1 = red_fac(SF_s2,rf_i);
        toc

        prevcalc = RF_s1.prevcalc;

    else % following runs use "prevcalc" to save computational time

        tic
        RF_s1 = red_fac(SF_s2,rf_i,prevcalc);
        toc
    end

    % find the number of high water marks within the flooded area
    hwms_in                 = RF_s1.FloodDepthRF(pos_hwms);
    hwms_in(isnan(hwms_in)) = [];
    hwms_in_f(i)            = length(hwms_in);
end

% Plot
figure;
plot(rf,hwms_in_f,'.-k','MarkerSize',20)
grid minor
ylabel('High water marks within flooded area (#)')
xlabel('Reduction factor')
set(gca,'FontName','Times','FontSize',12)

%% 1.3 Find the optimal value of the reduction factor as the maximum 
% reduction factor (highest reduced flooded area) that includes the most 
% high water marks 
% *************************************************************************
fopt    = find(hwms_in_f== max(hwms_in_f),1,'last');
rf_opt  = rf(fopt);

% obtain the flooded area for the optimal reduction factor
RF_s1 = red_fac(SF_s2,rf_opt,prevcalc);

% Plot
figure;
subplot(1,2,1); pcolor(SF_s2.lon,SF_s2.lat,SF_s2.FloodDepth); shading flat; colorbar
title('Flood before applying reduction factor')

subplot(1,2,2); pcolor(SF_s2.lon,SF_s2.lat,RF_s1.FloodDepthRF); shading flat; colorbar
title(['Applying reduction factor ' num2str(rf_opt)])



%% Section 2: using several reduction factors (spatially varying reduction factor)

% Set of reduction factors
rf = [-71.0514   42.3409    -0.00018
    -71.0475   42.3214    0.0012
    -71.0397   42.3494    0.002
    -71.0113   42.3154    0.0008
    -71.0107   42.3366    -0.0002];

% Apply spatially varying reduction factor
tic
[SF_red_s2,rfgrid] = red_fac(SF_s2,rf,prevcalc);
toc

% check how many high water marks are within the flood area
hwms_in_f_sv = SF_red_s2.FloodDepthRF(pos_hwms);

% Plot
figure; pcolor(SF_s2.lon,SF_s2.lat,SF_red_s2.FloodDepthRF); shading flat; colorbar
title('Flood depth reduced by spatially varying rf')
hold all
h3 = plot(hwms(:,1),hwms(:,2),...
    'o','MarkerSize',5,'MarkerFaceColor','r','MarkerEdgeColor','k');