Next Article in Journal
Two Dimensional Model for Backwater Geomorphology: Darby Creek, PA
Next Article in Special Issue
Geospatial Modelling of Watershed Peak Flood Discharge in Selangor, Malaysia
Previous Article in Journal
Identification of Hydrological Models for Enhanced Ensemble Reservoir Inflow Forecasting in a Large Complex Prairie Watershed
Previous Article in Special Issue
Hydrodynamic-Based Numerical Assessment of Flood Risk of Tamuín City, Mexico, by Tampaón River: A Forecast Considering Climate Change
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Flood-Prone Area Assessment Using GIS-Based Multi-Criteria Analysis: A Case Study in Davao Oriental, Philippines

1
Institute of Computing and Engineering, Davao Oriental State College of Science and Technology, Guang-guang, Dahican, Mati City 8200, Philippines
2
Graduate School for International Development and Cooperation, Hiroshima University, 1-5-1 Kagamiyama, Higashi-Hiroshima 739-8529, Japan
*
Author to whom correspondence should be addressed.
Water 2019, 11(11), 2203; https://doi.org/10.3390/w11112203
Submission received: 24 August 2019 / Revised: 15 October 2019 / Accepted: 19 October 2019 / Published: 23 October 2019
(This article belongs to the Special Issue Flood Modelling: Regional Flood Estimation and GIS Based Techniques)

Abstract

:
Flooding is one of the major destructive natural disasters in Davao Oriental, Philippines, and results primarily from a high incidence of typhoons and heavy rainfalls. The main objective of this study was to identify flood-prone risk areas by mapping them based on the integration of multiple indicators, including rainfall, slope, elevation, drainage density, soil type, distance to the main channel and population density. For this purpose, a GIS-based flood risk spatial assessment was conducted by using analytic hierarchy process (AHP), weights by rank (WR) and ratio weighting (RW) frameworks to determine the relative importance of each indicator against another in the province of Davao Oriental. The resulting flood-prone areas by the three methods are validated by comparing with the estimated flood map based on ground truthing points from a field survey. The comparison results show that AHP is the most appropriate method among them to assess flood hazard. The result of the AHP flood risk map shows that 95.99% (5451.27 km2) of Davao Oriental is under low and moderate flood risk. The high and very high flood risk area covers approximately 3.39% (192.52 km2) of the province, primarily in the coastal areas. Thirty-one out of the one hundred eighty-three (31/183) barangays (towns) are at a high to very high risk of flooding at current climate, calling for the immediate attention of decision-makers to develop mitigation strategies for the future occurrence of flooding in Davao Oriental.

1. Introduction

Flooding is considered one of the most serious and widespread natural hazards due to its devastating effects that endanger lives and cause property damage in the affected areas [1,2,3,4,5,6]. The changing climate behavior of extreme rainfalls and typhoons are the primary contributors to this problem. Human activities such as urbanization and the growth of settlements and assets in flooding areas likewise contribute to the increasing impacts of floods [1,2,3,4].
Phenomena commonly contributed to flooding are the limited capacity of river channels [7], human settlements in low-lying areas [8] and rapid growth of human settlements without upgrading the drainage infrastructures [1]. The problems related to flooding have increased extensively, and there is a need for efficient modeling to help mitigate the worst effects of flood disasters [9]. Thus, effective modeling will help in proper flood risk management planning and provides various insights into addressing the hazard and disaster problems. It is necessary to establish and implement a systematic approach to define areas that may have disastrous flood events.
Due to the prevailing occurrences of natural disasters that regularly hit the Philippines, the government created a law to reinforce preparation for possible calamities. The Republic Act 10121, also known as An Act Strengthening the Philippine Disaster Risk Reduction and Management (DRRM) System, was implemented, providing the national disaster risk reduction and management framework and institutionalizing the national disaster risk reduction and management plan [10]. Despite the enactment of the DRRM law, the local government units were not ready to face natural calamities such as typhoons and heavy rainfalls that cause massive flooding. Thus, effective flood risk assessment and strategies should be established, especially at the barangay (local term for town) level since each barangay area has its own physical and social characteristics.
Even though flooding is a severe hazard in Davao Oriental due to the typhoon, heavy rainfalls and storms, insufficient attention has been paid to flood hazard assessment. Recent scientific work undertaken in the province of Davao Oriental by Lagmay et al. [11] was concentrated on flood mapping in major river basins, which is a significant research gap identified by this study. They applied only in river run-off, and their study did not provide further insight into the entire province. Another study on the vulnerability assessment after Typhoon Pablo was conducted by Ross et al. [12] using light detection and ranging (LiDAR) as the primary tool. Ross et al. [12] used a high-resolution digital elevation model (DEM) to generate two-dimensional flood simulations of the floodplains using FLO-2D. Ross et al. [12] determined that 12%–28% of each of the eight municipalities is susceptible to flooding: mostly residential, business and community centers. It was also revealed that a storm surge with a height of approximately 1 m reached 500 m inland from the coastline during Typhoon Pablo. Unfortunately, this study focused only on the east coast municipalities, and Typhoon Pablo affected the entire province.
All the flooding-related studies in Davao Oriental are beneficial to providing a spatial presentation of the distribution of the flooded areas. However, no studies have yet to undertake the evaluation and flood risk mapping of the entire province of Davao Oriental using multiple datasets. Thus, the objective of this study is to identify flood-prone risk areas and map them based on several factors that are relevant to the study area. For this purpose, a GIS-based process for the spatial assessment of flood hazard was conducted by using the concepts of the analytic hierarchy process (AHP), weights by rank (WR) and ratio weighting (RW) frameworks with the multi-criteria indicators including rainfall, slope, elevation, drainage density, soil type, distance to the main channel and population density. The methodological approach used in this paper was inspired by various case studies [1,2,6,13,14,15,16] that used AHP and other GIS-based quantitative approaches. These previous works supported the study of risk of flooding as a combination of different factors. Moreover, adding more factors such as socio-economic variables (gender, income, educational attainment and housing) into the vulnerability map and applying the multi-criteria decision analysis (MCDA) methods can provide more insights in vulnerability assessments.
In the following, literature reviews on flood risk assessment methodology are summarized in Section 2, and the study area and datasets used in the analysis are described in Section 3. In Section 4, the MCDA methodology by AHP, WR and RW used in this study is described. The analysis results are presented in Section 5, followed by discussions on the sensitivity of flood hazard map to varying consistency ratio in Section 6. Finally, the conclusions are described in Section 7.

2. Flood Risk Assessment Methodology in Literature

Assessing the risk of flood is a widely conducted research and applied several strategies and modeling. This is categorized into three ways: (1) hydrological models such as SWAT [17], HEC-RAS [18,19,20,21] with digital terrain model (DTM) generation and vertical uncertainty assessment [22], ADCIRC model [23], FLO-2D [24]; (2) quantitative approaches such as frequency ratio [25], logistic regression [25], weights-of-evidence [25], genetic algorithm [26], differential evolution [26] and analytic hierarchy process [1,4,13] and (3) non-linear machine learning algorithms such as decision tree [27], artificial neural network [28], support vector machine [29] and random forest [30].
The abovementioned methods successfully provided flood susceptibility assessment, but every method faced some limitations. In hydrological modeling, the preparation and calibration of parameters are time-consuming [26], the quality of DEM as input will affect the performance of the model [21,23,27], and it needs high requirements of computer resources. The most common and challenging problem in hydrological modeling during the analysis is also the scarcity or lack of hydro-meteorological data. In particular, data issues are critical in ungauged catchments. The area proximity, ordinary least squares and regression analysis are among the commonly used methods to solve this problem. However, as stated by Haddad et al. [31], the use of regression analysis has limitations, which involves the trading space for time. Yet, Haddad et al. [31] compared the two regression-based methods that employ a Bayesian generalized least squares framework: (i.e., quantile regression technique and parameter regression technique) to understand more the contributing factors in flooding. In quantitative modeling like statistical and data-driven methods increased the subjectivity as it requires an expert to select the flood conditioning factors. In the non-linear machine learning algorithms, it may lead to poor projections due to the large and inconsistent value ranges in the datasets [26,32]. In this study, the quantitative approach, AHP, was used for flood risk assessment due to data scarcity in the study area. Furthermore, other quantitative approaches such as WR, and RW were also applied for flood risk assessment, and their results were compared with that from AHP. These approaches have the capability to include multiple datasets and can be used on limited computer resources. Using multiple data sources potentially also allows for more reliable disaster risk reduction plan, especially in the context of flood management and mitigation. For multi-data sources, the multi-criteria analysis approach has become widely used to solve complex phenomena and to assess flood risk [1,4,30,33].
AHP is a popular and widely used method developed by Saaty [34] to solve a wide range of multi-criteria decision-making problems, combined with the use of a pairwise comparison matrix that calculates the weights of every criterion considered. Recently, researchers have given considerable attention to using AHP in flood risk assessment [1,4,13,35,36]. In those case studies, it has been shown that AHP can assess and map flood risk areas with reliable accuracy. However, the result is based on expert opinions and thus may be subjected to intellectual limitations due to uncertainty and subjectivity.
In Danumah et al. [1], AHP was used to integrate several elements under two criteria, hazard and vulnerability, for flood risk assessment and mapping. The slope, drainage density, type of soil, and isohyet were the elements used to generate their hazard map. The land use, population density, and drainage system were the elements used to create their vulnerability map. Danumah et al. [1] evidenced the reliability and the powerful role of geoinformation techniques in natural disaster assessment that requires the contribution of multi-source data. In Elkhrachy [2], a flash flood map was generated by using satellite images and the AHP was used to determine the relative impact or weight of flood-causative factors: the runoff, soil type, surface slope, surface roughness, drainage density, distance to the main channel and land use. Gigović et al. [4] conducted a study of the hazard zone mapping of urban flood-prone areas using a GIS multi-criteria methodology by considering six factors relevant to urban flooding: the height, slope, distance to the sewage network, distance from the water surface, water table and land use. Their results indicated that the scenario in which the methodology was used provides the highest level of compatibility with historical flooding data. The case studies of Danumah et al. [1], Gigović et al. [4] and Elkhrachy [2] were, however, limited to show the resulting flood hazard map or flood-prone area without validations of their resulting maps with ground trothing points or historical flood events. In this study, the resulting flood-prone area by AHP is validated by comparing with ground trothing points of historical flood events from a field survey. Moreover, the sensitivity test of consistency ratio is conducted to investigate the effect of varying ratio on flood-prone area.

