A Visualization Tool for Flood Dynamics Monitoring Using a Graph-Based Approach

Insights into flood dynamics, rather than solely flood extent, are critical for effective flood disaster management, in particular in the context of emergency relief and damage assessment. Although flood dynamics provide insight in the spatio-temporal behaviour of a flood event, to date operational visualization tools are scarce or even non-existent. In this letter, we distil a flood dynamics map from a radar satellite image time series (SITS). For this, we have upscaled and refined an existing design that was originally developed on a small area, describing flood dynamics using an object-based approach and a graph-based representation. Two case studies are used to demonstrate the operational value of this method by visualizing flood dynamics which are not visible on regular flood extent maps. Delineated water bodies are grouped into graphs according to their spatial overlap on consecutive timesteps. Differences in area and backscatter are used to quantify the amount of variation, resulting in a global variation map and a temporal profile for each water body, visually describing the evolution of the backscatter and number of polygons that make up the water body. The process of upscaling led us to applying a different water delineation approach, a different way of ensuring the minimal mapping unit and an increased code efficiency. The framework delivers a new way of visualizing floods, which is straightforward and efficient. Produced global variation maps can be applied in a context of data assimilation and disaster impact management.


Introduction
The emergence of several new Synthetic Aperture Radar (SAR) missions throughout the past decade has boosted the field of flood remote sensing. However, most studies focus on the development of new or improved flood mapping algorithms, which result in single date flood extents. These products fail to provide information on flood dynamics, which is crucial for several applications. For example, flood duration can be linked directly to the degree of economic losses and damages to structures [1][2][3][4] in the context of disaster impact assessment (emergency relief, damage assessment, insurance claims . . . ). Moreover, data assimilation for hydrodynamic forecasting could be improved by preserving information-dense regions, while removing redundant points [5][6][7][8]. In this way it is possible to determine where, when, and how often remote sensing data should be acquired for flood extent assimilation to be most effective. [6]. With this research we aim at producing a unified flood map, not limited to just flood extent, but incorporating insights on dynamics.
In [9], a visualization/analysis framework was established for flood dynamics, combining objects and graphs. This approach was based on the research of [10,11]. In this letter, we have extended our previous work to cover larger areas. This allows for a new way of visualizing flood dynamics. In order to test the behavior of the workflow proposed in [9] at a larger scale, two floods were chosen which The second case is a flood which occurred between the 27 and 30 September 2019, caused by heavy rain in the East of India, in the State of Bihar, near the city of Patna ( Figure 2). While India's monsoon season normally starts to retreat by early September, heavy rainfall has continued here, triggering floods in many low-lying areas, resulting in waist-deep floodwaters [15]. For this case, Sentinel-1 imagery is available on eight timesteps in the period between 6 September and 30 October 2019.

Reference Data
While a ground truth is not available for these flood events, different "flood authorities" have produced flood extent maps for the regions of interest.
For the first case, Mozambique, we compare our results with flood maps from two external sources. The Luxembourg Institute of Science and Technology (LIST) provided us with a flood extent map for the Sofala region in Mozambique. They collaborate with the European Space Agency Research and Service Support Team (ESA RSS) for the production of automated Sentinel-1 based flood delineation products. The flood map is obtained by applying an automatic unsupervised change detection algorithm to pairs of Sentinel-1 images (pre-flood and during flood) [16]. The second source of flood maps for Mozambique is Disaster Charter [17], an initiative from ESA, The French National Centre for Space Studies (CNES) and The Canadian Space Agency (CSA), harvesting Earth observation assets from different space agencies. They provide flood maps at three timesteps for the Mozambique flood, based on Sentinel-1 imagery between 13 and 26 March 2019 over the Sofala province and its surroundings.
The second case, India, was compared with one external flood map, originating from Sentinel-Asia [18], a voluntary basis initiative led by the Asia-Pacific Regional Space Agency Forum (APRSAF) to support disaster management activity in the Asia-Pacific region. They provide a single-date flood map over the region of interest with a spatial resolution of 30 m × 30 m. The map is derived from SAR data acquired by Sentinel-1 satellites before (21 May 2019) and during (30 September 2019) the event.

