Post-Disaster Recovery Monitoring with Google Earth Engine

Post-disaster recovery is a complex process in terms of measuring its progress after a disaster and understanding its components and influencing factors. During this process, disaster planners and governments need reliable information to make decisions towards building the affected region back to normal (pre-disaster), or even improved, conditions. Hence, it is essential to use methods to understand the dynamics/variables of the post-disaster recovery process, and rapid and cost-effective data and tools to monitor the process. Google Earth Engine (GEE) provides free access to vast amounts of remote sensing (RS) data and a powerful computing environment in a cloud platform, making it an attractive tool to analyze earth surface data. In this study we assessed the suitability of GEE to analyze and track recovery. To do so, we employed GEE to assess the recovery process over a three-year period after Typhoon Haiyan, which struck Leyte island, in the Philippines, in 2013. We developed an approach to (i) generate cloud and shadow-free image composites from Landsat 7 and 8 satellite imagery and produce land cover classification data using the Random Forest method, and (ii) generate damage and recovery maps based on post-classification change analysis. The method produced land cover maps with accuracies >88%. We used the model to produce damage and three time-step recovery maps for 62 municipalities on Leyte island. The results showed that most of the municipalities had recovered after three years in terms of returning to the pre-disaster situation based on the selected land cover change analysis. However, more analysis (e.g., functional assessment) based on detailed data (e.g., land use maps) is needed to evaluate the more complex and subtle socio-economic aspects of the recovery. The study showed that GEE has good potential for monitoring the recovery process for extensive regions. However, the most important limitation is the lack of very-high-resolution RS data that are critical to assess the process in detail, in particular in complex urban environments.


Introduction
Natural disasters have devastating impacts on infrastructure, business sectors, and people in the affected region. Between 1998 and 2017, more than 5.7 billion people were affected by disasters, and more than one million people were killed by such events, while a total loss of USD 2.9 trillion was reported [1]. Recovery starts after the immediate post-disaster response phase, mainly search and rescue operations, have concluded. Post-disaster recovery is the process of reconstructing communities in all their aspects (e.g., physical, economic, social, and environmental) in order to return life, livelihoods, and the built environment to their pre-impact [2], or even better, states, as per the Sendai Framework [3]. Recovery can take years or even decades, and is the least studied phase of the disaster management high-resolution images to extract RS-based proxies. They concluded that although RS can provide useful information in addition to socio-economic survey data, the high cost of very high-resolution images and the need for high computational processing power are the limitations of utilizing RS data. Moreover, deep learning, in particular convolutional neural networks (CNN), has become the state-of-the-art method in computer vision and RS data analysis, and has also been used for post-disaster damage and recovery assessments [9,23,24]. For instance, a new image-patch generation approach was developed to train a CNN-based network simultaneously with multi-temporal satellite images and assess recovery [25]. In addition, free OpenStreetMap building footprints were employed to automatically generate training data from very high-resolution satellite images for a pixel-level deep learning method for post-disaster building database updating [26]. The study showed that the proposed approach significantly decreased the manual work of training area collection, while maintaining the accuracy of the detected damaged, reconstructed, and newly constructed buildings at a high level. However, extensive computational power is needed to execute deep learning methods.
Most of the developed RS-based approaches in post-disaster damage and recovery assessments focus on the use of costly very high-resolution data that require extensive digital storage and computing capacity to make use of them. In recent years, freely available RS images made available through platforms such as Google Earth have attracted the attention of researchers to extract useful information [27,28]. In addition, cloud-based platforms such as Google Earth Engine (GEE) provide free RS data and computing power with a coding environment to develop and implement user-defined methods and process the data [29,30]. Freely available, low-to medium-resolution images by data providers (e.g., Copernicus) are mainly collected in these platforms. The potential of such systems has been studied for vegetation cover change analysis [31], land cover and land use classification [32,33], change detection/analysis [34,35], wetland map generation [36,37], rangeland and crop monitoring [38,39], and in the DRM domain for flood prevention and emergency response [40], drought assessment [41], and wildfire progress mapping [42]. However, they have not yet been studied to monitor the post-disaster recovery process. Furthermore, there is a need for higher spatial resolution images for detailed analysis (e.g., functional assessments), especially in urban areas.
The aim of this study was to test the suitability of GEE for a large-scale post-disaster recovery assessment. The big time series Landsat 7 and 8 data available on GEE were used to first generate image composites for the before, event/just after, and post-disaster times (i.e., three time-steps after the disaster) to obtain cloud-free images, by adapting the method developed by De Alban et al. [43]. Then, the Random Forest classifier available on the cloud computing platform was employed to classify the generated images and monitor the land cover changes during the recovery process. The entire Leyte island (with an area of ca. 7300 km 2 ) in the Philippines was selected to test the developed approach, and municipality-based damage and recovery maps for the entire island were generated for the selected time epochs. Leyte region was hit by Super Typhoon Haiyan on 8 November 2013. As one of the strongest typhoons on record worldwide, it resulted in more than 6000 fatalities, and total damages were estimated at ca. USD 2.2 billion [44].