3. Study Region and Dataset

3.1. Case Study Area

In the Philippines, the occurrence of natural hazards and disasters is frequent due to the physical environment of the country, which faces the Pacific Ocean where catastrophic typhoons originate from. Furthermore, the Philippines is located in a part of the Pacific Ring of Fire. The Philippines is vulnerable to typhoons, floods, earthquakes, storm surges and tsunamis [37]. The primary causes of the disasters in this area are typhoons and flooding because of their frequent occurrence and their magnitude of impact on society. Approximately 20 typhoons per year approach or make landfall, and this is the highest frequency of typhoon events in the entire world [9].
On the other hand, Mindanao Island, as described, is known to be rather typhoon free regions with less flood and storm surge risk in the Philippines. However, disasters due to tropical storm Sendong (international name, Washi) in 2011 and Super Typhoon Pablo (international name, Bopha) in 2012 changed that; Mindanao is now considered to be in one of the typhoon paths. These two deadly storms killed more than 1000 persons and led to approximately 100 missing persons because of heavy rains that triggered the rivers to overflow and that caused flash floods and landslides [38,39]. Davao Oriental is located in the southeast of Mindanao Island, and disasters due to floods and storm surge in this region are anticipated to be more frequent under global warming. Therefore, Davao Oriental is considered for the study area.
Davao Oriental is a province of Davao Region, Philippines. Davao Oriental is the easternmost province of the country and is located between 6°20’ and 7°10’ N latitude and 125°0’ and 126°20’ E longitude (Figure 1). The province is composed of 183 barangays (towns) and two congressional districts. District 1, also known as the East Coast, consists of five municipalities; Tarragona, Manay, Baganga, Cateel and Boston. District 2, also called Davao Gulf, and includes four municipalities and one city: San Isidro, Governor Generoso, Banaybanay, Lupon and Mati City. Davao Oriental covers an area of approximately 5679 km2. The population is approximately 558,958 with a population density of 98 km2.
The topographic condition of Davao Oriental is characterized by a widespread chain of mountain ranges with an uneven distribution of plateaus, swamps and lowlands. The Mount Hamiguitan range is the newly enlisted UNESCO Heritage located at the administrative boundaries between the municipality of San Isidro, Governor Generoso and Mati City. This province occupies the largest land area of the provinces of Region XI (Davao Region): approximately 516,446 hectares or 26% of the total land area of Davao Region.
Davao Oriental is unique in that it is the only province in the country where all municipalities have a coastline. The coastline of the province measures 513.2 km from the municipality of Boston in the northern part of the province to the municipality of Banaybanay in the southwestern part of the province; it is the longest section of coastline in the country and approximately 3% of the total coastline of the country.

3.2. Dataset

3.2.1. Rainfall

The observed rainfall data are obtained from the Hinatuan and DOST-RXI Stations (Figure 2). The Hinatuan Station has observed daily rainfall data from 1973 to 2015. On the other hand, the DOST-RXI Station has observed daily rainfall data from 2006 to 2015. These two stations are far from the boundary of Davao Oriental. The distances of the Hinatuan and DOST-RXI Stations from the nearest border of the province are 43.9 km and 33.4 km, respectively, as shown in Figure 1D. Therefore, utilizing the rainfall data from these two stations will give an unreliable result because of their geographic locations.
To address the geographic limitation of the weather stations, rainfall data were obtained from the National Climatic Data Center (NCDC) and the Global Precipitation Climatology Centre (GPCC; see Table S1 in Supplementary Material 1). The NCDC rainfall data are global daily precipitations with a spatial grid resolution of 0.5° latitude × 0.5° longitude. In contrast, the GPCC rainfall data are global daily precipitations on a regular grid with a spatial resolution of 1.0° latitude × 1.0° longitude. The rainfall data were extracted within the boundary of Davao Oriental. Additionally, the nearest grid points were compared to the Hinatuan Station and DOST-RXI Station to evaluate the reliability of the data, as shown in Figure 1C,D. Figure 2 reveals that the NCDC rainfall data show good agreement with the observed rainfalls at the two stations. The GPCC rainfall, however, depicts lower values at Hinatuan while the rainfall at DOST-RXI shows good agreement with the observed and NCDC rainfalls. The rainfall data were interpolated as most of the related studies of flood hazard mapping used an interpolation method to create the rainfall distribution map [1,5,15,40].

3.2.2. Digital Elevation Model

The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) is one of the most widely used DEM datasets. ASTER GDEM has been applied in many fields, such as soil erosion, topography, geomorphology and hydrology [41,42]. ASTER GDEM is also widely used in developing flood hazard maps to extract drainage networks, flow accumulations and directions, basin boundaries, watershed boundaries, slopes and elevations [43,44,45,46].
The first version of ASTER GDEM, released in June 2009, was generated using stereo-pair images collected by the ASTER instrument onboard the Terra satellite. ASTER GDEM coverage spans from 83° north latitude to 83° south, encompassing 99% of the landmass of Earth [47]. The latest version of the ASTER GDEM V2 data was released on October 17, 2011, as shown in Figure 1D for the study area, Davao Oriental. The improved ASTER GDEM V2 dataset adds 260,000 additional stereo-pairs, improving coverage and reducing the occurrence of artifacts. The refined production algorithm provides improved spatial resolution, increased horizontal and vertical accuracy and superior water body coverage and detection [47]. Therefore, the ASTER GDEM V2 dataset is of better quality than the first version. ASTER GDEM V2 has a 30 m spatial resolution in the GeoTIFF image format with decimal degrees and WGS84 datum.

3.2.3. Administration Boundary

The Administrative boundaries of Davao Oriental include provincial, municipal and barangay boundaries. Davao Oriental is the easternmost province of the country. On the west side of Davao Oriental is the province of Compostela Valley, and the provinces of Surigao del Sur and Agusan del Sur are to the north. The Philippine Sea, part of the Pacific Ocean, is to the east of Davao Oriental. The administrative boundaries include ten municipal boundaries and 183 barangay boundaries. The administrative boundaries were provided as geographical information system (GIS) data from the global administrative areas and Philippine GIS organization, as indicated in Table S1. These GIS data are in decimal degrees and have a WGS84 datum. The data were then projected to the UTM coordinate system zone 51 N.

3.2.4. Population and Socio-Economic Data

Based on the 2015 census of population and housing [48] as shown in Table 1, the province of Davao Oriental has a total population of 558,958 in 2015. The 2015 census includes 41,340 more persons than were counted in the 2010 census, which determined a total population of 517,618 persons. This increase in the population from 2010 to 2015 translates into an average annual population growth rate of 1.47%.
Mati City has the highest population of all municipalities, with 25.3% of the total provincial population. The municipality of Lupon is the second largest, with 11.8% of the total provincial population, followed by the municipalities of Baganga and Governor Generoso, with 10.1% and 9.9%, respectively. The rest of the municipalities contribute 43% of the total provincial population. Visayan migrants inhabit the province. Ethnic groups include the Mandaya, Mansaka, Manobo and Kalagan. The native languages spoken in the province are Bisaya, Mandaya, Dabawenyo and Chavacano. Due to the geographic location of Davao Oriental, the primary source of income is farming and fishing. The coastline of the province plays a vital role in the fishing activities of its fisherfolk.

3.2.5. Soil Type

The Davao Oriental soil cover is mainly loam and sandy clay loam, and a section of rough mountainous land has an unidentified soil type. The area of Davao Oriental is classified into four sets of input parameters based on United States Department of Agriculture Natural Resources Conservation Service (USDA-NRCS) [49]. The Camasan sandy clay loam and undifferentiated mountain soil are classified as sandy clay loam. The San Manuel silty clay loam and San Miguel clay loam are classified as clay loam. The Malalag loam is classified as loam, and the Bolinao clay is classified as clay.

4. Methodology