Method/Method Adaptations
Most of the workflow, presented in [9], has been retained for upscaling this research, yet some modifications were still necessary ( Figure 3).
The Active contour modelling (ACM) sufficed for delineating water bodies in a small area [9], yet application of ACM on the current regions of interest leads to big over-and underestimations of the water extent. A more suitable way of obtaining the binary water maps was found by application of tiled thresholding using the Otsu algorithm [19]. This method selects a threshold by minimizing the within-class variance of two groups of pixels separated by the thresholding operator. It does not depend on modeling the probability density functions, and it assumes a bimodal distribution of gray-level values.
Permanent water bodies are then eliminated from the binary water maps by classifying pixels which are recognized as water in all timesteps, to permanent water, and excluding them from further evaluation.
In [9], the minimal mapping unit (MMU) was ensured by dropping all water polygons smaller than 16 pixels after polygonization. This however, results in a net loss of water area. For the current research, the MMU was obtained in a different way, i.e., by applying morphological operators [20]. First, a morphological closing reassigns non-water areas of maximum 4 by 4 pixels to water, if they are completely surrounded by water. Subsequently, morphological opening reassigns isolated water bodies smaller than 4 by 4 pixels non-water. In this way, the net loss of water area is compensated. At this point, objects are created ('Polygonize' in Figure 3), reducing all detected water bodies to their individual extent and mean backscatter. Subsequently, they are grouped into entities ('Graph construction' in Figure 3) if they display spatial overlap on consecutive timesteps. Their evolution is then studied, based on the evolution of area and mean backscatter. This results in a different perception of the flood compared to a pixel-based analysis. Moreover, the analysis is simplified, which improves calculation time and ease of interpretation.
The evolution of each entity o * is described by an evolution graph G o * , which is a Directed Acyclic Graph (DAG), G 0 * = (V 0 * , E 0 * ). The nodes V represent the water bodies and the edges E represent their spatio-temporal overlap (from Khiali et al. [11]).
Dynamics are quantified as follows [10] : Equation (4) describes the variation within one entity (graph) between timestep i and timestep i + 1 as the weighted sum of variations of all polygons at timestep i, where these polygon variations are defined by the difference in backscatter with overlapping polygons from the next timestep (i + 1), weighted by their amount of spatial overlap.
The global variation (GlobVar) for a graph represents the sum of the variations on each timestep (Equation (5)). The value of GlobVar reflects the temporal behaviour of the entity, where high values indicate high spatio-temporal dynamics. Variation and global variation are unitless (Equations (4) and (5)).

Pre-Processing and Thresholding
Results of pre-processing and Otsu tiled thresholding can be found in the Supplementary Materials. Since the external flood extent maps, presented in Section 2.2, are also associated with uncertainty, a statistical validation has been omitted. Only a visual comparison with external flood extent maps is executed. Based on this visual comparison, we conclude that the produced flood maps correspond well with the external flood maps and that these can be used for testing our algorithm. The results of these comparisons are only presented in the Supplementary Materials because the focus of this tool is on combining flood extent maps over a period of time and not on producing single date flood maps.

Global Variation Versus Maximum Flood Extent
Variation is calculated for each timestep transition (available in Supplementary Materials) and summed to obtain global variation values. All maps presented throughout this article, are produced using EPSG:4326, the geodetic coordinate system for the world (with Datum "World Geodetic System 1984", Ellipsoid "WGS 84" and prime meridian "Greenwich"). When

Graphs
Contruction of graphs and variation calculation takes about 42 minutes for Case 1 (Mozambique, 7 timesteps) with an Intel(R) Core(TM) i7-4770 CPU @ 3.40GHz, 3401 Mhz, 4 Core(s), 8 Logical Processor(s). It took 6h24min for the second case. This is because of the dispersed nature of the India flood, forming a multitude of graphs (32505 graphs in Case 2 versus 2192 graphs in Case 1). Because of the object-based approach, an upscaling in area does not necessarily imply an increase of calculation time. However, the calculation time can be linked to the degree of dispersion. A dispersed flood will result in more entities, meaning more graphs and a higher calculation time. The Mozambique flood reaches a very large extent at the third timestep, combining previous separate entities in one big entity (one graph).
In Figure ??, the temporal profile of one entity is depicted (an additional example can be found in the supplementary materials). Each node in the graph represents a (flood) water polygon, while an edge depicts a spatial overlap between the connecting nodes. In this way, a water body is described by the evolution of the number of polygons and their backscatter over the timesteps. In Figure ?? the last timestep of Figure ?? is shown in detail. Five polygons are present at this timestep (TS8) with