Case Study and Google Earth Engine Data
Leyte province, located in the Eastern Visayas, is the seventh largest island in the Philippines (Figure 1). With a population of approximately 250,000, Tacloban is the biggest city and the capital of Leyte, which in total houses ca. 1.7 million people. The region has a tropical climate with two slightly different rainfall patterns during a year. While most of Leyte island has an even annual rainfall distribution, some eastern parts have a pronounced maximum rainfall from November to January. Typhoon Haiyan (also known as Typhoon Yolanda in the region) passed over Leyte island close to Tacloban City on 8 November 2013. Tacloban was hit by the full force of the typhoon, which caused  Since Leyte island is a tropical region and is covered by clouds most of the time, we first aimed at using Synthetic Aperture Radar (SAR) data. However, only Landsat images were available and found suitable for our study in the GEE data sets for the time of the disaster and recovery processes. Accordingly, atmospherically corrected surface reflectance Landsat 7 and 8 satellite images with 30 m spatial resolution were used to map the selected land cover types and implement the developed approach for post-disaster recovery assessment (Table 1). In total, 446 Landsat images were employed in this study. In particular the potential of big geospatial time-series data (i.e., Landsat 7 and 8 images) available on the cloud was used to generate the best-available-pixel image composites (explained in the methodology section) for the pre-disaster, event time, and post-disaster time steps over Leyte island to minimize/remove the cloud and shadow areas.  Since Leyte island is a tropical region and is covered by clouds most of the time, we first aimed at using Synthetic Aperture Radar (SAR) data. However, only Landsat images were available and found suitable for our study in the GEE data sets for the time of the disaster and recovery processes. Accordingly, atmospherically corrected surface reflectance Landsat 7 and 8 satellite images with 30 m spatial resolution were used to map the selected land cover types and implement the developed approach for post-disaster recovery assessment (Table 1). In total, 446 Landsat images were employed in this study. In particular the potential of big geospatial time-series data (i.e., Landsat 7 and 8 images) available on the cloud was used to generate the best-available-pixel image composites (explained in the methodology section) for the pre-disaster, event time, and post-disaster time steps over Leyte island to minimize/remove the cloud and shadow areas.

Reference Data and Land Cover Classes
We defined the land cover classes based on the local knowledge of the authors about the study area, and by considering the capacity/potential of the Landsat images (i.e., spatial and spectral resolutions) Appl. Sci. 2020, 10, 4574 5 of 15 to identify different classes through field verification and visual interpretation of high-resolution satellite imagery acquired from various platforms (i.e., WorldView1-3, GeoEye 1, Pleiades) and images available from Google Earth Pro. We identified five land cover classes: forest/trees, built-up, crop land, water body, and other (i.e., non-tree vegetation, bare land, and debris/rubble) ( Table 2), and the corresponding regions of interests (ROIs) were collected for training and testing of the classification method using the tools available within GEE. Table 2. Description of the selected land cover classes.