In this paper, the methodology is based on a GIS-based spatial assessment process for flood hazard and conducted by using multi-criteria analysis concepts and an AHP, WR and RW frameworks. This approach used the spatial data management capabilities of GIS and the flexibility of MCDA to combine factual evidence with value-based information. MCDA is a systematic procedure for designing, evaluating and selecting decision alternatives on the basis of conflicting criteria. The motivation of integrating GIS and MCDA stems from the need to make the geographic information technology more relevant for analyzing decisions. Spatial or GIS-based MCDA can be thought of as a collection method for transforming and combining geographic data and preferences (value judgments) to obtain information for decision-making. The factual evidences considered in this paper were slope, elevation, soil type, rainfall, drainage density and distance to the main channel, and defined as indicators. The value-based information approach, following Saaty [34], uses an expert decision to identify which indicators are the most crucial in flood hazard map. Figure 3 presents the flowchart of the methodology applied in this paper.

4.1. Indicators and Criteria

An important step of this analysis is the criteria selection for evaluating the flood risk map. The criteria considered are two, flood hazard and vulnerability. The barangay population density was used as the indicator in the vulnerability assessment. There are many indicators affecting flood hazard identification and modeling, and they vary from one study area to another. For instance, urban flood modeling is incredibly complex compared to rural flood modeling due to the interactions with manmade structures such as buildings, roads, channels, tunnels and underground structures. This paper used a composite flood hazard map index based on following six indicators. These indicators were selected based on various case studies [1,2,15,16,33] and based on the available data in the study area.
(1) Rainfall: At any location, the chance of flooding increases as the amount of rain increases. Higher rainfall intensity can result in more runoff because the ground cannot quickly absorb the water. In this study, the 5-day continuous maximum rainfall in every year was determined for 1991–2016 as shown in Figure 2. Then, the average of the 5-day continuous maximum rainfall over the 25 years was used in the analysis. Due to the limited number of weather stations in the study area, data from the NCDC was used in this study (Table S1). A spatial interpolation was performed to obtain the data for the study area by using the Kriging interpolation method, as in other case studies [5,50,51,52]. The spatial distribution of the 5-day continuous rainfall is illustrated in Figure 4A.
(2) Slope: The slope is one of the crucial elements in flooding. The danger of flooding increases as the slope decrease. The slope is a reliable indicator of flood susceptibility [53]. The slope is presented in Figure 4B.
(3) Elevation: In contrast, the elevation of an area is a major factor in flooding. Low elevation is a good indicator of areas with a high potential for flood accumulation. Water flows from higher to lower elevations; therefore, slope influences the amount of surface runoff and infiltration [5]. Flat areas at low elevations may flood more quickly than areas at higher elevations with steeper slopes [5]. The elevation is appeared in Figure 4C.
(4) Soil: Soil type is a significant factor in determining the water holding and infiltration characteristics of an area and consequently affects flood susceptibility [54,55]. Generally, runoff from intense rainfall is likely to be faster and greater in clay soils than in sand [56]. Additionally, rain runoff from intense rainfall is likely to be faster and greater in loam than in sand. The soil type of Davao Oriental is shown in Figure 4D.
(5) Drainage density: Drainage density is the length of all channels within the basin divided by the area of the basin [2]. A dense drainage network is a good indicator of flow accumulation pathways and of areas with a high potential for flooding [53]. The drainage density and distance to the main channel are exhibited in Figure 4E.
(6) Distance to the main channel: Distance to the main channel significantly impacts flood mapping. Areas located close to the main channel and flow accumulation path are more likely to flood [4,53].
Once the six indicators were defined, the next step was to build the spatial database. Each indicator was converted into raster data with a 30 m × 30 m grid resolution. The ASTER GDEM data were registered and projected to the UTM coordinate system zone 51N. The slope and elevation were obtained using the 3D Analyst algorithm based on a DEM. The drainage density and distance to the main channel were obtained using Arc Hydro, which is a set of data models that operate within ArcGIS to support geospatial and temporal data analyses. All data were integrated into the GIS environment using the AHP, WR and RW methods. Finally, the weighted overlay was used to calculate the flood hazard, which was then combined with the flood vulnerability to assess a flood risk in the study area.

4.2. Analytic Hierarchy Process

The AHP is a systematic approach developed in the 1970s to give decision-making based on experience, intuition and heuristics derived from sound mathematical principles [57]. It helps structure the decision-maker’s thoughts and can help in organizing the problem in a way that is simple to follow and analyze. Basically, the AHP helps in structuring the complexity, measurement and synthesis of rankings [58]. Therefore, these features make it suitable for a wide variety of applications, including alternative selection, resource allocation, forecasting, business process re-engineering, quality function development, benchmarking, public policy decisions, healthcare, hazard mapping, risk analysis and many more [57].
The AHP provides a means of decomposing the problem into a hierarchy of sub-problems, which can be more easily comprehended and subjectively evaluated. The subjective evaluations are converted into numerical values and processed to rank each alternative (indicator in this study) on a numerical scale. It can be explained in the following steps:
Step 1: The problem is decomposed into a hierarchy of goal, criteria and indicators. Structuring the decision problem as a hierarchy is fundamental and the most important part in the process of the AHP. The hierarchy illustrates a relationship between elements of one level with those of the level immediately below. Figure 3 shows a tree structure representing an inverted hierarchy structure. At the root of the inverted hierarchy is the goal (to produce a flood risk map) of the problem being studied and analyzed. The leaf nodes are the six indicators to be compared.
Step 2: Data are collected from experts or decision-makers corresponding to the hierarchy structure, in the pairwise comparison of indicators based on a qualitative scale between 1 (equally important) and 9 (extremely important) proposed by Saaty [34] (see, Table S2 for the Saaty scale in the Supplementary Material). For each pairing within each criterion, a higher number means that the chose indicator is considered more important than the other indicator used in the comparison. Four local experts from the province of Davao Oriental in the field of natural disaster provided their judgments of the relative importance of one indicator against another (Table S3).
Step 3: The pairwise comparison of the six indicators generated at step 2 above are organized into a square matrix, the pairwise comparison matrix (Table 2). The diagonal elements of the matrix are 1. The indicator in the ith row is more important that the indicator in jth column if the value of element (i,j) is larger than 1; otherwise the indicator in the jth column is more important than in the ith row. The (j,i) element of the matrix is the reciprocal of the (i,j) element.
Step 4: The pairwise comparison matrix from step 3 is then normalized by Equation (1). Then, the priority vector (PV) is derived by dividing the sum of the normalized column of the matrix by the number of indicators, n, as shown in Equation (2).
X ij = C ij / i = 1 n C ij ,
PV ij = j = 1 n X ij   / n ,
where Cij is the value of an indicator in the pairwise comparison matrix, Xij is the normalized score and PVij is the priority vector. The PVs give the relative importance of the indicators being compared and represent the weights of the indicators. Table 3 depicts the normalized comparison matrix and the PVs.
Step 5: The consistency of the normalized comparison matrix is evaluated as described below. Since the numeric values of the matrix are derived from the subjective judgments of experts, it is impossible to avoid some inconsistency in the final matrix of their judgments [59]. Then, the question is how much inconsistency is acceptable. For this purpose, the AHP calculates a consistency ratio (CR) comparing the consistency index (CI) of the final normalized comparison matrix with that of a random-like matrix (RI). Saaty [34] provides the calculated RI value for matrices of different sizes (Table S4).
Firstly, use the PV as weights of each column of the pairwise comparison matrix as shown in Table S5. Secondly, multiply each value in the first column of the pairwise comparison matrix in Table S5 by the PV of first indicator as shown in the first column of Table S6; continue this process for all the columns. Table S6 shows the resulting matrix. Thirdly, add the values in each row to obtain the weighted sum as shown in Table S7. Fourthly, divide the weighted sums by the corresponding PVs of each indicator to obtain the consistency measure (CM) as shown in Table S8. Fifthly, calculate the average of CMs, λ max . Lastly, calculate the CI and CR by the followings.
CI =   ( λ max     n )   /   ( n     1 ) ,
CR = CI / RI ,
The λmax was 6.121 obtained from Table S8 and then the CI was 0.024 with n = 6. Therefore, by dividing the CI by RI (1.24 for n = 6 from Saaty), the CR was obtained as 0.02. Saaty suggests the value of CR should be less than 1 to have an acceptable consistency. Since the CR was 0.02 in this case, it can be assumed that the comparison matrix was reasonably consistent so that we could continue the decision-making process using AHP.
Step 6: By multiplying the six indicators by the corresponding weights and aggregating them, the hazard index (HI) is obtained by:
HI =   St   ×   0.041   +   Sl   ×   0.238   +   Dd   ×   0.061   +   Dc   ×   0.095   +   E   ×   0.147   +   R   ×   0.418 ,
where St, Sl, Dd, Dc, E and R represent the six indicators, soil type, slope, drainage, distance to main channel, elevation and rainfall, respectively.

4.3. Weights by Rank

