The procedure of developing the urban 3D terrain is shown on the left side of Figure 3. The methods of generating DEM from a set of aerial images or videos are quite mature (Zhou et al., 2004; Pollefeys et al., 2008), which are based on the fundamental principle of collinearity condition (Figure 4) expressed by the following equations: x p = − f [ m 11 ( X P − X L ) + m 12 ( Y P − Y L ) + m 13 ( Z P − Z L ) m 31 ( X P − X L ) + m 32 ( Y P − Y L ) + m 33 ( Z P − Z L )   ] (1) y p = − f [ m 21 ( X P − X L ) + m 22 ( Y P − Y L ) + m 23 ( Z P − Z L ) m 31 ( X P − X L ) + m 32 ( Y P − Y L ) + m 33 ( Z P − Z L ) ] , (2) where x p and y p are the image coordinate of any point p ; X P ,   Y P ,  and  Z P represent the ground coordinate of point p; X L ,   Y L ,  and  Z L represent ground coordinate of the projection (optical) center; f is the focal length; m 11 … m 33 are the coefficients of a 3 × 3 rotation matrix defined by the angles ω , ϕ , and κ that transforms the ground coordinate system to the image coordinate system (Lillesand and Kiefer, 1999). The six parameters X L ,   Y L ,  and  Z L and ω , ϕ , and κ are used for exterior orientation of an image, which can be determined through the process of space resection, as illustrated in Figure 4. In the figure, the X , Y , and Z coordinate of any point in the matched stereo pair of tilted images can be determined. Using a set of images taken by UAV, the ground coordinate of any point in the overlap of tilted images can be determined by the image matching method of structure from motion (Remondino and Fraser, 2006; Westoby et al., 2012). Finally, the coordinate accuracies at checkpoints are examined, and the urban 3D terrain model is obtained by resampling the coordinates to a regular gird system by the nearest neighboring method (Wu, 2013).

The UAV used in this study is DJI Phantom 2 Vision+ (Da-Jiang Innovations) which weighs 1.2 kg and has a camera with 4,384 × 2,466 pixels. The focal length of the camera is 3.3 mm, and the field-of-view is 110°. The UAV campaign was conducted in the early morning of July 22, 2015, in a cloudless condition from 06:00 to 06:40 to reduce the disturbances from weather and traffic. The camera on the UAV is set in a nadir-looking orientation for image data collection. In total, 589 positioned images were acquired with an overlap ratio of 75–85% and a mean spatial resolution of 2.84 cm. After space intersection, the average ground sampling distance of point cloud is 0.03 m. The UAV images were processed to generate orthomosaic image and digital surface model (DSM) with Pix4Dmapper Pro version 1.4.46 (Pix4D). To derive absolute coordinates, three ground control points (GCPs) were distributed in the study area (Figure 5). The coordinates of the three GCPs were acquired using the static positioning of the global navigation satellite system (GNSS) with positional accuracy at centimeter level. The absolute positions of the images captured through the GNSS receiver in the UAV were recorded to establish the coordinate system constrained on the three GCPs.