ID
Land Cover Class Description 1 Forest/Trees All types of trees (e.g., forest, palm and banana), tree canopy coverage >30% 2 Built-up Any type of developed lands such as buildings, roads, impervious surfaces 3 Crop land Agricultural fields with any non-tree crop type plantation 4 Water body Bodies of water including lakes, oceans, rivers and flooded areas 5 Other All other land cover classes (i.e., non-tree vegetation, bare land, rubble and debris)

Methods
The developed approach for post-disaster recovery assessment has three main steps ( Figure 2): (i) land cover classification of the pre-disaster, event time, and post-disaster images, (ii) change detection of pre-disaster and event time classified images to obtain the damage map, and (iii) change analysis of the post-disaster classified images and the damage map to obtain the recovery maps at T2-T4. Generating the land cover maps from the Landsat 7 and 8 satellite images with GEE consisted of three main stages: image composite generation, image classification, and accuracy assessment.
We defined the land cover classes based on the local knowledge of the authors about the study area, and by considering the capacity/potential of the Landsat images (i.e., spatial and spectral resolutions) to identify different classes through field verification and visual interpretation of highresolution satellite imagery acquired from various platforms (i.e., WorldView1-3, GeoEye 1, Pleiades) and images available from Google Earth Pro. We identified five land cover classes: forest/trees, builtup, crop land, water body, and other (i.e., non-tree vegetation, bare land, and debris/rubble) ( Table  2), and the corresponding regions of interests (ROIs) were collected for training and testing of the classification method using the tools available within GEE.

Methods
The developed approach for post-disaster recovery assessment has three main steps ( Figure 2): (i) land cover classification of the pre-disaster, event time, and post-disaster images, (ii) change detection of pre-disaster and event time classified images to obtain the damage map, and (iii) change analysis of the post-disaster classified images and the damage map to obtain the recovery maps at T2-T4. Generating the land cover maps from the Landsat 7 and 8 satellite images with GEE consisted of three main stages: image composite generation, image classification, and accuracy assessment.

Image Composite Generation
We generated the best-available-pixel image composites for the T0-T4 time steps by extracting the best observations from several Landsat 7 and 8 images for our case study (Leyte island), using pixel-based image compositing within the GEE platform, and implemented using rule-based criteria such as dates of acquisition and exclusion of clouds and shadows [45][46][47][48]. To do so, we adapted a script within the GEE environment provided by De Alban et al. [43]. The generated Landsat image composites were exported for subsequent use in the image classification. The developed script needs user-defined inputs and parameters: the coordinates and extent of Leyte island, the number of years, starting and ending date range for the selected region (i.e., Leyte), the thresholds for cloud detection and masking (10% or less) and shadows (z-score = −1), and the image collection used, (i.e., Landsat 7 ETM+ and Landsat 8 OLI). Then the areas that were masked out in the previous step were filled by the median of the selected image pixels [43,49]. The scripts developed for this study are provided as supplementary materials in this paper.
In addition, we calculated five well-known indices from the Landsat composite images and added to the available Landsat bands to be used in the classification step as follows: Normalized Difference Vegetation Index (NDVI; [50]) using Equation (1): The Soil-Adjusted Total Vegetation Index (SATVI; [51]) using Equation (2): Enhanced Vegetation Index (EVI; [52,53]) using Equation (3): Land Surface Water Index (LSWI; [54,55]) using Equation (4): Normalized Difference Tillage Index (NDTI; [56]) using Equation (5): where Blue is the blue band, Red is the red band, NIR is the near-infrared band, SWIR1 is the shortwave infrared 1 band, and SWIR2 is the shortwave infrared 2 band of the Landsat image. NDVI was shown to be a useful tool for distinguishing green vegetation areas from other land cover types; LSWI is for detecting water bodies and when combined with NDVI improves the performance of separating crop lands and forests [57]. SATVI, EVI, and NDTI have been shown to be effective tools for identifying forest and crop types from other land cover classes [43,58].