In this method, the weight of indicator is arranging through a rank order. In that way, every indicator is ranked as per decision-maker preference. The decision-makers were the same four local experts considered in the AHP analysis. The most important is 1, the next is 2 and so on. The ranking of the indicators from the experts are shown in Table S9. Once ranking is established for a set of indicators, the numerical weights are calculated using Equation (6).
W i = ( n r i + 1 )   /   j = 1 n ( n r j + 1 ) ,
where Wi is the normalized weight for ith indicator, n is the total number of indicators under consideration (j = 1, 2, …, n) and rj is the rank position of jth indicator. Each indicator is weighted (n − ri + 1) and then normalized by the sum of all weights ∑(n − rj + 1). Therefore, the ranking method estimated weight should be considered as an approximation [60]. The results are given in Table 4. The average rank (AR) values are the average results from the four experts in the study area. Then, the resulting weight (W) in Table 4 was used to calculate the HI as shown in Equation (7).
HI =   St   ×   0.07   +   Sl   ×   0.25   +   Dd   ×   0.11   +   Dc   ×   0.11   +   E   ×   0.20   +   R   ×   0.26.

4.4. Ratio Weighting

In this method, it also uses the pairwise comparisons to establish the relative importance among indicators. It is also likely that the pairwise comparisons will be inconsistent. Therefore, this approach requires more information to compensate for the inconsistencies of judgments. Efficient techniques to retrieve weights from pairwise comparison data include the simplified prioritization method as follows [61]:
Step 1: A decision-maker assesses importance (weight) ratios between indicators using pairwise comparisons. In this study, the same local experts in Davao Oriental ranked the indicators. Table S10 shows the comparison matrix based on the expert’s judgments among the indicators.
Step 2: Compute the geometric mean of each row of the comparison matrix in Table S10, and then normalize the resulting numbers. The geometric mean and the numerical weights are calculated using Equations (8) and (9). Table 4 depicts the resulting geometric mean and weights of the judgments.
GM i = ( i = 1 n x i   ) 1 / n ,
W i = ( GM i )   /   i = 1 n ( GM i ) ,
where GMi is the geometric mean for ith indicator, n is the total number of indicators and Wi is the normalized weight for ith indicator. The HI is computed by using the average weight (AW) in Table 4 as in Equation (10).
HI =   St   ×   0.04   +   Sl   ×   0.31   +   Dd   ×   0.06   +   Dc   ×   0.09   +   E   ×   0.17   +   R   ×   0.33.
Finally, the HI obtained by the three models, AHP, WR and RW, was further processed using a weighted overlay analysis in ArcGIS. The results of weighted overlay analysis were classified using equal interval into four levels such as very low, low, moderate and high.

4.5. Validation

Super Typhoon Bopha was considered as one of the world’s damaging storms in 2012. Bopha entered the Philippine Area of Responsibility (PAR) as a rapidly intensified Category 5 Super Typhoon and made its landfall at Baganga, Davao Oriental, of 4 December 2012. It has reached an average speed of 185 kph and gusts reaching 210 kph. Typhoon Pablo, as what Typhoon Bopha is called in the Philippines, was the most powerful storm to have hit the island of Mindanao, southern Philippines, in more than 100 years of recorded storms [60]. Its torrential rains generated massive debris flow in the Mayo River watershed in the Andap village in New Bataan municipality, causing some areas to be buried under a rubble as thick as 9 m and killing 566 people [60]. Moreover, it was also revealed that a storm surge with a height of approximately 1 m has reached a distance of 500 m inland away from the coastline [12]. Unfortunately, the government has no accounts in recording the flooded areas during the said event.
Due to the hazards caused by storm surge, storm surge prediction and forecasting currently is actually progressing. In fact, in a study conducted by Ross, et al. [60] regarding the projected susceptibility to storm surge modeling using historical records from Philippine Atmospheric, Geophysical and Astronomical Services Administration (PAGASA), it showed that there was a storm surge event as high as 3.66 m. The said storm surge happened in the city of Bislig, located about 60 km from the northern coastal boundary of the Davao Oriental province. Hence, to address such scenarios, this study was conducted using a GPS-based field survey to determine the flooded areas during Typhoon Bopha.
In order to validate the resulting flood-prone areas by the AHP, WR and RW methods, a GPS-based field survey to local people was carried out along the east coast of Davao Oriental to investigate ground true flooded points using historical flooding events. This survey was carried out to determine the flooded area during the Typhoon Bopha. It is to note that the field survey was also not carried out in the municipalities located near the Davao Gulf thus, the points do not cover all barangays specified in the study area. Hence, only 70 ground truthing points were collected and were used for performance evaluation (see Figure 1D for the locations and Supplementary Material 3 for the survey result).
The performance evaluation was done using the accuracy assessment of the flood classification. One of the commonly used methods is to apply a confusion matrix or error matrix. This method can be used to compute several assessment elements such as overall accuracy (ACC), true positive rate (TPR), true negative rate (TNR), false positive rate (FPR) and false negative rate (FNR) using the Equations (11)–(15), respectively.
ACC = TP + TN P + N = TP + TN TP + TN + FP + FN ,
TPR = TP TP + FN ,
TNR = TN TN + FP ,
FPR = FP FP + FN = 1 TNR ,
FNR = FN FN + TP = 1 TPR ,
where P, N, TP, TN, FP and FN are condition positive, condition negative, true positive, true negative, false positive and false negative, respectively.
To evaluate the performance, an accuracy assessment was performed using the confusion matrix, and it was calculated based on the 70 ground points collected from the field survey. Elevations were extracted from these truthing points, and it showed that it is in the range between 10 m and 20 m above mean sea level. After spatially interpolating the points, it was identified that all areas in Davao Oriental that belongs to the below 20 m elevation were considered as susceptible flood-prone areas.
Based on the survey, a resulting map was created (as seen in Figure 1C) and was overlaid with the prediction maps obtained by the AHP, RW and WR methods to compare and assess the accuracy of the prediction maps. In the accuracy assessment, the TP and FP are the numbers of pixels correctly classified and the numbers of pixels mistakenly classified, respectively. With this, the results showed that the number of correctly predicted pixels by the model (TP) are plotted to the number of incorrect predicted pixels (FP), as illustrated in Table 5 and Figure 5 , indicating that the AHP-based flood area is truly flooded in the ground truthing-based result.
Data scarcity and sparsity are the main problems solved in this approach. In this paper, data scarcity means limited historical flood events, or the recorded floods are limited as compared to the amount needed in the validation of the model (see Figure 1D). Moreover, as shown in Figure 1D, the survey data of flood are not evenly distributed (i.e., sparsity) to the study area. In this study, the flood data used was derived from the flooding caused by typhoon Bopha and did not cover all flood events happened in Davao Oriental. This approach has uncertainty and limitation, same with other hydrological modeling. Fortunately, the GIS-based flood validation approach is a good start for the researcher to conduct flood susceptibility mapping in data-scarce and sparse region.

5. Results

5.1. Vulnerability Map

Vulnerability expresses the level of inability to resist a hazard or to respond when disaster has occurred [62]. For example, people who live in low-lying areas are more vulnerable to floods than people who live at higher elevations. Moreover, vulnerability is the most crucial component of flood risk because vulnerability determines if exposure to a hazard constitutes a threat. Flood vulnerability mapping is the process of determining the degree of susceptibility and exposure of a given place to flooding. The susceptibility and exposure issues include several factors, such as the age and health of the population, the socio-economic activities, the quality of buildings and their location with respect to flooding. In this study, due to the limited socio-economic data from the provincial office of the Philippine Statistics Authority in Davao Oriental, the only indicator used for the vulnerability map is the population density from barangay (town). Thus, the vulnerability map obtained from the barangay population density. Vulnerability map is divided into five classifications (shown in Figure 4F). Barangays with high populations are mostly proximal to coastal areas and riverside.

5.2. Hazard Map

Hazard is the function of scale and probability of natural disaster to cause harm or loss in a specified location in a period of time. Thus, it reflects the natural properties of the disaster [36], and the physical and statistical aspects of the flooding process [63]. Flood hazard maps contain information about the probability and/or magnitude of an event [63], the return period of the flood, the extent, and depth of inundation [64]. Due to the exclusion of the hydrological modeling of this study, the depth of inundation, and the return period of the flood were not included. Thus, in this study, the flood hazard map only covered the extent or the possibility of the spatial distribution of flood in the study area as shown in Figure 4G. Figure 4G shows that the hazard map was a combination of maps: rainfall (Figure 4A), slope and elevation (Figure 4B,C), soil (Figure 4D) and drainage density and distance to the main channel (Figure 4E). Each map had different weights that, in this paper, the rainfall was emphasized more. The rainfall map (Figure 4A) shows that each east coast municipality had a higher incidence of rainfall compared to that of the Davao Gulf municipalities. The slope and elevation maps (Figure 4B,C) indicate that most of the areas in Davao Oriental were at low elevations, less than 217 m and have slopes up to 16 degrees. All maps were classified based on equal interval classification.
The resulting flood hazard maps by the AHP, WR and RW methods are illustrated in Figure 6. Table 5 shows the confusion matrix of the three models. The data in Table 5 was the number of pixels from the models (i.e., predicted) in high classification, and the observed data from the flood and non-flood area that are shown in Figure 1C. In the overall assessment, AHP had superior results in ACC, TNR and TPR as shown in Figure 5A. Moreover, Figure 5B shows the accuracy of the positive prediction (TPR) compared to the wrong prediction (FPR). The upper part of the red line in the Figure 5B means the model was accurate, otherwise not accurate. Therefore, AHP shows higher accuracy compared to the other two methods. In overall, the analysis reveals that AHP was much more accurate and reliable for flood risk analysis in this study.
Finally, in AHP, the very low and low classes covered approximately 1% and 22% of the total area of Davao Oriental, respectively. These are areas with high elevations, high slopes and low drainage densities. The categories of moderate, and high classes were estimated to cover 55%, and 22% of the total area of Davao Oriental, respectively, and are mostly in the east coast municipalities (Boston, Cateel and Baganga). For the high hazard risk areas of 22.15% of Davao Oriental, the rainfall (41.8%) and slope (23.8%) were the most significant causative factors of flood occurrence. In the flood-prone area, 64 out of 183 barangays (towns), and 30% (168,034/558,958) of the population were at a high of flooding at the current situation as shown in Table 6. Additionally, these were areas that had high occurrences of rainfall, high drainage densities, low elevations and low slopes, and they were proximal to river channels and shorelines. In the Davao Gulf municipalities, most of the areas of Mati City and Lupon belonged to the moderate hazard class due to their geophysical characteristics. Lupon had one of the widest rivers of all the municipalities in the province. On the other hand, Mati City was a low-lying city, especially in the city centre where the urbanized area was located.