Graphs
Construction of graphs and variation calculation takes about 42 minutes for Case 1 (Mozambique, 7 timesteps) with an Intel(R) Core(TM) i7-4770 CPU @ 3.40 GHz, 3401 Mhz, 4 Core(s), 8 Logical Processor(s). It took 6h24min for the second case. This is because of the dispersed nature of the India flood, forming a multitude of graphs (32505 graphs in Case 2 versus 2192 graphs in Case 1). Because of the object-based approach, an upscaling in area does not necessarily imply an increase of calculation time. However, the calculation time can be linked to the degree of dispersion. A dispersed flood will result in more entities, meaning more graphs and a higher calculation time. The Mozambique flood reaches a very large extent at the third timestep, combining previous separate entities in one big entity (one graph).
In Figure 6, the temporal profile of one entity is depicted (an additional example can be found in the Supplementary Materials). Each node in the graph represents a (flood) water polygon, while an edge depicts a spatial overlap between the connecting nodes. In this way, a water body is described by the evolution of the number of polygons and their backscatter over the timesteps. In Figure 7 the last timestep of Figure 6

153
In Figure 6, the temporal profile of one entity is depicted (an additional example can be found in the 154 supplementary materials). Each node in the graph represents a (flood) water polygon, while an edge 155 depicts a spatial overlap between the connecting nodes. In this way, a water body is described by    Figure 6 showing the 5 different polygons and their respective mean backscatter.

Discussion
This method was adapted to be operational on a larger scale: a different water delineation approach was applied (Otsu tiled thresholding), the minimal mapping unit was ensured in a different way (by application of morphological operators) and code efficiency was improved. This results in an innovative framework to visualize floods that is not only straightforward but also efficient. Depending on the flood type (disperseness), calculation times differ between 42 min and 6 h24 min.
While the Mozambique flood occurs in a coastal zone and the India flood takes places inland, their global variation maps (Figures 4a and 5b) show similarities. The majority of the high dynamic zones are located in or near the riverbeds (of the Buzi river in Mozambique and of the Ganges in India). This is not surprising because these areas are prone to flooding if the water level in the river rises. However the high dynamic zones are not restricted to the riverbeds only. These zones can point out interesting areas in the context of damage assessment and data assimilation.
Flood duration is a parameter often used for damage assessment in an agricultural [2,3] or urban context [4], for emergency relief or insurance claims. Because flood duration fails to incorporate spatio-temporal dynamics, we propose to combine both flood duration and (global) variation. If the flood duration parameter is not available, (global) variation values could be used as an operational parameter since they reflect the flood duration. In this way, variation values can be combined with land-use maps, road networks and population density maps in the context of emergency relief.
Calculated variation values could also be linked to the sensitivity of observation locations in a data assimilation framework, thereby improving flood predictions [5][6][7]. Observation targeting is a process in which supplementary observations are assimilated to improve the analysis in selected regions, and thereby reduce the uncertainty of a (flood) model [21]. When remotely sensed observations are used for data assimilation, thinning strategies are often applied to eliminate redundant information before the data assimilation [8]. Several authors determine their observation locations (whether or not for data assimilation purposes) on the basis of variances, be it model variance or observation variance [5,21,22]. In this way, observed (global) variation, calculated by the current framework, could function as an estimation for model variation in data scarce regions.
Finally, global variation results can be combined with the maximum flood extent in order to increase the understanding of the flood dynamics. It might even be interesting to add both maps, combining maximum flood extent (binary) and the global variation map in one unified flood extent/variation map. In this way, the area component can be combined with a dynamics component. The result of this approach is displayed in Figure 8 for the India flood.

Conclusions
In this letter we have successfully distilled a flood dynamics map from a radar SITS. The resulting global variation map evaluates floods, not only in terms of extent and duration, but also in terms of spatio-temporal variation. These insights into flood dynamics can therefore serve in operational frameworks for multiple applications. For example, this method can point out locations to be monitored to improve a hydraulic model (in the context of data assimilation) or the locations in need of damage relief (in the context of disaster management). Moreover, the framework is fast, effective and simple. By using 2 exemplary flood cases, further upscaling is expected to go smoothly. Computation time increases with the number of entities, rather then with an increase in area. Since the method is generically deployable, it could be applied to different types of (natural) phenomena (e.g., snow, vegetation). Future research will include further upscaling to country/continent scale, analysis of longer timeseries, and inclusion of different sensors and platforms (e.g., Sentinel-2).