Image Classification
The final stacked/combined images for each time step consisted of seven Landsat bands (i.e., Blue, Green, Red, NIR, SWIR1, SWIR2, and TIR) and the computed indices (i.e., NDVI, EVI, SATVI, NDTI, and LSWI). The combined images were stored under the "Asset" tab of the GEE code editor platform to be later used for training and testing ROI selection/extraction.
We collected regions of interest (ROIs)/training areas from the images, using the tools within GEE in a way that each class had at least 40 polygons. These selections were based on field verifications and visual interpretations from high-resolution Google Earth Pro historical images and high-resolution satellite images acquired by different sensors (i.e., WorldView1-3, GeoEye1, Pleiades), and Landsat true-and false-color images to obtain final robust training and testing datasets. From the collected ROIs, 70% were randomly selected for training the classifier method and the rest were used for testing purposes. We used Random Forest, as an ensemble machine learning algorithm, to classify the images based on the collected training samples. Random Forest has been widely used in RS image processing for different applications [59], including land cover and land use classification [60,61], landslide detection [62], and hyperspectral image classification [63], due to providing efficient and accurate results. It is also one of the machine learning-based methods available and ready to use in GEE. This approach has also been used by other researchers and its ability to produce accurate classification results has been demonstrated [64][65][66][67]. In addition, Random Forest was selected as the best classifier among several advanced algorithms, including SVM and neural networks, for classifying hundreds of different real-world datasets [68].
Random Forest grows multiple decision trees on randomly selected training subsets. In addition, it takes the advantage of the bootstrapping approach to construct the decision trees from training samples and input variables at every node. The constructed numerous decision tress then go through the majority voting mechanism to obtain the final classification result. This strategy contributes to overcoming the limitations of single tree-based classifiers such as overfitting, and makes it less sensitive to noise. In the current study we executed Random Forest with 100 decision trees per class, a 0.5 fraction number to bag per tree, and 10 for the minimum size of a terminal node.

Accuracy Assessment
The overall user's and producer's accuracies were calculated for each time step classified maps (i.e., T0-T4) to assess the accuracy of the produced land cover classification results. A proportion of 30% of the collected ROIs for each class were selected using stratified random sampling for testing and computing the accuracy measures [51]. To do so, we adapted the code script used by De Alban et al. [43] and executed it within the GEE environment. Selecting the training and ground truth ROIs was challenging for some classes due to similarities between the land cover classes (e.g., crop land and vegetation classes). We tried to minimize the effect of these inaccuracies using different sources (e.g., very high-resolution satellite images) in selecting the ROIs. In addition, considering debris/rubble as the damaged area may cause inaccuracies, e.g., in the case of storm surge that may relocate/wash debris to other regions covering the intact built-up area [13].

Change Analysis
Final damage and recovery maps were generated based on the amount of change in the forest, built-up, and crop land classes. The changes were computed at the municipality level, and the damage map was generated by comparing the T0 and T1 land cover classification results; also, the recovery maps for each post-disaster time step (T2-T4) were computed by comparing the T1 and T2-T4 land cover classification results.