5.3. Risk Map

In the study of Zhou et al. [36], the risk is the product of hazard and vulnerability as defined by the United Nations, Department of Humanitarian Affairs [65,66] and World Meteorological Organization [67]. A study by Thieken et al. [63] also defined flood risk as to the product of hazard and vulnerability.
The flood risk assessment has been widely used in understanding flood disaster evaluation, warning and awareness [36]. It is essential in management and decision making. In this study, flood risk map is a combination of hazard map and vulnerability map. The result of the flood risk map was generated by a weighted overlay of the hazard and vulnerability maps with equal weights. The flood risk map was classified into five classes (Figure 4H). The very low, low and moderate classes covered 0.63%, 60.64% and 35.35% of the total area, respectively. The categories of high and very high were estimated at 3.34% and 0.05% of the total area, respectively. These areas were barangays with high population densities and were mostly in low-lying areas near the shoreline and rivers. On the other hand, the moderate and very low classes were areas at high elevations with smaller populations. Thirty-one out of one hundred eighty-three (31/183) barangays were at high risk of flooding. Only 1 barangay had a very high risk of flooding at the current climate, calling for the immediate attention of the government to develop strategies or flood protection for the future occurrence of flooding within the province of Davao Oriental.

6. Effects of Varying Consistency Ration (CR) on Flood Hazard Map

Integrating multiple indicators (data sources) in the multi-criteria analysis presents a real advantage, and results show that the AHP approach allowed a better understanding of all the indicator contributions in the flood assessment because a weight was given to each indicator. However, data from different sources with different resolutions were factors of bias during processing and analysis. On the other hand, the addition of weights reduced the bias and uncertainty in the result.
However, the AHP method shows some limitations due to the subjectivity in choosing the value of the indicator weighting from the arbitrary judgments of experts [1,6,13]. The CR is the crucial element when generating HI for the multi-criteria analysis using AHP method. According to Saaty [34], this subjectivity problem was reduced by the CR test. The CR threshold must be less than 0.10 or 10% to make a coherent judgment. The lower the CR, the better, and if CR is 0, it means the HI is perfect. Unfortunately, the range of CR between 0 and 9.9% is a rough estimation. The HI is coherent if the CR is 9.9% or CR is 0%, but the weights of the indicators vary in every CR generation.
Therefore, the sensitivity of HI was tested in different three CR scenarios; (1) the lowest CR was 2.0%, (2) the highest CR was 8.6% and (3) the middle CR was 5.8%. Figure 7 shows the flood hazard maps from the sensitivity analysis in different CR scenarios. Table 7 shows the weights of every indicator according to the ratio of indicators. Table 8 shows the rank of every classification in different CR.
The result illustrates that the ranking of classifications is the same in the three CR scenarios. Average change according to classification was also ranging from −0.77% to 1.21%, and the overall average change was only 0.02%. Therefore, it was found that the changes of the CR did not significantly affect the result of the hazard map, and the changes were acceptable. Thus, the lowest criteria ratio was more appropriate to use as according to Saaty [34] so that the lower CR made the judgment more coherent. Therefore, the CR with 0.02 was used for to obtain the flood risk map.

7. Conclusions

Flooding is a natural hazard that poses a risk to both population and physical structures within the affected areas. Flood risk occurrence can be better understood with the use of a spatial assessment. Thus, in this study, the main objective was to identify flood-prone and categorize risk areas to produce a flood risk assessment based on several factors. For this purpose, a GIS-based flood risk assessment was conducted by using the concepts of AHP, WR and RW frameworks in the province of Davao Oriental, and determines which model has a superior result.
The analysis was based on the integration of multiple indicators including rainfall, slope, elevation, drainage density, soil type and distance to the main channel. The AHP model depicts a superior result in the computation of relative importance of indicators compared to WR and RW methods. The result of the AHP model shows that Davao Oriental is under high flood-prone areas approximately 22.15% of the province, and these are in low-lying areas near the coastline and rivers. The multi-criteria analysis allowed the integration of the several indicators under two criteria for flood risk assessment: hazard and vulnerability. The results showed that the province of Davao Oriental is moderately exposed to the risk of flooding, which requires immediate attention from the decision-makers to develop strategies for the future occurrence of flooding within the province of Davao Oriental. An effective strategy should be drawn up at the barangay level since each barangay has its own physical and social characteristics.
Immediate strategies like installing weather measurement instruments in at least one station per municipality will give us precise weather forecasting especially during extreme events like drought, heavy rainfalls, storms and typhoon. Since the national government has no specific flood-risk mitigation plan, the municipal government in coordination with the office of the disaster risk reduction and management council can design a flood-risk education program in every barangay. The flood-risk education program is an education campaign on how to act before, during and after the event of flooding.
Additionally, the results may be improved by adding land use/land cover through image analysis methods that use high-resolution images such as Landsat 8+. Hydrologic modeling in 2D/3D for efficient analysis in determining the depth of flood [19], the extent of flood and the damage estimate using damage curve of the flood is also recommended as shown in literatures [20,21,23,24,27].
Lastly for future works, in the selection of methods, the subjectivity of the quantitative approaches such as AHP, WR and RW can be improved using data mining and machine learning techniques.

Supplementary Materials

The following are available online at https://www.mdpi.com/2073-4441/11/11/2203/s1, The Tables S1–S10 can be found in the Supplementary Material 1. The Supplementary Material 2 is the Excel file demonstrating the AHP procedure to calculate the weighting factors used in this study. The Supplementary Material 3 includes the field survey result.

Author Contributions

J.S.C. developed the concept of this study under the supervision of H.S.L. The analysis was carried out by both authors. Both authors drafted the first version of the manuscript and worked on improving and finalizing the manuscript.

Funding

This research is partly supported the Grant-in-Aid for Scientific Research (17K06577) from JSPS, Japan.

Acknowledgments