The lens distortion of the camera was calibrated by flexible and powerful self-calibrating bundle adjustment (Remondino and Fraser, 2006). The calibration relies on the continuous overlapped images to do the aerotriangulation adjustment. After the interior orientation and the coordinates of the object points are calculated, the correction Δ x and Δ y can be revised as follows: { Δ x = − Δ x 0 − x i f Δ f + K 1 x p r 2 + K 2 x p r 4 + K 3 x p r 6 + P 1 ( r 2 + 2 x p 2 ) + 2 P 2 x p y p Δ y = − Δ y 0 − y i f Δ f + K 1 y p r 2 + K 2 y p r 4 + K 3 y p r 6 + P 2 ( r 2 + 2 y p 2 ) + 2 P 1 x p y p , (3) where Δ f is the principal distance error; Δ x 0 and Δ y 0 are the displacements of the principal point; K 1 , K 2 , and K 3 are the parameters of the radial distortion; P 1 and P 2 are the parameters of the decentering distortion; and r is the distance between the image point and the principal point. The calibrated parameters are listed in Table 1.

For CFS application, the DSM was converted to DEM by removing the vegetation and the viaduct. The vegetation such as shrubs and grasses is detected by the normalized difference vegetation index (NDVI) in the range of 0.2–0.3 (Candiago et al., 2015). Since the UAV images only observed red, green, and blue bands, the near-infrared band was built on a specific linear combination of the three bands with a lower-pass filter (Rabatel et al., 2011). To remove the viaduct, the DSM points with elevations higher than 9 m on the Keelung Road are removed, and the DEMs beneath the viaduct were obtained through interpolation using the neighboring ground elevation data obtained by the UAV. This treatment did not deteriorate flood simulation because the area beneath the viaduct occupied about only 10% of the entire study domain, and the viaduct is located along the road centerline with higher ground elevations where floodwater less accumulated. Despite the viaduct, all the buildings are maintained so that flood simulation can more truly reflect the field conditions. Based on the NDVI and elevation thresholds, the vegetation and viaducts were filtered out in the DEM so that floodwater can transport smoothly on the ground surface. However, unlike traditional DEM, the elevations of buildings were not removed to reflect the blocking effect on water transverse.

Ethical and legal concerns are big issues for collecting and using VGI (Foody et al., 2017). Fortunately, the Copyright Act and Relative Laws in Taiwan allow researchers to quote, within a reasonable scope, publicly released works for reports, comment, teaching, research, or other legitimate purposes (Intellectual Property Office., 2019). Based upon the act, the VGI data used in this study are collected from the most famous bulletin board system (BBS) in Taiwan named PTT. There are 8 photos collected from PTT posted during 15:20–16:30 on June 14, 2015. From these photos, we visually identified 8 locations in the study area, as shown in Figure 6. Since the VGI photos downloaded from social media were without EXIF attributes, an obvious landmark in each photo was identified to mark the GNSS and visually determine the relative flood depth with reference to neighboring objects with generally fixed size such as wheels of bikes or elevations of sidewalks. This GNSS information can be obtained based on post-disaster field surveys and is invariant as long as the landmarks were not destroyed during flooding. The time stamp and the virtual water depths in these photos were used to validate the CFS model. Although the time stamp when photos were posted on the internet may not always be the acquisition time and the flood depth estimated from photos may be subject to experts’ experience, these uncertainties can be reduced by the functions of “live stream” and “image recognition” on social media.

The CFS model used in this study was developed by Jang et al. (2018), in which a 2D overland flow model (OFM) is coupled with a 1D sewer flow model (SFM) for sophisticated flood simulation in cities with drainage systems. The OFM and SFM are established based on shallow water equations and finite difference numerical methods. The alternate direction explicit scheme and implicit backward Euler algorithm are used to solve the OFM and SFM, respectively. The SFM adopted the Preissmann slot method (Cunge and Wegner, 1964) to calculate the full and partially full flow conditions at the same time. For the OFM simulation, the elevation at each grid center is extracted from the DEM data. During rain, the OFM is first initiated for surface water routing and the SFM is then initiated by the water that flows into the sewer pipes via street inlets. When the sewer pipes get full, the sewer water surcharges back onto the ground surface via manholes. In the simulation process, the water exchanged between the two models is determined by weir and orifice functions via a one-to-one relationship. To deal with wetting and drying, a tolerance water depth of 0.001 m is used to distinguish wet and dry cells. If a cell has a water depth lower than 0.001, all velocities in this cell are set to zero to ensure numerical stability caused by sudden drying in the next time step. The CFS model has been validated in urban and coastal areas with an overall hit rate of 0.8, meaning that 80% of the observed food extents are correctly predicted. The details of the CFS model can be referred to Jang et al. (2018); Jang et al. (2019).

The DEM accuracy is examined at the six GCPs. The errors in X , Y , and Z directions range from −0.027 to 0.025 m, −0.055–0.061 m, and −0.024–0.026 m, respectively. The root mean square errors (RMSEs) of the absolute coordinate at the three GCPs are less than 0.053 m (Table 2). The RMSEs of SRTM were less than 10.343 m, and the RMSEs of 20 m DEM were less than 3.953 m. The comparison of the RMSE values showed that the UAV terrain model was with the smallest RMSE to achieve high precision topography. The orthoimage and the DEM resampled under spatial resolutions of 0.5, 5, 10, and 20 m are shown in Figure 7. In the figure, the elevation distributes from 5.5 to 55.8 m, and the buildings are displayed in warm colors.

To discover the influence of DEM resolution on flood simulation, the grid meshes of the CFS model are established under four grid sizes ( 0.5  m × 0.5  m , 5  m × 5  m , 10  m × 10  m , 20  m × 20  m ), in which the elevation at each grid center is extracted from the DEM with accordant resolution. The flood extents simulated under the four DEM resolutions at different times are displayed in Figure 8, in which the VGI points out of the 8 locations are marked if the simulated flood depths exceed 0.05 m. The value of 0.05 m is selected as the threshold because it is the shallowest water depth observed in the VGI photos (see point #6 in Figure 6). In the case with 0.5 m DEM resolution, flooding starts around 14:00 at points #4, #7, and #8, peaks around 15:00 at all points, and retreats at 17:00 back to original points. With coarser DEM resolutions, flooding occurs earlier and retreats later with more water trapped between buildings. When the grid resolutions increase to 10 and 20 m, the flood areas become more discontinuous with staggered distribution of wet and dry cells, and the flooding at points #1, #3, #5, and #6 are not predicted. Observing the flood maps between 15:00 and 18:00, water inundated on the rooftops of buildings can be properly simulated in the case with 0.5 m DEM resolution but not in the case with coarser resolutions.

The time series of flood volume simulated under the four DEM resolutions are displayed in Figure 9. Compared with the case with 0.5 m DEM resolution, the flood volume simulated under coarser DEM resolutions arises faster but descends slower with higher and earlier appearance of flood peak. This implies that when DEM resolution decreases, the topography becomes rugged, the friction increases, and the floodwater travels slower. The flood volumes for 10 and 20 m DEM resolutions are much larger than those for 0.5 and 5 m DEM resolutions, indicating that more water will be trapped in the depressions when the topography becomes more rugged. Table 3 compares the simulation and observation results in terms of flood occurrence time and water depth. The time stamps and estimated water depths (WD) are obtained from the VGI photos in Figure 6, and the flood durations at the eight VGI points when the water depth exceeds 0.05 m are determined based on the CFS results. It is seen that under 0.5 and 5 m DEM resolutions, the time stamps of VGI photos all lie within the simulated flood duration at the points with observed WD larger than 0.05 m (points #1, #2, #4, #7, and #8). At the rest points, the simulated and observed WDs are both smaller than 0.05 m. This good agreement between observation and simulation reveals that the flood model is accurate in rebuilding the process of flood transport under both DEM resolutions.

When the DEM resolutions become larger than 5 m, the periods with simulated WDs larger than 0.05 m do not coincide with the time stamps of VGI photos for 10 m DEM resolutions (at points #1 and #8) and 20 m DEM resolutions (at points #1, #2, and #8). Thus, the spatial resolution of DEM in urban areas should be at least finer than road width so that road profiles can be clearly displayed; otherwise, runoff transportation around buildings and on roads cannot be correctly simulated. The resolution of 5 m can be regarded as a threshold for flood modeling in urban areas, above which the geometries of buildings and roads cannot be clearly identified and the simulation performance deteriorates. Sub-meter resolution DEM data generated by UAVs are adequate for assessing the impact of localized flooding on transportation in cities. Some open-access topographic datasets, such as the 30 m resolution DEM by STRM from NASA (https://www2.jpl.nasa.gov/srtm/) and the open DEM in Taiwan with 20 m resolution (https://data.gov.tw/dataset/35430), are too coarse to serve the purpose.

Flood simulation under high grid resolutions is usually more time-consuming as compared to that under coarse grid resolutions due to the increase in grid numbers. The choice of grid resolution for flood simulation is a trade-off between accuracy and efficiency. For the study area in this research, there are 573,000 and 5,730 grids for the mesh with grid size 0.5 m and the mesh with grid size 5 m, which consumes 1,127 and 16 min of computational time, respectively (with Intel Core i7-7700 K CPU @ 4.2GHz and 64 GB RAM). For disaster emergency response on a regional scale, flood simulation under coarse grid resolution is enough to gain a fast and overall understanding of flood patterns. However, for evaluating the flood impact on critical infrastructures such as metro stations, power facilities, schools, government agencies, and hospitals, high-resolution flood simulation is required.