Results and Discussion
The proposed approach produced overall accuracies of 92.8%, 88.6%, 93.8%, 94.6%, and 89.0% for land cover classifications of the T0, T1, T2, T3, and T4 images, respectively (Table 3). In addition, Figure 3 shows the original true-color image composites and the results of the land cover classification for the T0-T4 time steps. Visual interpretation of the results and the high accuracy rates demonstrate the robustness of the approach and the Random Forest classifier for land cover classification for the defined classes in such a case study. However, the producer's accuracy values of the built-up and the other (including non-tree vegetation, bare land, and debris) classes for T1 and T4 images were 76.5% and 75.4%, respectively, which shows the relatively high omission errors for those classes. Relocation of palm tree fronds by the typhoon made distinguishing trees, crop land, and other (including debris) classes more difficult for the event time, and resulted in the lowest overall accuracy of 88.6% for the T1 image classification. Furthermore, the generated composite images may still include cloud coverage. For example, the presence of clouds in the T4 image even after image composition generation resulted in some inaccuracies that produced the second-lowest overall accuracy of 89.0% among other images.
The developed method produced accurate results in classifying built-up, water body, and forest/trees classes in most of the images. Table 3. The land cover classification accuracies for T0, T1, T2, T3, and T4 time epochs for Leyte island. PA-producer's accuracy; UA-user's accuracy; OA-overall accuracy; Non-tv-Non-tree vegetation.  Figure 4 shows the land cover class change trends extracted from the T0 to T4 images. Typhoon Haiyan clearly had a significant impact on the forest/trees and built-up areas, reducing their coverage by about 40% and 30% from T0 to T1 images, respectively. The increase of the other class from T0 to T1 is mostly due to the debris and rubble caused by Haiyan and denuded trees, which led to an increase in non-tree vegetation coverage even in (formerly) forested areas. The storm surge during the typhoon resulted in the increase in water bodies in the T1 image. Monitoring three years after Typhoon Haiyan for Leyte island demonstrates that most of the land cover changes returned to their pre-disaster situation. However, there is an increase in the built-up class that is due to the construction of resettlement sites (as safe zones) mostly around the eastern coastal cities (e.g., Tacloban city- Figure 3). Figure 5 shows the damage and recovery maps for the post-Typhoon Haiyan situation generated at the municipality level for Leyte island. The damage map was produced based on the change analysis of the land cover classes for the T0 and T1 times, and recovery maps (RC) (i.e., RC1, RC2, and RC3) were produced based on change analysis of the damage map and the land cover classes of T2, T3, and T4 times, respectively.  T1 is mostly due to the debris and rubble caused by Haiyan and denuded trees, which led to an increase in non-tree vegetation coverage even in (formerly) forested areas. The storm surge during the typhoon resulted in the increase in water bodies in the T1 image. Monitoring three years after Typhoon Haiyan for Leyte island demonstrates that most of the land cover changes returned to their pre-disaster situation. However, there is an increase in the built-up class that is due to the construction of resettlement sites (as safe zones) mostly around the eastern coastal cities (e.g., Tacloban city- Figure 3).  Figure 5 shows the damage and recovery maps for the post-Typhoon Haiyan situation generated at the municipality level for Leyte island. The damage map was produced based on the change analysis of the land cover classes for the T0 and T1 times, and recovery maps (RC) (i.e., RC1, RC2, and RC3) were produced based on change analysis of the damage map and the land cover classes of T2, T3, and T4 times, respectively.  pre-disaster situation. However, there is an increase in the built-up class that is due to the construction of resettlement sites (as safe zones) mostly around the eastern coastal cities (e.g., Tacloban city- Figure 3).  Figure 5 shows the damage and recovery maps for the post-Typhoon Haiyan situation generated at the municipality level for Leyte island. The damage map was produced based on the change analysis of the land cover classes for the T0 and T1 times, and recovery maps (RC) (i.e., RC1, RC2, and RC3) were produced based on change analysis of the damage map and the land cover classes of T2, T3, and T4 times, respectively.  In Figure 5, the darker red color shows more damage while, for the recovery maps, the darker green color shows a higher level of recovery. The most damaged/impacted municipalities are located in the center and northern part of the island. However, since Haiyan passed through the northern part of the island, the municipalities in the southern parts have the lowest damage ratios, close to zero. The most affected land cover classes were the forest/trees and built-up classes ( Figure 4). Consequently, the municipalities with more forest/trees and built-up coverage compared to other land cover types have higher damage ratios. The change in the recovery maps (from RC1 to RC3) shows that the municipalities with more forest/tree coverage recovered more slowly than the others, and eventually, after a period of three years, most of the forest trees had recovered (RC3) (Figure 5d). Moreover, the results reported in a previous study [4] show almost the same pattern of recovery in terms of change in land cover classes at a small scale (i.e., Tacloban) in the post-Haiyan time. The only difference is in built-up class, where our study shows growth in built-up areas three years after Haiyan, while they reported that Tacloban city returned to almost the same level of built-up coverage four years after the disaster. This is mainly because their focus was on central urban areas of Tacloban, not including the northern part of Tacloban in their analysis, in which the primary resettlement sites for the period after Haiyan were built (Figure 3).