The first author is supported by The Human Resource Development Scholarship (JDS), Japan.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Danumah, J.H.; Odai, S.N.; Saley, B.M.; Szarzynski, J.; Thiel, M.; Kwaku, A.; Kouame, F.K.; Akpa, L.Y. Flood risk assessment and mapping in Abidjan district using multi-criteria analysis (AHP) model and geoinformation techniques, (cote d’ivoire). Geoenviron. Disasters 2016, 3, 10. [Google Scholar] [CrossRef]
  2. Elkhrachy, I. Flash Flood Hazard Mapping Using Satellite Images and GIS Tools: A case study of Najran City, Kingdom of Saudi Arabia (KSA). Egyp. J. Remote Sens. Space Sci. 2015, 18, 261–278. [Google Scholar] [CrossRef] [Green Version]
  3. Forkuo, E.K. Flood hazard mapping using aster image data with GIS. Int. J. Geomat. Geosci. 2010, 1, 932–950. [Google Scholar]
  4. Gigović, L.; Pamučar, D.; Bajić, Z.; Drobnjak, S. Application of GIS-Interval Rough AHP Methodology for Flood Hazard Mapping in Urban Areas. Water 2017, 9, 360. [Google Scholar] [CrossRef]
  5. Kazakis, N.; Kougias, I.; Patsialis, T. Assessment of flood hazard areas at a regional scale using an index-based approach and Analytical Hierarchy Process: Application in Rhodope–Evros region, Greece. Sci. Total Environ. 2015, 538, 555–563. [Google Scholar] [CrossRef] [PubMed]
  6. Papaioannou, G.; Vasiliades, L.; Loukas, A. Multi-Criteria Analysis Framework for Potential Flood Prone Areas Mapping. Water Resour. Manag. 2015, 29, 399–418. [Google Scholar] [CrossRef]
  7. Pedersen, A.; Mikkelsen, P.; Arnbjerg-Nielsen, K. Climate change-induced impacts on urban flood risk influenced by concurrent hazards, J. Flood Risk Manag. 2002, 5, 203–214. [Google Scholar] [CrossRef]
  8. Vino, P. Flood hazard assessment of Vamanapuram River Basin, Kerala, India: An approach using Remote Sensing & GIS techniques. Adv. Appl. Sci. Res. 2013, 4, 263–274. [Google Scholar]
  9. Otieno, J.A. Scenario study for Flood Hazard Assessment in the lower Bicol Floodplain Philippine using A 2D Flood Model; International Institute for Geo-information Science and Earth Observation: Enschede, The Netherlands, 2004. [Google Scholar]
  10. Republic of the Philippines. Republic Act 101211: An act strengthening the Philippine disaster risk reduction and management system, providing for the national disaster risk reduction and management framework and institutionalizing the national disaster risk reduction and management plan, appropriating funds therefor and for other purposes. Available online: http://www.ndrrmc.gov.ph/attachments/article/45/Republic_Act_10121.pdf (accessed on 8 July 2018).
  11. Lagmay, A.M.F.A.; Racoma, B.A.; Aracan, K.A.; Alconis-Ayco, J.; Saddi, I.L. Disseminating near-real-time hazards information and flood maps in the Philippines through Web-GIS. J. Environ. Sci. 2017, 59, 13–23. [Google Scholar] [CrossRef]
  12. Ross, M.; Santiago, J.; Lagmay, A.M. Integrating and Applying Technology in Response to the Super Typhoon Bopha Disaster. Procedia Eng. 2015, 107, 100–109. [Google Scholar] [CrossRef] [Green Version]
  13. Cabrera, J.; Lee, H. Impacts of Climate Change on Flood-Prone Areas in Davao Oriental, Philippines. Water 2018, 10, 893. [Google Scholar] [CrossRef]
  14. Chakraborty, A.; Joshi, P.K. Mapping disaster vulnerability in India using analytical hierarchy process. Geomat. Nat. Hazards Risk 2016, 7, 308–325. [Google Scholar] [CrossRef]
  15. Wang, Y.; Li, Z.; Tang, Z.; Zeng, G. A GIS-Based Spatial Multi-Criteria Approach for Flood Risk Assessment in the Dongting Lake Region, Hunan, Central China. Water Resour Manag. 2011, 25, 3465–3484. [Google Scholar] [CrossRef]
  16. Di Baldassarre, G.; Castellarin, A.; Montanari, A.; Brath, A. Probability-weighted hazard maps for comparing different flood risk management strategies: A case study. Nat. Hazards 2009, 50, 479–496. [Google Scholar] [CrossRef]
  17. Yu, D.; Xie, P.; Dong, X.; Hu, X.; Liu, J.; Li, Y.; Peng, T.; Ma, H.; Wang, K.; Xu, S. Improvement of the SWAT model for event-based flood simulation on a sub-daily timescale. Hydrol. Earth Syst. Sci. 2018, 22, 5001–5019. [Google Scholar] [CrossRef] [Green Version]
  18. Veleda, S.; Martínez-Graña, A.; Santos-Francés, F.; Sánchez-SanRoman, J.; Criado, M. Analysis of the Hazard, Vulnerability, and Exposure to the Risk of Flooding (Alba de Yeltes, Salamanca, Spain). Appl. Sci. 2017, 7, 157. [Google Scholar] [CrossRef]
  19. Arrighi, C.; Brugioni, M.; Castelli, F.; Franceschini, S.; Mazzanti, B. Flood risk assessment in art cities the exemplary case of Florence (Italy). In Proceedings of the Italian National Conference on Hydraulics, Hydrology and Water Works, Bari, Italy, 7–10 September 2014. [Google Scholar]
  20. Macchione, F.; Costabile, P.; Costanzo, C.; De Lorenzo, G.; Razdar, B. Dam breach modelling: Influence on downstream water levels and a proposal of a physically based module for flood propagation software. J. Hydroinform. 2016, 18, 615–633. [Google Scholar] [CrossRef]
  21. Ali, A.; Solomatine, D.P.; Di Baldassarre, G. Assessing the impact of different sources of topographic data on 1-D hydraulic modelling of floods. Hydrol. Earth Syst. Sci. 2015, 19, 631–643. [Google Scholar] [CrossRef] [Green Version]
  22. Zazo, S.; Rodríguez-Gonzálvez, P.; Molina, J.L.; González-Aguilera, D.; Agudelo-Ruiz, C.A.; Hernández-López, D. Flood hazard assessment supported by reduced cost aerial precision photogrammetry. Remote Sens. 2018, 10, 10. [Google Scholar] [CrossRef]
  23. Falcão, A.P.; Mazzolari, A.; Gonçalves, A.; Araújo, M.; Trigo-Teixeira, A. Influence of elevation modelling on hydrodynamic simulations of a tidally-dominated estuary. J. Hydrol. 2013, 497, 152–164. [Google Scholar] [CrossRef]
  24. Huang, Y.; Qin, X. Uncertainty analysis for flood inundation modelling with a random floodplain roughness field. Environ. Syst. Res. 2014, 3, 9. [Google Scholar] [CrossRef] [Green Version]
  25. Tehrany, S.M.; Shabani, F.; Neamah Jebur, M.; Hong, H.; Chen, W.; Xie, X. GIS-based spatial prediction of flood prone areas using standalone frequency ratio, logistic regression, weight of evidence and their ensemble techniques. Geomat. Nat. Hazards Risk 2017, 8, 1538–1561. [Google Scholar] [CrossRef]
  26. Hong, H.; Panahi, M.; Shirzadi, A.; Ma, T.; Liu, J.; Zhu, A.X.; Chen, W.; Kougias, I.; Kazakis, N. Flood susceptibility assessment in Hengfeng area coupling adaptive neuro-fuzzy inference system with genetic algorithm and differential evolution. Sci. Total Environ. 2018, 621, 1124–1141. [Google Scholar] [CrossRef] [PubMed]
  27. Tehrany, M.; Pradhan, B.; Jebur, M. Spatial prediction of flood susceptible areas using rule based decision tree (DT) and a novel ensemble bivariate and multivariate statistical models in GIS. J. Hydrol. 2013, 504, 69–79. [Google Scholar] [CrossRef]
  28. Kia, B.M.; Pirasteh, S.; Pradhan, B.; Mahmud, R.A.; Sulaiman, W.N.A.; Moradi, A. An artificial neural network model for flood simulation using GIS: Johor River Basin, Malaysia. Environ. Earth Sci. 2012, 67, 251–264. [Google Scholar] [CrossRef]
  29. Hong, H.; Tsangaratos, P.; Ilia, I.; Liu, J.; Zhu, A.; Chen, W. Application of fuzzy weight of evidence and data mining techniques in construction of flood susceptibility map of Poyang County, China. Sci. Total Environ. 2018, 625, 575–588. [Google Scholar] [CrossRef]
  30. Wang, Z.; Lai, C.; Chen, X.; Yang, B.; Zhao, S.; Bai, X. Flood hazard risk assessment model based on random forest. J. Hydrol. 2015, 527, 1130–1141. [Google Scholar] [CrossRef]
  31. Haddad, K.; Rahman, A.; Stedinger, J.R. Regional flood frequency analysis using Bayesian generalized least squares: A comparison between quantile and parameter regression techniques. Hydrol. Process. 2012, 26, 1008–1021. [Google Scholar] [CrossRef]
  32. Bui, T.D.; Pradhan, B.; Nampak, H.; Bui, Q.-T.; Tran, Q.-A.; Nguyen, Q.-P. Hybrid artificial intelligence approach based on neural fuzzy inference model and metaheuristic optimization for flood susceptibilitgy modeling in a high-frequency tropical cyclone area using GIS. J. Hydrol. 2016, 540, 317–330. [Google Scholar]
  33. Yahaya, S.; Ahmad, N.; Abdalla, R.F. Multicriteria Analysis for Flood Vulnerable Areas in Hadejia-Jama’are River Basin, Nigeria. Eur. J. Sci. Res. 2010, 42, 71–83. [Google Scholar]
  34. Saaty, T.L. The Analytic Hierarchy Process: Planning, Priority Setting, Resource Allocation; McGraw-Hill: NewYork, NY, USA, 1980. [Google Scholar]
  35. Ouma, Y.; Tateishi, R. Urban Flood Vulnerability and Risk Mapping Using Integrated Multi-Parametric AHP and GIS: Methodological Overview and Case Study Assessment. Water 2014, 6, 1515–1545. [Google Scholar] [CrossRef]
  36. Zou, Q.; Jianzhong Zhou, B.; Chao Zhou, B.; Lixiang Song, B.; Jun Guo, B. Comprehensive flood risk assessment based on set pair analysis-variable fuzzy sets model and fuzzy AHP. Stoch. Environ. Res. Risk Assess. 2013, 27, 525–546. [Google Scholar] [CrossRef]
  37. Gurenko, E.; Lester, R. Rapid onset natural disasters: The role of financing in effective risk management. World Bank Policy Res. Pap. 2004, 3278, 34. [Google Scholar]
  38. NDRRMC. Sitrep No. 36 re Effects of Typhoon "Pablo" (Bopha); National Disaster Risk Reduction and Management Council: Philippines, 2012; p. 47.
  39. NDRRMC. Final report on the effects and emergency management re Tropical Storm Washi; National Disaster Risk Reduction and Management Council: Philippines, 2012.
  40. Zhang, K.; Manning, T.; Wu, S.; Rohm, W.; Silcock, D.; Choy, S. Capturing the Signature of Severe Weather Events in Australia Using GPS Measurements. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2015, 8, 1839–1847. [Google Scholar] [CrossRef]
  41. Chatterjee, S.; Krishna, A.P.; Sharma, A.P. Geospatial assessment of soil erosion vulnerability at watershed level in some sections of the Upper Subarnarekha river basin, Jharkhand, India. Environ. Earth Sci. 2014, 71, 357–374. [Google Scholar] [CrossRef]
  42. Frey, H.; Paul, F. On the suitability of the SRTM DEM and ASTER GDEM for the compilation of: Topographic parameters in glacier inventories. Int. J. Appl. Earth Obs. Geoinf. 2012, 18, 480–490. [Google Scholar] [CrossRef]
  43. Huong, D.T.V.; Nagasawa, R. Potential flood hazard assessment by integration of ALOS PALSAR and ASTER GDEM: A case study for the Hoa Chau commune, Hoa Vang district, in central Vietnam. J. Appl. Rem. Sens. 2014, 8, 083626. [Google Scholar] [CrossRef]
  44. Othman, N.; Jafri, M.Z.M.; Lim, H.S.; Tan, K.C. Using ASTER GDEM and SRTM digital elevation models to generate contour lines over rugged terrain of Makkah. In Proceedings of the 2011 IEEE International Conference on Space Science and Communication (IconSpace), Penang, Malaysia, 12–13 July 2011; pp. 11–13. [Google Scholar]
  45. Ramani, S.; Balasubramaniam, R. Micro level erosion suseptibility zonation using ASTER GDEM : A case study of hill sub—Watershed in Kodaikanal, South India. In Proceedings of the International Conference on Civil Engineering and Infrastructural Issues in Emerging Economies, Thanjavur, India, 27–28 February 2013. [Google Scholar]
  46. Reddy, G.P.O.; Kumar, N.; Sahu, N.; Singh, S.K. Evaluation of automatic drainage extraction thresholds using ASTER GDEM and Cartosat-1 DEM: A case study from basaltic terrain of Central India. Egypt. J. Remote Sens. Space Sci. 2018, 21, 95–104. [Google Scholar] [CrossRef]
  47. Tachikawa, T.; Gesch, D.B. ASTER Global Digital Elevation Model Version 2—Summary of Validation Results. Jpn. Sp. Syst. 2011. Available online: https://ssl.jspacesystems.or.jp/ersdac/GDEM/ver2Validation/Summary_GDEM2_validation_report_final.pdf (accessed on 8 July 2018).
  48. Philippine Statistics Authority. Available online: https://www.psa.gov.ph/sites/default/files/attachments/ hsd/pressrelease/R11.xlsx (accessed on 6 June 2017).
  49. USDA-NRCS. Urban Hydrology for Small Watersheds. Available online: https://www.nrcs.usda.gov/Internet/FSE_DOCUMENTS/stelprdb1044171.pdf (accessed on 30 June 2018).
  50. Basconcillo, J.; Lucero, A.; Solis, A.; Sandoval, J.R.; Bautista, E.; Koizumi, T.; Kanamaru, H. Statistically Downscaled Projected Changes in Seasonal Mean Temperature and Rainfall in Cagayan Valley, Philippines. J. Meteorol. Soc. Jpn. Ser. Ii 2016, 94A, 151–164. [Google Scholar] [CrossRef] [Green Version]
  51. Shahabi, H.; Salari, M.; Bin Ahmad, B.; Mohammandi, A. Soil erosion hazard mapping in Central Zab Basin using Epm model in GIS environment. Int. J. Geogr. Geol. 2016, 5, 224–235. [Google Scholar] [CrossRef]
  52. Zhao, S.; Wang, L.; Cheng, W.; Liu, H.; He, W. Rectification methods comparison for the ASTER GDEM V2 data using the ICESat/GLA14 data in the Lvliang mountains. China. Environ. Earth Sci. 2015, 74, 6571–6590. [Google Scholar] [CrossRef]
  53. Islam, M.; Sado, K. Flood hazard map and land development priority map developed using NOAA AVHRR and GIS data. Asian J. Geoinform. 2000, 45, 3. [Google Scholar]
  54. Ebaid, H.; Farag, H.; El Falaky, A. Using GIS and remote sensing approaches to delineate potential areas for runoff management applications in Egypt. Int. J. Environ. Sci. Eng. 2016, 7, 85–93. [Google Scholar]
  55. Nyarko, B. Application of a rational model in GIS for flood risk Aassessment in Accra, Ghana. J. Spat. Hydrol. 2002, 2, 1–14. [Google Scholar]
  56. The University Corporation for Atmospheric Research. Flash Flood early Warning System Reference Guide; The University Corporation for Atmospheric Research: Boulder, CO, USA, 2010; p. 204. [Google Scholar]
  57. Bhushan, N.; Rai, K. Strategic Decision Making: Applying the Analytic Hierarchy Process; Springer: London, UK, 2004. [Google Scholar]
  58. Forman, E.; Gass, S. The Analytic Hierarchy Process: An Exposition. Oper. Res. 2001, 49, 469–486. [Google Scholar] [CrossRef]
  59. Mu, E.; Pereyra-Rojas, M. Practical Decision Making: An Introduction to the Analytic Hierarchy Process (AHP) Using Super Decisions; Springer International Publishing: New York, NY, USA, 2016. [Google Scholar]
  60. Rodolfo, K.S.; Lagmay, M.A.; Eco, N.; Herrero, T.L.; Mendoza, J.E.; Minimo, L.G.; Santiago, J.T. The December 2012 Mayo River debris flow triggered by Super Typhoon Bopha in Mindanao, Philippines: Lessons learned and questions raised. Nat. Hazards Earth Syst. Sci. 2016, 16, 2683–2695. [Google Scholar] [CrossRef]
  61. Malczewski, J. GIS and Multicriteria Decision Analysis, 1st ed.; John Wiley & Sons, Inc.: New York, NY, USA, 1999. [Google Scholar]
  62. UN Office for Disaster Risk Reduction. What Does Vulnerability Mean? Available online: https://www.unisdr.org/2004/campaign/booklet-eng/Pagina8ing.pdf (accessed on 20 September 2019).
  63. Thieken, A.; Cammerer, H.; Dobler, C.; Lammel, J.; Schöberl, F. Estimating changes in flood risks and benefits of non-structural adaptation strategies—A case study from Tyrol, Austria, Mitig. Adapt. Strateg. Glob. Chang. 2016, 21, 343–376. [Google Scholar] [CrossRef]
  64. De Moel, H.; Van Alphen, J.; Aerts, J. Flood maps in Europe—Methods, availability and use. Nat. Hazards Earth Syst. Sci. 2009, 9, 289–301. [Google Scholar] [CrossRef]
  65. Nations United Affairs Department of Humanitarian. Mitigating Natural Disasters: Phenomena, Effects and Options—A Manual for Policy Makers and Planners; United Nations: New York, NY, USA, 1991. [Google Scholar]
  66. United Nations Department of Humanitarian Affairs. Internationally Agreed Glossary of Basic Terms Related to Disaster Management, DNA/93/36; United Nations: New York, NY, USA, 1992. [Google Scholar]
  67. World Meteorological Organization. Comprehensive Risk Assessment for Natural Hazards, WMO/TD No 955; World Meteorological Organization: Geneva, Switzerland, 1999. [Google Scholar]
