1. Data availability and completeness
To compare different methods, we used open access data on the road network and health facilities, demographic, and administrative data as well as remote sensing-based flood masks. For the raster-based accessibility analysis, a friction surface representing the travel time per cell of a tessellation was provided by the Malaria Atlas Project (40). It was calculated as a function of street presence and street characteristics, land use, and terrain.
The flood extents were provided by two different datasets which focused on complementary areas. Both were combined in a vector layer and used as the flood mask for the following analysis steps. To lower the computational burden, the geometries were simplified using the Visvalingam algorithm (41) with only one percent of points retained. The first flood extent dataset was based on results published by the World Food Program (42) on the 21st of March, four days after the end of the cyclone. The flood extent was retrieved from a 22 km resolution microwave radiation satellite imagery from SSM/I, AMSR2, and GMI sensors and downscaled with a DEM to a 90 m resolution raster dataset (43). The images were acquired on a daily basis between the 12th and the 20th of March 2019 and aggregate the extent of the detected floods. The second dataset was published in vector format on the 26th of March 2019 by UNOSAT (44) and was acquired by the Sentinel-1 radar sensor with a resolution of 10 m. The product presents the cumulative floods detected on the 13th, 14th, 19th, 20th, and 26th of March. The whole analysis was restricted to the extent of Mozambique, using the administrative boundary from geoboundaries.org (45).
To assess quantitatively the vulnerabilities of the population regarding healthcare accessibility, we extracted population counts from a raster dataset by worldpop (46). Worldpop (47) provides several datasets of gridded population counts with a 100 m resolution. We used the 2020 constrained top-down method adjusted with UN population estimates for Mozambique.
Road network and health facility data were obtained from OpenStreetMap. As the completeness of the dataset is dependent on volunteer contributions and increases over time, we used the latest extract available at the moment of the analysis in April 2021. The road network was extracted from Geofabrik.de (48). For the analysis of the flood event, we used the tool Osmosis (49) to simulate the road network alteration due to the flood. All nodes of the road network that were intersecting with the flood mask were removed. The healthcare facilities were obtained via the ohsome API (50). We considered 1059 healthcare facilities recorded in OpenStreetMap for Mozambique, including both primary and non-primary healthcare facilities.
The acquisition of the OSM dataset at the time of the analysis instead of the time of the event represents a potential bias. To assess this bias and verify the data quality, we performed an intrinsic data quality approach. We accessed the OSM history to assess the evolution of the dataset over time on relevant features. We used Ohsome API (50) to filter the contributions on specific highways (highways = motorway, trunk, primary, secondary, tertiary, unclassified, track, path) and all health facilities that can provide primary health care (amenity=clinic, health_post, doctors, healthcare=doctors, clinic, health_post, midwife, nurse, center) including also hospitals (amenity=hospital, healthcare=hospital, building=hospital). Additionally, we analyzed the evolution of user activities on these features. A sudden increase in the mapped features in the flooded area in the wake of the cyclone would give an insight into the lack of data completeness at this time and the activation of remote mapping campaigns to address it. On the contrary, a sudden decrease in user activities would illustrate a drop in the interest in mapping the area.
Moreover, we used another indicator to estimate the completeness of OSM features mapped through remote mapping campaigns such as roads and buildings (51). We compared the density of these features in OSM and population density from an external dataset (worldpop) and used it as an estimation of feature completeness. The correlation between roads and population density is not straightforward, thus we used this indicator on building completeness as a proxy for road completeness. High scores of the indicators imply that feature density is well aligned with population density which was taken as an indication of a relatively high completeness. On the contrary, a low score was interpreted as an indication that many features were not mapped and thus taken as an indication for a low feature completeness. The indicator was calculated for the whole country using a regular tessellation of 346 hexagons. Each hexagon covered 2903 km2. We used the Ohsome Quality Analyst (52) to calculate the indicator, taking a snapshot of the OSM database from the 20th of September 2021.
Furthermore, we interviewed a humanitarian field worker who participated in the emergency response in the wake of the cyclone in the area of Dondo, Sofala Province. Both his observations about health facilities and the GPS coordinates collected by him were used as auxiliary data to judge the quality of the OSM data. The field worker visited seven health facilities in the Dondo District which were part of the study region. His information allowed us to assess whether all the health facilities were mapped in OSM in the sample region and if their location was accurate. A lack of location precision would affect the analysis as it focused on the accessibility and the criticality of the roads to access these sites.
2. Healthcare accessibility analysis
Our analysis estimated the potential access to health facilities, the loss of accessibility and the impacted population by a loss of accessibility. The potential access to healthcare analysis was run independently in both flooded and normal conditions. Both were then compared based on their geometrical difference (see Additional file 1: Figure S1 for detailed workflows). Thereby, we identified the areas which experienced a loss of accessibility to health facilities. Using the worldpop dataset, we finally estimated the number of inhabitants by access time under both conditions and in the areas experiencing a loss of accessibility. The results allowed us to measure the impact of the floods on the accessibility of the population to health facilities.
The healthcare accessibility analysis was done with both network- and raster-based methods and their results were compared to identify the most appropriate approach depending on contextual characteristics. The literature on healthcare accessibility tends to focus on a transportation mode based on motorized vehicles. Nevertheless, in many Sub-Saharan countries, including Mozambique (53), walking remains the main means of transport especially in disaster situations and in rural areas. Therefore, we considered both transportation profiles.
The estimation of accessibility with a network-based approach was based on isochrones (54) which represent the area accessible from one point within a given time. The isochrones were calculated for the 1059 health centers of Mozambique using the openrouteservice API (55). Openrouteservice uses OSM road network data to build a weighted graph according to road categories, speed limits, and traffic restrictions. Hence, the calculation was restricted to road segments and did not consider off-road travel. The isochrones calculated were based on an interval of time to reach the health facility. We used 10 minute intervals capped to a maximum of one hour for the driving-only profile. For the walking-only profile, we also considered 10 minute intervals which were then aggregated by hour to a maximum of six hours. As the loss of accessibility calculations rely on the difference in travel times and the geometric differences between the respective isochrone polygons, it resulted in potential precision error. This potential error was equal to the time interval (10 minutes).
The raster-based method calculates a least-cost path on every pixel of the raster grid to access the same healthcare facilities. The algorithm is based on a friction layer that attributes a time cost to cross a pixel. We used a walking-only friction layer as well as a friction layer combining walking with motorized transportation where road and train information is available. Contrary to the network-based method, the calculation was not restricted to the road network and considered travel on any surface with an adjusted speed. The gDistance (see Additional file 1: Table S1) R package was used to compute the least-cost path. In this approach, we altered the cost of the flooded areas to NA values to forbid travel. The continuous result was categorized into time intervals (10 minutes) to allow a comparison with the network-based approach.
Populations experiencing a loss of potential access may be in the close surroundings of the floods or further away. In both cases, it would highlight a lack of resilience of one or several road segments which inhibit access to health facilities along the shortest paths. A new network centrality analysis allows for identifying such vulnerable axes.
3. Network Criticality analysis
The analysis of the road network criticality requires a graph representation of the network. We used openrouteservice (55) to export the weighted graph. We developed a novel criticality indicator, the targeted edge centrality (TEC) indicator, as an extension of edge betweenness centrality (56). Instead of considering the shortest paths between all possible pairs of nodes in the network (as for edge betweenness centrality), the TEC calculates the road segment usage frequency among the shortest paths between targeted destination points (in our application healthcare facilities) and departure points of the area of interest (equation 1). The departure points were set as population nuclei and were computed from the centroid of the cells of the worldpop dataset after aggregation to a 1 km grid. For each population nucleus, the targeted destination points were selected within a buffer zone of 20 km. This filter implied the assumption that a person has an uniform probability to visit any of the destination points within the buffer zone. To avoid any bias due to the heterogeneity of the destination density, the frequency was weighted by the number of destinations accessible from each departure point.
Segments with higher TEC scores are more frequently crossed. A similar TEC score for all segments within a network indicates a lower network criticality level as there are many alternatives to any road disruptions. On the contrary, a wider range of TEC scores indicates a stronger vulnerability of the network as the loss of one main axis could isolate parts of the network.
The TEC was calculated for the road network for both normal and flooded conditions. For the flooded case, road segments inside the flood masks were not considered. The scores for both cases were then compared. For a single road segment, a decrease in the TEC score after flooding indicates higher criticality as it is not able to provide access due to the floods. In a network with a hub-and-spoke topology (57), these criticalities could result in a lack of resilience if no mitigation action is taken, such as flood walls. On the other hand, an increase in the score would identify the segment as an important backup road. Such roads should get special attention in disaster preparedness plans as they have major importance in disasters.
4. Software and data
The data processing was done with R 4.0.5 (58). The scripts and documentation are available in open access. In addition, the complete list of used R packages is available in the Additional file 1 : Table S1. Osmosis was used to remove the OpenStreetMap data within the extent of the floods. Openrouteservice was installed locally using Docker (59). The maps were created with QGIS 3.18 (60).