Discussion and Conclusions
The aim of this study was to assess the suitability of GEE for post-disaster recovery assessment for a large-scale region after Typhoon Haiyan, which hit Leyte island in the Philippines in 2013. We proposed an approach to monitor and quantify the recovery process using GEE, in particular with Landsat images, based on land cover classification and subsequent change analysis, and then producing municipality-based damage and recovery maps. The land cover classification results show the robustness of the employed approach, including image composite generation and Random Forest-based image classification. Furthermore, the change analysis and generated recovery maps clearly illustrate the damage and recovery processes during the selected time epochs, i.e., three time-steps for recovery assessment for three years after Haiyan. However, among the big geo dataset available in the platform and considering the time of Typhoon Haiyan (i.e., November 2013), we only found Landsat images suitable for this study. Hence, due to the limited spatial resolution of the data, detailed recovery analysis was not possible (e.g., extracting more land-use classes, such as building use and crop types). As also concluded by Sheykhmousa et al. [4], land cover maps can provide the main information for decision-makers to have an overview of the damage and reconstruction/recovery processes, which are critical in the early stages of recovery. However, a comprehensive recovery assessment needs more details that can be obtained using land use maps extracted from very high-resolution images. Little work has been done to date in damage and recovery mapping of the selected area, and none of this previous research investigated the entire Leyte island. In addition, the available studies that focused on some selected small areas were based on very high-resolution images with different class definitions. Hence, the validation is only based on the accuracy assessment of the land cover classification results. However, since GEE also allows to upload user images for processing, it might be used to upload a very high-resolution image for validation purposes.
Moreover, the change in types of crops and trees, or any detail about the recovered trees, can provide more information about the recovery process. For example, it is essential from the environmental and ecological points of view to understand whether the same trees that were denuded by the typhoon recovered, or whether the trees detected were newly planted after the disaster. More detailed information on the age of the trees can provide significant additional insights into the post-disaster recovery process.
In conclusion, further studies are needed to investigate the potential of the more recently available higher resolution optical and SAR images (Sentinel 1 and 2 images) for most of the case studies (i.e., disasters) that occurred after those platforms became operational. In addition, SAR data have been shown to be a useful tool for extracting land cover and use classes such as distinguishing different crop and tree types [43]. Hence, using higher resolution images and SAR data available in GEE will allow the researchers to detect damage to the crops, contribute to detecting changes in tree types, and extract more details for urban areas after a disaster [69][70][71].
Computer vision and image processing/classification methods deal with the actual image-based pixel values and try to classify them to achieve the best accuracies. However, in post-disaster remote sensing images, especially in water-related disasters, seeing actual damage and visually distinguishing damaged and intact areas is nearly impossible in some cases, as some intact areas are often buried under debris. Hence, assuming the debris area as the damaged class is not entirely correct, there is a need for further studies to identify the debris types (e.g., [13]) or use proxies for distinguishing between damaged and other classes [5] instead of only using direct image classification methods.
Evaluating the recovery process with more details (e.g., building use detection, crop type classification) and several time epochs at such a large scale, and providing municipality level recovery scores similar to those proposed in our study, can be further used for per-municipality resilience assessment at large scales, assuming that the speed of recovery is a proxy for evaluating resilience. However, it has been shown that this assumption is not correct for all cases, and further investigation is needed to clearly determine the relationship between the speed of the post-disaster recovery and resilience [6].
Although the Random Forest-based image classification produced a high accuracy rate, using more advanced methods (i.e., deep learning-based approaches) within the Google Earth API could help to improve the accuracy and would allow executing more advanced image classification and object detection tasks.