Figure 1. Study area of Davao Oriental, Philippines, and its DEM dataset. (A) Map of the Philippines, (B) map of Mindanao where Davao Oriental is located, (C) Davao Oriental map with National Climatic Data Center (NCDC) precipitation data points and (D) Davao Oriental Map with its DEM and two weather stations (DOST-RXI and Hinatuan), and 70 field survey points at the historical flooded areas due to Typhoon Haiyan in 2013. The elevations of 70 field survey points flooded due to the typhoon are all between 10 m and 20 m, and thus the area below 20 m above mean sea level are considered as flood-prone (1: red zones) and the other area as non-flood-prone (0: grey zone) in Figure 1C.
Figure 1. Study area of Davao Oriental, Philippines, and its DEM dataset. (A) Map of the Philippines, (B) map of Mindanao where Davao Oriental is located, (C) Davao Oriental map with National Climatic Data Center (NCDC) precipitation data points and (D) Davao Oriental Map with its DEM and two weather stations (DOST-RXI and Hinatuan), and 70 field survey points at the historical flooded areas due to Typhoon Haiyan in 2013. The elevations of 70 field survey points flooded due to the typhoon are all between 10 m and 20 m, and thus the area below 20 m above mean sea level are considered as flood-prone (1: red zones) and the other area as non-flood-prone (0: grey zone) in Figure 1C.
Water 11 02203 g001
Figure 2. Comparisons of the observed annual rainfalls at the (A) Hinatuan and (B) DOST-RXI Stations with the NCDC and Global Precipitation Climatology Centre (GPCC) rainfall data at DO63 and DO20 points, respectively. Annual maximums of 5-day continuous rainfall at DO63 and DO20 are also presented in green. The locations of DO20 and DO63 are 7.083° N, 125.95° E, and 8.00° N, 126.33° E, respectively.
Figure 2. Comparisons of the observed annual rainfalls at the (A) Hinatuan and (B) DOST-RXI Stations with the NCDC and Global Precipitation Climatology Centre (GPCC) rainfall data at DO63 and DO20 points, respectively. Annual maximums of 5-day continuous rainfall at DO63 and DO20 are also presented in green. The locations of DO20 and DO63 are 7.083° N, 125.95° E, and 8.00° N, 126.33° E, respectively.
Water 11 02203 g002
Figure 3. Flowchart of multi-criteria data analysis for flood risk assessment. I1-I6 mean the six indicators and C1–C2 are the two criteria considered in the analysis.
Figure 3. Flowchart of multi-criteria data analysis for flood risk assessment. I1-I6 mean the six indicators and C1–C2 are the two criteria considered in the analysis.
Water 11 02203 g003
Figure 4. Multi-source datasets for Davao Oriental: (A) 25-yr averaged (1990–2015) 5-day continuous rainfall distribution, (B) slope, (C) elevation, (D) soil type, (E) drainage density and distance to the main channel and (F) Population density. (G) The resulting flood hazard map by analytic hierarchy process (AHP) method, and (H) flood risk map are also shown.
Figure 4. Multi-source datasets for Davao Oriental: (A) 25-yr averaged (1990–2015) 5-day continuous rainfall distribution, (B) slope, (C) elevation, (D) soil type, (E) drainage density and distance to the main channel and (F) Population density. (G) The resulting flood hazard map by analytic hierarchy process (AHP) method, and (H) flood risk map are also shown.
Water 11 02203 g004
Figure 5. Results of accuracy assessment: (A) Comparison of the accuracy assessment for the three methods using accuracy (ACC), true negative rate (TNR) and true positive rate (TPR). (B) Model evaluation using TPR and false positive rate (FPR). Redline show the line of discrimination or random guess. The upper part of redline indicates the model is accurate, otherwise not accurate.
Figure 5. Results of accuracy assessment: (A) Comparison of the accuracy assessment for the three methods using accuracy (ACC), true negative rate (TNR) and true positive rate (TPR). (B) Model evaluation using TPR and false positive rate (FPR). Redline show the line of discrimination or random guess. The upper part of redline indicates the model is accurate, otherwise not accurate.
Water 11 02203 g005
Figure 6. Flood hazard maps obtained by AHP, WR and RW methods.
Figure 6. Flood hazard maps obtained by AHP, WR and RW methods.
Water 11 02203 g006
Figure 7. AHP-based flood hazard maps with different CR scenarios, (A) CR was 2.7%, (B) CR was 5.8% and (C) CR was 8.6%.
Figure 7. AHP-based flood hazard maps with different CR scenarios, (A) CR was 2.7%, (B) CR was 5.8% and (C) CR was 8.6%.
Water 11 02203 g007
Table 1. Census of population and housing by municipality/city [48].
Table 1. Census of population and housing by municipality/city [48].
MunicipalityDistrictPopulationArea (km2)Density (/km2)No. of Barangay.
Ratio in Total (%)20102015Annual Growth Rate (%)
Baganga1st10.153,42656,2410.98945.505918
Banaybanay2nd7.439,12141,1170.95408.5210014
Boston1st2.412,67013,5351.27357.03388
Caraga1st7.236,91240,3791.72642.706317
Cateel1st7.338,57940,7041.03545.567516
Gov. Gen.2nd9.950,37255,1091.73365.7515020
Lupon2nd11.861,72365,7851.22886.397421
Manay1st7.640,57742,6900.97418.3610017
Mati City2nd25.3126,143141,1412.16588.6324026
San Isidro2nd6.432,42436,0322.03220.4416016
Tarragona1st4.725,67126,2250.41300.768710
Total 517,618558,9581.475679.6498183
Table 2. Pairwise comparison matrix based on the judgments of local experts.
Table 2. Pairwise comparison matrix based on the judgments of local experts.
IndicatorsRSlEDcDdSt
Rainfall (R)1.002.004.005.006.007.00
Slope (Sl)0.501.002.003.004.005.00
Elevation (E)0.250.501.002.003.004.00
Distance to main channel (Dc)0.200.330.501.002.003.00
Drainage (Dd)0.170.250.330.501.002.00
Soil type (St)0.140.200.250.330.501.00
Sum2.264.288.0811.8316.5022.00
Table 3. Normalized comparison matrix with priority vectors (weights). Note that the values in the tables are rounded off at 3 digits from decimal points.
Table 3. Normalized comparison matrix with priority vectors (weights). Note that the values in the tables are rounded off at 3 digits from decimal points.
IndicatorsRSlEDcDdStSum of RowsPV
Rainfall (R)0.440.470.490.420.360.322.510.418
Slope (Sl)0.220.230.250.250.240.231.430.238
Elevation (E)0.110.120.120.170.180.180.880.147
Distance to main channel (Dc)0.090.080.060.080.120.140.570.095
Drainage (Dd)0.070.060.040.040.060.090.370.061
Soil type (St)0.060.050.030.030.030.050.240.041
Sum1.001.001.001.001.001.00-1.000
Table 4. The results of the weights by rank (WR) and ratio weighting (RW) derived from the ranking of the indicators from the same local experts considered in the AHP method. The W and AW are the resulting weights of each indicator used in the calculations of hazard index in Equations (7) and (10) (refer to Tables S9 and S10 in the Supplementary Material for derivation procedure).
Table 4. The results of the weights by rank (WR) and ratio weighting (RW) derived from the ranking of the indicators from the same local experts considered in the AHP method. The W and AW are the resulting weights of each indicator used in the calculations of hazard index in Equations (7) and (10) (refer to Tables S9 and S10 in the Supplementary Material for derivation procedure).
IndicatorsWeights by Rank (WR)Ratio Weighting (RW)
ARWGMWAW
Rainfall (R)1.50.26[2.99,1.98,1.98,2.99][0.38,0.27,0.25,0.41]0.33
Slope (Sl)1.750.25[1.98,2.99,2.99,1.47][0.25,0.40,0.38,0.20]0.31
Elevation (E)2.750.20[1.26,1.26,1.26,1.26][0.16,0.17,0.16,0.17]0.17
Distance to main channel (Dc)4.750.11[0.79,0.44,0.79,0.79][0.10,0.06,0.10,0.11]0.09
Drainage (Dd)4.750.11[0.51,0.40,0.51,0.51][0.06,0.05,0.04,0.05]0.06
Soil type (St)5.50.07[0.33,0.33,0.33,0.33][0.04,0.05,0.04,0.05]0.04
Sum-1.00[7.87,7.40,7.87,7.35][1,1,1,1]1.00
Table 5. Confusion matrix to calculate the accuracy (ACC), true positive rate (TPR), true negative rate (TNR) and false positive rate (FPR) for performance evaluation of three methods, AHP, WR and RW.
Table 5. Confusion matrix to calculate the accuracy (ACC), true positive rate (TPR), true negative rate (TNR) and false positive rate (FPR) for performance evaluation of three methods, AHP, WR and RW.
ObservedPredicted
FloodNon-floodSum
Floodtrue positive (TP)false negative (FN)P
Non-floodfalse positive (FP)true negative (TN)N
AHP
Flood228,388298,332526,720
Non-flood358,439701,2471,059,686
RW
Flood197,607329,113526,720
Non-flood719,362340,3241,059,686
WR
Flood209,746316,974526,720
Non-flood592,507467,1791059,686
Table 6. Flood risk assessment (affected population and the number of affected barangay) in each municipality by the AHP method.
Table 6. Flood risk assessment (affected population and the number of affected barangay) in each municipality by the AHP method.
MunicipalityPopulation AffectedNo. of Barangay.MunicipalityPopulation AffectedNo. of Barangay
Boston13,5388Tarragona83085
Cateel40,70416Mati20,2765
Baganga51,83716Lupon86603
Caraga47044Banaybanay62523
Manay13,7554---
Total124,53848Total43,49616
Table 7. Hazard index weights of the six indicators according to CR values in the AHP method.
Table 7. Hazard index weights of the six indicators according to CR values in the AHP method.
CR (%)RainfallSlopeElevationDistanceDrainageSoil Type
8.6039.3223.7717.2710.445.833.36
5.8045.7325.5413.458.034.322.94
2.0042.1722.4614.9210.216.154.09
Table 8. Ranking of flood hazard classifications according to the CR values in the AHP method.
Table 8. Ranking of flood hazard classifications according to the CR values in the AHP method.
Classifications CR 2.0CR 5.8CR 8.6Average Change (%)
%Rank%Rank%Rank
Very Low0.6840.9841.44−0.36
Low22.23223.27222.2320
Moderate54.93155.92152.5211.21
High22.15319.84323.683−0.77
Total100-100-100-average: 0.02

Share and Cite

MDPI and ACS Style

Cabrera, J.S.; Lee, H.S. Flood-Prone Area Assessment Using GIS-Based Multi-Criteria Analysis: A Case Study in Davao Oriental, Philippines. Water 2019, 11, 2203. https://doi.org/10.3390/w11112203

AMA Style

Cabrera JS, Lee HS. Flood-Prone Area Assessment Using GIS-Based Multi-Criteria Analysis: A Case Study in Davao Oriental, Philippines. Water. 2019; 11(11):2203. https://doi.org/10.3390/w11112203

Chicago/Turabian Style

Cabrera, Jonathan Salar, and Han Soo Lee. 2019. "Flood-Prone Area Assessment Using GIS-Based Multi-Criteria Analysis: A Case Study in Davao Oriental, Philippines" Water 11, no. 11: 2203. https://doi.org/10.3390/w11112203

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop