Wildﬁre Damage Assessment over Australia Using Sentinel-2 Imagery and MODIS Land Cover Product within the Google Earth Engine Cloud Platform

: Wildﬁres are major natural disasters negatively affecting human safety, natural ecosystems, and wildlife. Timely and accurate estimation of wildﬁre burn areas is particularly important for post-ﬁre management and decision making. In this regard, Remote Sensing (RS) images are great resources due to their wide coverage, high spatial and temporal resolution, and low cost. In this study, Australian areas affected by wildﬁre were estimated using Sentinel-2 imagery and Moderate Resolution Imaging Spectroradiometer (MODIS) products within the Google Earth Engine (GEE) cloud computing platform. To this end, a framework based on change analysis was implemented in two main phases: (1) producing the binary map of burned areas (i.e., burned vs. unburned); (2) estimating burned areas of different Land Use/Land Cover (LULC) types. The ﬁrst phase was implemented in ﬁve main steps: (i) preprocessing, (ii) spectral and spatial feature extraction for pre-ﬁre and post-ﬁre analyses; (iii) prediction of burned areas based on a change detection by differencing the pre-ﬁre and post-ﬁre datasets; (iv) feature selection; and (v) binary mapping of burned areas based on the selected features by the classiﬁers. The second phase was deﬁning the types of LULC classes over the burned areas using the global MODIS land cover product (MCD12Q1). Based on the test datasets, the proposed framework showed high potential in detecting burned areas with an overall accuracy (OA) and kappa coefﬁcient (KC) of 91.02% and 0.82, respectively. It was also observed that the greatest burned area among different LULC classes was related to evergreen needle leaf forests with burning rate of over 25 (%). Finally, the results of this study were in good agreement with the Landsat burned products.


Introduction
Natural disasters are generally divided into two main categories: geophysical (e.g., wildfires, floods, earthquakes, and droughts) and biological (e.g., biotic stresses) [1]. Among natural geophysical disasters, wildfire occurs more than a thousand times annually over different regions of the globe [2][3][4] and can cause considerable economic and ecological damages, as well as human casualties. Additionally, wildfires have significant longterm impacts on both the environment and people, such as changing the structure of the ecosystems, soil erosion, destructing wildlife habitats, and increasing potential of flooding [4]. Moreover, wildfires emit a considerable amount of greenhouse gases (e.g., carbon dioxide and methane) which results in global warming [5]. Therefore, reliable, timely, and detailed information on burned areas in a wildfire event is necessary.
Wildfire mapping using traditional methods (e.g., field surveys and manual digitization) has many limitations. Although traditional methods can provide the highest accuracies, they are time-consuming and have limitations in terms of automation and repeating measurements over a long period of time. However, automated RS applications 2 of 30 are more efficient in terms of cost, time, and coverage [6][7][8], and facilitate accurate wildfire mapping with many advanced machine learning algorithms (e.g., deep learning).
RS has been widely used for wildfire mapping and monitoring. For instance, Dragozi, et al. [9] investigated the effects of derivative spatial and spectral features to map burned areas using IKONOS imagery in Greece. To this end, they utilized a Support Vector Machine (SVM) and fuzzy complementary criterion for classification and feature selection, respectively. Moreover, Oliva and Schroeder [10] evaluated the Visible-Infrared Imaging Radiometer Suite (VIIRS) active fire products with the resolution of 375 (m) to estimate burned areas in ten different areas worldwide. Based on this research, the burned areas depended on the environment features and behavior of fire. Additionally, Chen, et al. [11] detected forest burned areas based on multiple methods using Landsat TM data. Their proposed method utilized four spectral indices to determine burned areas. Afterwards, the damaged areas were obtained by selecting an optimum threshold based on histogram analysis. Furthermore, Hawbaker, et al. [12] proposed a burned area estimation algorithm based on dense time-series of Landsat data and generated the Landsat Burned Area Essential Climate Variable (BAECV) products. This method utilized surface reflectance and multiple spectral indices as inputs into the gradient boosted regression models to generate burned probability maps. Finally, the burned areas were extracted from the generated burned probability maps using a combination of pixel-level thresholding and region growing process. Pereira, et al. [13] also proposed an automatic sample generation method for burned area mapping using VIIRS active fire products over the Brazilian Cerrado savanna. They also used one-class SVM algorithm for burned area mapping. Moreover, Roteta, et al. [14] developed a locally adapted multitemporal burned area algorithm based on Sentinel-2 and Moderate Resolution Imaging Spectroradiometer (MODIS) datasets over the sub-Saharan Africa. This method had two main parts: (1) burned areas were detected based on a fixed threshold, and (2) tile dependent statistical thresholds were estimated for each predictive variable in a two-phase strategy by overlaying the MCD14ML products. Furthermore, Ba et al. [15] developed a burned area estimation method using a single MODIS imagery based on back-propagation neural network and spectral indices. Their method had three main phases: (1) extraction of samples for five classes (cloud, cloud shadow, vegetation, burned area and bare soil), (2) feature selection, and (3) classification based on back-propagation neural network. They also evaluated the performance of the proposed method at three different areas in Idaho, Nevada, and Oregon. Woźniak and Aleksandrowicz [16] also developed an automatic burned area mapping framework based on medium resolution optical Landsat images. This method was implemented in four main steps: (1) extraction of the spectral indices and differencing them, (2) multiresolution segmentation and masking, (3) detection of the core burned areas based on an automatically adjusted threshold, and (4) region growing procedure and neighborhood analysis for detection of final burned areas. Additionally, Otón, et al. [17] presented a global product of burned areas using Advanced Very High Resolution Radiometer (AVHRR) imagery. A synthetic burned area index was proposed using a combination of spectral indices and surface reflectance data in the form of a single variable. The final burned area product was obtained based on the time series classification of synthetic burned area indices using the Random Forest (RF) algorithm. Finally, Liu, et al. [18] evaluated the performance of ICESat-2 photon-counting LiDAR data in estimating burned and unburned areas using both RF and logistic regression methods. They tested recent fires in 2018 in northern California and western New Mexico. This study confirmed the feasibility of employing ICESat−2 data for burned forest classification.
RS application for wildfire mapping and change analysis over large areas has several challenges, such as efficiently processing big RS data (e.g., thousands of satellite images). To address this issue, several cloud computing platforms have been so far developed, one of the most widely used among which is Google Earth Engine (GEE) [19,20]. This platform has brought a unique opportunity for undertaking Earth Observation (EO) research studies. GEE has been designed for parallel processing, storing, and mapping different types of RS Remote Sens. 2021, 13, 220 3 of 30 datasets at various spatial scales. GEE has been extensively used for various applications, including wildfire mapping and trend analysis due to its key advanced characteristics, such as being free to use, providing high-speed parallel processing without downloading data, and being user-friendly. For example, Long, et al. [21] developed an automated framework to map burned areas at a global-scale using dense time-series of Landsat-8 images within GEE for 2014 to 2015. This method was based on feature generation from various spectral indices with a RF classifier. Moreover, Zhang, et al. [22] generated a global burned area product using Landsat-5 satellite imagery and a RF classifier within GEE. They also analyzed the spatial distribution pattern of the global burned areas. Furthermore, Barboza Castillo et al. [23] investigated forest burned areas using a combination of Sentinel-2 and Landsat-8 satellite images within the GEE cloud platform over Uttarakhand, Western Himalaya. They mapped forested burned area by differencing the spectral indices for post and pre fire events as well as combining unsupervised clustering and supervised classification methods.
As discussed, although many methods have been so far proposed to map burned areas, a few of them were developed within GEE. Moreover, although different algorithms have been developed for mapping burned areas using RS datasets, there are still several limitations as follows: (1) Many burned area products contain moderate spatial resolution. However, with the increasing availability of higher spatial resolution satellite imagery (e.g., Sentinel-2 with 10 (m) spatial resolution), there is potential to produce more detailed products in terms of spatial resolution. (2) Many RS methods for burned area mapping which are based on high-resolution imagery are complex and do not support the mapping of burned areas over large regions. (3) Although multi sensor-based methods for burned area mapping provide promising results, most of them are not computationally efficient. (4) A thresholding method which is applied in many research studies for discriminating burned from unburned areas does not provide accurate results for large regions, because burned areas at different regions depend on the characteristics of ecosystem and behavior of fire. Thus, a dynamic thresholding method should be developed to obtain high accuracies over different regions. (5) Many methods are based on the binary mapping (i.e., burned vs. unburned). However, estimation of the Land Use/Land Cover (LULC) over burned areas is necessary for many applications. (6) Many methods only use some specific spectral features. However, the potential of spatial features should be also comprehensively investigated to improve wildfire mapping and monitoring.
To address these issues, a framework was proposed in this study by leveraging GEE big data processing platform. The proposed method was implemented in two phases. The first phase detected burned area and was implemented in five main steps: (1) preprocessing, (2) spatial and spectral feature extraction; (3) change detection by image differencing on spatial and spectral features; (4) feature selection using the Harris's Hawk Optimization (HHO) algorithm; and (5) mapping burned area by applying a supervised classifier to the selected features. The second phase estimated the type of burned LULC using MODIS LULC products.
The main contributions of this study are: (1) mapping Australian burned areas and damage assessment at a higher spatial resolution compared to other studies in GEE; (2) estimating burned area for different types of LULC using an innovative idea (most previous methods were only focused on binary burned mapping); (3) implementing a novel feature selection algorithm (i.e., HHO) for determining the most optimal spectral and spatial features for burned area mapping, (4) investigating the importance of features for burned area mapping; (5) utilizing a combination of spatial and spectral features for burned area mapping; and (6) comparing the performance of different classifiers for burned area mapping.

Study Area
The study area is mainland Australia with an area of 7.692 (Mkm 2 ). This study uses Australian wildfire as a case study to evaluate the accuracy of the proposed method to map the burned area between September 2019 to February 2020. Australian wildfire started in the southeast of the country in June 2019 and rapidly spread throughout the continent [24]. Based on several reports, at least 28 people were killed, and more than 3000 homes were destroyed or damaged during this incident. Moreover, during the Australian wildfires an enormous amount of CO 2 was released into the air which negatively impacted the quality of air. Due to the large extent of the Australian continent, it has six major climatic zones (i.e., desert, temperate, equatorial, subtropical, tropical, and grassland climates). Most areas experience the annual average rainfall of less than 600 (mm), and in some areas the annual average rainfall can be more than 1200-8000 (mm) [25].

Satellite Data
In this study, Sentinel-2 optical satellite images were utilized. Sentinel-2 is an Earth Observation (EO) mission launched by the European Space Agency (ESA) to continue ESA's global services on multispectral high spatial resolution observations. This mission includes the Sentinel-2A and Sentinel-2B satellites. The temporal resolution of the constellation is 5 days. The MultiSpectral Instrument (MSI) is the main sensor of Sentinel-2 and is based on the pushbroom concept. MSI provides 13 spectral bands with a wide spectral coverage over the visible, Near Infrared (NIR) and Shortwave Infrared (SWIR) domains at different spatial resolutions from 10 (m) to 60 (m) [26]. The data is freely available at Sentinel Scientific Data Hub (https://scihub.copernicus.eu). In this study, the level 1C products of Sentinel-2, which are available in GEE (Dataset ID: ee.ImageCollection("COPERNICUS/S2_SR"), were used. All maps are presented in a reference geographic coordinate system of the World Geodetic System 1984 (WGS 1984).

Reference Data
Reference samples are necessary to train any supervised classification algorithm. In this study, reference samples for the burned areas were generated using multiple reports about the locations of wildfires (https://www.ktnv.com/) [27]. Figure 1a illustrates the location of burned areas generated from these reports. Furthermore, the samples for the burned and unburned areas were generated by visual interpretation of the Sentinel-2 time series imagery, acquired during September 2019 to February 2020. In addition, Figure 1a was used for generation of sample data for several burned areas. Figure 1b illustrates the spatial distribution of the reference samples for both burned and unburned areas and some regions for visual accuracy assessment, and Table 1 provides the number of samples. Finally, all samples were randomly divided into three groups: training (50%), validation (17%), and test (33%). The validation dataset was used as an intermediate dataset for the fine tuning of the utilized classifiers [28].

MODIS LULC Product
The MODIS LULC product was used in this study to determine the LULC types over the burned areas. Figure 2 demonstrates the MODIS LULC map over Australia (https://modis.gsfc.nasa.gov/). This map contains 17 classes (see Table 2 and Figure 2). The original spatial resolution of the MODIS LULC map is 500 (m) which was resampled to 10 (m) to conform the spatial resolution of Sentinel-2. As is clear, the dominant classes are open Shrublands and Grasslands. Forests are also mainly found in south east and south west of the country. In this study, the MODIS LULC products available in GEE (Dataset ID: ee.ImageCollection("MODIS/006/MCD12Q1")) were used.

MODIS LULC Product
The MODIS LULC product was used in this study to determine the LULC types over the burned areas. Figure 2 demonstrates the MODIS LULC map over Australia (https://modis.gsfc.nasa.gov/). This map contains 17 classes (see Table 2 and Figure 2). The original spatial resolution of the MODIS LULC map is 500 (m) which was resampled to 10 (m) to conform the spatial resolution of Sentinel-2. As is clear, the dominant classes are open Shrublands and Grasslands. Forests are also mainly found in south east and south west of the country. In this study, the MODIS LULC products available in GEE (Dataset ID: ee.ImageCollection("MODIS/006/MCD12Q1")) were used.

Landsat Burned Area Product
Burned area product of Landsat was used to assess the results of the proposed method in this study. The Landsat level 3 burned area product is designed to detect burned areas across all ecosystems [30]. This product is available at spatial resolution of 30 (m) from 2013 to present. The Landsat-8 burned product is generated based on surface reflectance and top of atmosphere brightness temperature data. In this study, the burned area products of Landsat-8 generated from September 2019 to February 2020 were used after aggregating them based on a burned probability threshold of 0.8 ( Figure 3).

Proposed Method
The flowchart of the proposed method is provided in Figure 4. The proposed framework was implemented in two main phases: (1) binary burned area mapping (i.e., burned vs. unburned), (2) estimation of burned areas across different LULC types. More details on the proposed framework are provided in the following subsections.

Phase 1: Binary Burned Areas Mapping
The first phase of the proposed method for binary burned areas mapping included five steps: (1) preprocessing; (2) spatial-spectral feature extraction; (3) change detection to predict burned areas by differencing pre-fire and post-fire layer-staked images (i.e., original spectral bands, spectral indices, and texture features); (4) feature selection using the HHO optimizer and supervised classifiers; and (5) binary burned area mapping using the selected features and the optimum classifier algorithm. These steps are discussed in more details in the following subsections.

Preprocessing
In this study, all preprocessing steps were performed in the GEE big data processing platform. The atmospheric correction was performed using the module Py6S that is available at https://github.com/robintw/Py6S. Moreover, the cloudy regions in the images were masked using the F-Mask algorithm with a threshold of 10%. Finally, all the spectral bands were resampled to 10 (m).

Feature Extraction
Although most wildfire studies conducted using RS methods have only utilized spectral features, spatial features could also provide valuable information to improve the results. Therefore, both spectral and spatial features were extracted and investigated in this study.
The spectral indices have many applications in RS image analysis [31]. The main characteristic of these indices is that they make some features more discernible compared to the original data (i.e., spectral bands). In this study, 90 different spectral indices were used (see Appendix A). These features could improve the result of burned area mapping. Since the proposed method was based on change detection, several areas were unwantedly included in changes, such as the boundary of water bodies and rapidly growing vegetation. Therefore, several water indices were used to reduce false alarms originated from undesired changes in water bodies. Vegetation indices were also necessary due to diversity of vegetation over the entire study area and the similarity of the burned vegetation response with those of croplands. Many research studies have argued that the burned indices improved the accuracy of burned area mapping [15,32,33].

Landsat Burned Area Product
Burned area product of Landsat was used to assess the results of the proposed method in this study. The Landsat level 3 burned area product is designed to detect burned areas across all ecosystems [30]. This product is available at spatial resolution of 30 (m) from 2013 to present. The Landsat-8 burned product is generated based on surface reflectance and top of atmosphere brightness temperature data. In this study, the burned area products of Landsat-8 generated from September 2019 to February 2020 were used after aggregating them based on a burned probability threshold of 0.8 ( Figure 3).

Proposed Method
The flowchart of the proposed method is provided in Figure 4. The proposed framework was implemented in two main phases: (1) binary burned area mapping (i.e., burned The spatial features, such as texture indices, can also improve the results of LULC classification [34]. Texture refers to the relationship between an individual pixel with its neighboring pixels and provides valuable information for LULC classification. Different features, such as density, equality, non-roughness, and size uniformity, can be generated by texture analysis. In this study, 17 texture features were generated from the Grey Level Cooccurrence Matrices (GLCM) (seen Appendix B). It is worth noting that textural features are usually extracted from a panchromatic band. However, since Sentinel-2 does not contain a panchromatic band, three visible bands were combined based on Equation (1) to generate the panchromatic band.
where B pan is the panchromatic band and B 2 , B 3 , and B 4 are the blue, green, and red bands, respectively. The window size for texture analysis was set to 3 × 3.  The first phase of the proposed method for binary burned areas mapping included five steps: (1) preprocessing; (2) spatial-spectral feature extraction; (3) change detection to predict burned areas by differencing pre-fire and post-fire layer-staked images (i.e., original spectral bands, spectral indices, and texture features); (4) feature selection using the HHO optimizer and supervised classifiers; and (5) binary burned area mapping using the selected features and the optimum classifier algorithm. These steps are discussed in more details in the following subsections.

Change Detection
Change detection is one of the most important applications of RS and has been widely used for estimating burned areas using pre-and post-fire datasets [32]. In this study, image differencing was used as the change detection algorithm. Image differencing is a simple and popular change detection algorithm [6]. This algorithm is based on the band-to-band pixel subtraction of datasets from the first-and second-time datasets (Equation (2)): where, X b r,c and Y b r,c are the pixel values in the first and second images in the row r and column c for band b, respectively.

Feature Selection
The main purpose of feature selection is eliminating redundant features and using only the most optimal features during the classification [35]. In this study, feature selection was implemented in two steps: (1) defining the number of features, and (2) identifying the most useful features by an optimizer. To this end, the HHO algorithm was employed in this study. HHO is a new population-based and nature-inspired optimization algorithm, which has been widely used for optimization purposes in different fields [36]. HHO is based on the chasing-escaping patterns observed between the hawks and their prey (such as rabbit). In this method, the Harris' hawks are considered as candidate solutions and the intended prey is considered as the best candidate solution in every iteration. The HHO algorithm has three main steps, including amaze pounce (exploration), transforming from exploration to exploitation, and exploitation [37]. The main purpose of exploration phase is mathematically modeling the waiting, searching, and discovering the desired hunt. Transforming from exploration to exploitation is applied based on the external energy of a rabbit. The exploitation is considered the residual energy of the prey. In the HHO algorithm, three parameters need to be initialized: total population size (N), maximum iteration, and number of features. Figure 5 illustrates the proposed feature selection framework using the HHO algorithm. HHO is initialized with the number of agents, number of features and number of iterations. The initial values of the three HHO parameters were manually set after trial and error. Finally, 15, 500, and 60 were selected as the optimum values for the total population size, maximum iteration, and number of features, respectively. Then, the population of agents (solutions) was evaluated by a fitness function (i.e., Overall Accuracy (OA) of the classifier). The fitness of population of agents was calculated based on the performance of the classifiers on the test dataset. In this study, the supervised classifiers of the k-Nearest Neighbors (kNN), Support Vector Machine (SVM), and Random Forest (RF) were trained using the training samples and their performances were evaluated using the validation dataset. Since classifiers are sensitive to input features, it was required to tune the classifier parameters when the input features were changed [20,38]. In this study, the optimization of the tuning parameters of the classifiers was applied based on the grid search algorithm. The feature selection process was continued until all interactions were finished.

Classification
After selecting the most optimal features using the HHO algorithm, they were used as input for three classification algorithms (RF, kNN, and SVM) to produce burned areas maps. Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 31

Classification
After selecting the most optimal features using the HHO algorithm, they were used as input for three classification algorithms (RF, kNN, and SVM) to produce burned areas maps.

RF Classifier
RF is an ensemble machine learning algorithm that has been widely used for RS classifications [39][40][41]. This method combines a number of weak classifiers to construct a powerful classifier. RF combines a set of decision trees based on randomly selected subsets of data where each decision tree classifies input data. A bootstrap sample approach is also used to train each tree using a set of sample data. Then, the highest number of votes as was chosen as the classification result [42]. RF classifier has two main tuning parameters which need to be defined by users. These parameters are the number of features to split each node and number of trees to grow into a whole forest. More details about the RF classifier can be found in [42].

k-NN Classifier
The k-NN algorithm is one of the simplest non-parametric classifiers which uses an instance-based learning approach. kNN has been extensively used for image classifications due to the simplicity of the implementation of this algorithm [43]. This algorithm classifies an unlabeled pixel based on the class attributes of its k nearest neighbors. In fact, the kNN algorithm identifies a group of k training samples that have nearest distance to unlabeled pixels. Then, it assigns a label of class to unknown pixels by calculating the average of the response variables. Therefore, k plays an important role in the accuracy of

RF Classifier
RF is an ensemble machine learning algorithm that has been widely used for RS classifications [39][40][41]. This method combines a number of weak classifiers to construct a powerful classifier. RF combines a set of decision trees based on randomly selected subsets of data where each decision tree classifies input data. A bootstrap sample approach is also used to train each tree using a set of sample data. Then, the highest number of votes as was chosen as the classification result [42]. RF classifier has two main tuning parameters which need to be defined by users. These parameters are the number of features to split each node and number of trees to grow into a whole forest. More details about the RF classifier can be found in [42].

k-NN Classifier
The k-NN algorithm is one of the simplest non-parametric classifiers which uses an instance-based learning approach. kNN has been extensively used for image classifications due to the simplicity of the implementation of this algorithm [43]. This algorithm classifies an unlabeled pixel based on the class attributes of its k nearest neighbors. In fact, the kNN algorithm identifies a group of k training samples that have nearest distance to unlabeled pixels. Then, it assigns a label of class to unknown pixels by calculating the average of the response variables. Therefore, k plays an important role in the accuracy of the kNN classifier and should be selected during the identification of the optimized tuning parameters. More details about the kNN classifier can be found in [44].

SVM Classifier
SVM is a supervised non-parametric algorithm based on statistical learning theory that is commonly utilized for classification purposes [45]. The basic idea of the SVM classifier is to build a hyperplane that maximizes the margin among classes [46]. SVM maps the input image data from the image space to feature space, where the classes can be effectively separated by a hyperplane. The mapping can be applied by different types of kernel functions, such as the Radial Basis Function (RBF), which has proved to provide the best performance among kernel types [47]. This classifier has two main tuning parameters: the kernel parameter (e.g., gamma in RBF) and the penalty coefficient (C). One can refer to [44] for more details on this algorithm.
As mentioned above, each of these supervised classifiers has some parameters that should be tuned to obtain the highest possible classification accuracies. One of the simplest methods for tuning these parameters is the grid search. Based on this approach, a range of unknown parameters is initially defined, and the model is built by the training data and is assessed by the validation data. This process continues until the defined bound is finished. Finally, the optimum tuning parameters correspond to the highest accuracy obtained from the validation data.

Accuracy Assessment
In this study, the confusion matrix of classification was used for statistical accuracy assessment. As illustrated in Table 3, a confusion matrix has four components: True Positive (TP), False Positive (FP), True Negative (TN), and False Negative (FN). Various accuracy matrices were generated based on these four parameters and were used to report the accuracy levels. As demonstrated in Table 4, these matrices were Overall Accuracy (OA), Balanced Accuracy (BA), F1-Score (FS), False Alarm (FA), Kappa Coefficient (KC), Precision (PCC), Recall, Miss-Detection (MD), and Specificity [48,49]. Table 3. Confusion matrix [48,49].

Burned Unburned
Actual Burned TP FN Unburned FP TN Table 4. The metrics which were used for accuracy assessment of the burned map in this study [48,49].

Accuracy Index Formula
Overall Accuracy (OA) Besides statistical accuracy assessment, the wildfire map produced by the proposed method was visually evaluated over multiple burned areas with different LULC types.
Moreover, the produced burned area map was compared with the Landsat burned area product described in Section 2.5.

Phase 2: Mapping LULC Types of Burned Areas
After producing a binary burned area map through the phase 1 of the proposed method, the LULC types of the burned areas were identified using the MODIS LULC product available in GEE (see Section 2.4). To this end, the binary burned areas map was overlaid on the MODIS LULC product and then, the amount of burned areas for different LULC classes were calculated.

Parameter Setting and Feature Selection
As discussed in Section 2.6.1 the optimal tuning parameters for the classifiers should be initially selected to obtain the highest possible classification accuracies. Table 5 provides the optimum tuning parameters for the three classifiers (RF, kNN, and SVM) that were obtained using the grid search method. The HHO algorithm has three main parameters that need to be initialized. The number of hawks (population size) is set to 15, number of features for feature selection is set to 60, and the number of iterations is set to 500.
As discussed in Section 2.6.1 121 features (90 spectral features, 13 original spectral bands, 17 spatial features, and one panchromatic band) were initially extracted. Then, the HHO algorithm was applied to select the most optimal features. The convergence curves of the HHO optimizer for the three classifiers were obtained for the test data are illustrated in Figure 6. The kNN and SVM algorithms converged after 250 iterations, and the RF classifier converged after 450 iterations. The feature selection process of the SVM took relatively more time due to solving equations to build a hyperplane at each epoch. Table 6 provides the selected features for each classifier. As clear, the results of the HHO algorithm were different for various classifiers, indicating that each classifier was sensitive to different features. However, multiple common features were identified for all of them (e.g., NBR, DISS, and B2). Remote Sens. 2021, 13, x FOR PEER REVIEW 14 of 31  Table 6 provides the selected features for each classifier. As clear, the results of the HHO algorithm were different for various classifiers, indicating that each classifier was sensitive to different features. However, multiple common features were identified for all of them (e.g., NBR, DISS, and B2).

Results
The confusion matrices of the classifications obtained by different classifiers for two scenarios (using selected features and without feature selection) are illustrated in Figure  7. The statistical accuracy results are also provided in Table 7. It was observed that all classifiers resulted in better accuracies when the HHO feature selection algorithm was employed. Moreover, the processing time was decreased when HHO was used compared to the case when all features were utilized. For example, the OA of the kNN algorithm

Results
The confusion matrices of the classifications obtained by different classifiers for two scenarios (using selected features and without feature selection) are illustrated in Figure 7. The statistical accuracy results are also provided in Table 7. It was observed that all classifiers resulted in better accuracies when the HHO feature selection algorithm was employed. Moreover, the processing time was decreased when HHO was used compared to the case when all features were utilized. For example, the OA of the kNN algorithm was improved from 58.13% to 86.34% when HHO algorithm was applied to the selected features. Moreover, the OA of the SVM algorithm increased from 72.67% to 78.87% and MD rate reduced from 23.26 to 11.08 when feature selection was utilized. Generally, the HHO feature selection method improved the performance of the kNN and SVM classifiers in detection of Burned pixels which was reflected on the TP component in the confusion matrix and led to decreasing FN pixels. The performance of the RF classifier was also considerably improved by feature selection. For example, the detection rate of the unburned pixels (TN) was improved, which reduced the FP pixels. pixels (TN) was improved, which reduced the FP pixels.
Comparing the results of the three classifiers, the kNN algorithm had the poorest performance, especially in detecting burned areas without the feature selection method. The accuracies of kNN for both burned and unburned classes were under 60% and the error of classification was more than 40%. Although the performance of the SVM classifier was improved when feature selection was employed, it had low accuracy in detecting unburned pixels (lower than 70%). Overall, the RF had the highest accuracy and lowest errors in terms of MD and FA. RF had the highest performance in detecting both Burned and Unburned classes. Based on these results, the RF classifier was selected for wildfire mapping over mainland Australia.    Figure 8 shows the binary burned areas based on the RF classifier applied to the selected features. Based on this Figure, the coverage of burned and unburned areas were 249,358 (km 2 ) and 5,162,072 (km 2 ), respectively. As clear from Table 7, the OA and KC of this classification were 91.02% and 0.82, respectively.   Comparing the results of the three classifiers, the kNN algorithm had the poorest performance, especially in detecting burned areas without the feature selection method. The accuracies of kNN for both burned and unburned classes were under 60% and the error of classification was more than 40%. Although the performance of the SVM classifier was improved when feature selection was employed, it had low accuracy in detecting unburned pixels (lower than 70%). Overall, the RF had the highest accuracy and lowest errors in terms of MD and FA. RF had the highest performance in detecting both Burned and Unburned classes. Based on these results, the RF classifier was selected for wildfire mapping over mainland Australia. Figure 8 shows the binary burned areas based on the RF classifier applied to the selected features. Based on this Figure, the coverage of burned and unburned areas were 249,358 (km 2 ) and 5,162,072 (km 2 ), respectively. As clear from  After producing the binary burned area map, phase 2 of the proposed method was applied to detect the LULC types over the burned areas. Figure 9 illustrates the produced map and Table 8 provides the coverage burned areas for different LULC classes. Based on the results, Evergreen Needleleaf Forests have suffered a 25% decrease relative to their own distribution, but the actual area covered by this class is quite small (5629 km 2 ) compared to other LULC classes, such as Grasslands (91,106 km 2 ), Open Shrublands (71,511 km 2 ), and Savannas (27,878 km 2 ). After producing the binary burned area map, phase 2 of the proposed method was applied to detect the LULC types over the burned areas. Figure 9 illustrates the produced map and Table 8 provides the coverage burned areas for different LULC classes. Based on the results, Evergreen Needleleaf Forests have suffered a 25% decrease relative to their own distribution, but the actual area covered by this class is quite small (5629 km 2 ) compared to other LULC classes, such as Grasslands (91,106 km 2 ), Open Shrublands (71,511 km 2 ), and Savannas (27,878 km 2 ). Remote Sens. 2021, 13, x FOR PEER REVIEW 17 of 31

Discussion
In this study, the accuracy of the final wildfire map produced using the proposed method ( Figure 10) was also visually assessed over multiple burned areas identified with blue squares in Figure 1b. These areas were related to different LULC types, such as Evergreen Broadleaf Forests (sample regions 1, 2, and 4 in Figure 1b), Woody Savannas (sample region 3 in Figure 1b), and Grasslands (sample region 5 in Figure 1b). Moreover, the results were compared with those obtained from the Landsat burned area product (Figure 3). Most burned areas were correctly classified by the proposed method (central column in Figure 11), and the results were more accurate than those of the Landsat-8 products (right-hand side in Figure 10). Finally, it should be noted that the produced burned area map in this study had higher spatial resolution compared to the Landsat burned areas product.

Discussion
In this study, the accuracy of the final wildfire map produced using the proposed method ( Figure 10) was also visually assessed over multiple burned areas identified with blue squares in Figure 1b. These areas were related to different LULC types, such as Evergreen Broadleaf Forests (sample regions 1, 2, and 4 in Figure 1b), Woody Savannas (sample region 3 in Figure 1b), and Grasslands (sample region 5 in Figure 1b). Moreover, the results were compared with those obtained from the Landsat burned area product ( Figure  3). Most burned areas were correctly classified by the proposed method (central column in Figure 11), and the results were more accurate than those of the Landsat-8 products (right-hand side in Figure 10). Finally, it should be noted that the produced burned area map in this study had higher spatial resolution compared to the Landsat burned areas product. It is important to select the optimum spectral and spatial features to obtain accurate burned area maps. Based on the numerical results in this study, the feature selection by the HHO algorithm could considerably improve the performance of the classifiers by selecting optimum features for the classification. Moreover, feature selection improved the efficiency of the proposed method in terms of time processing. The sentinel-2 has high potential for burned area detection due to its higher temporal and spatial resolutions compared to other satellites which provide free data. For example, many research studies have been conducted to detect burned areas using Landsat imagery which has low temporal resolution (about 15 days) [3,11,12,[21][22][23]50,51]. The poor temporal resolution of Landsat negatively affects the timely mapping of burned areas. Recently, burned area detection using unnamed aerial vehicle (UAV) imagery were also investigated. However, these methods have their own limitations, especially over large areas such as Australian continent [52][53][54][55][56]. Overall, based on the numerical and visual analysis, the sentinel-2 has a high potential for burned area mapping.
Although selecting optimum features is an important step in burned area mapping, many studies have ignored it and only employed a few spectral features [56,[58][59][60]. Based on the experiments in this study, both spectral and spatial features play key roles to obtain accurate burned maps. For example, based on Figure 11, the spectral indices and spatial features had high importance compared to the original spectral bands. In summary, the spectral indices had the highest importance. For instance, the BAIS2 index had the highest impact on the results by providing the importance of more than 0.2. Additionally, band 9 and 12 of Snetinel-2 had the highest importance compared to other spectral bands. Finally, SAVG and shade were the most important spatial features for burned area mapping. Consequently, it is important to use both spatial and spectral features to produce accurate burned area maps.
It is important to select the optimum spectral and spatial features to obtain accurate burned area maps. Based on the numerical results in this study, the feature selection by the HHO algorithm could considerably improve the performance of the classifiers by selecting optimum features for the classification. Moreover, feature selection improved the efficiency of the proposed method in terms of time processing.
Many studies have proved the high potential of the deep learning methods for burned area detection [61][62][63]. However, these methods need a large amount of training samples, the collection/generation of which is usually challenging. One of the advantages of the proposed method is that it could produce accurate burned are maps using relatively lower numbers of training samples (e.g., 11,000 pixels) compared to the commonly used deep learning algorithms. Furthermore, most of the deep learning methods require significant time for training the network, tuning the hyperparameters, and designing efficient architecture. However, the proposed framework had few tuning parameters and improved the time efficiency by an optimum architecture.

Conclusions
This study estimated burned areas over mainland Australia using Sentinel-2 imagery at spatial resolution of 10 (m) from September 2019 to February 2020 within the GEE platform. This study produced an accurate burned area map by (1) estimating burned area for different types of LULC, (2) utilizing a novel feature selection method for determining the most optimal spectral and spatial features for burned area mapping, (3) utilizing a combination of spatial and spectral features for burned area mapping, and (4) comparing the performance of different classifiers for burned area mapping. The proposed framework was based on image differencing and supervised machine learning algorithms. Many spectral (vegetation indices, water indices, soil indices, and burned area indices) and spatial (texture features extracted from the GLCM matrix) features were initially generated. Then, the performance of three classifiers was investigated through two scenarios: (1) classification based on all spectral and spatial features, (2) classification based on the selected features by the HHO algorithm. It was observed that feature selection significantly improved the accuracy of the classifications. Additionally, the RF classifier had the highest accuracy (OA = 91.02%). It was also observed that the proposed method had a higher accuracy compared to the Landsat burned area product. Based on the burned area map produced by the proposed method, about 250,000 (km 2 ) across mainland Australia was damaged by wildfires. Based on the obtained results by feature selection, the BASI2 index was the most importance feature among other spectral and spatial features, as well as the original bands. Among different LULC classes, the Evergreen Needleleaf went through the most significant damage relative to its extent. Overall, it was concluded that the proposed algorithm within the GEE cloud platform had a high potential for burned area mapping over large areas in a costly, timely, and computationally efficient manner.
However, the main limitation of the proposed approach was using the MODIS LULC product within phase 2 of the method to identify the type of LULC over the burned areas. This product has low spatial resolution (0.5 km). Thus, future studies should utilize a higher resolution LULC product or perform the LULC classification using Sentinel-2 images in the phase 2 of the proposed method.

Acknowledgments:
The authors would like to thank the European Space Agency (ESA) for providing the Sen-tinel-2 Level-1C products and National Aeronautics and Space Administration (NASA) for providing Landsat burned product and MODIS land cover dataset. We thank the anonymous reviewers for their valuable comments on our manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Vegetation Index [83] 30 GBNDVI Green-Blue normalized difference vegetation index Vegetation Index [    Bare Soil Index [124]  Burned Index [102] Appendix B Table A2. Textural features computed from the gray level co-occurrence matrix (GLCM). i + j − µ x − µ y 4 * P(i, j) Asymmetry of image Where the P(i, j) is the (i, j)-th entry of the normalized GLCM, N is the total number of gray levels in the image, and x, y and x, y denote the mean and standard deviation of the row and column sums of the GLCM, respectively. HXY = − ∑ N i=1 ∑ N j=1 P(i, j)log 2 (P(i, j) + ε), HXY1 = − ∑ N i=1 ∑ N j=1 P(i, j) log 2 P x (i)P y (j) + ε , HX = − ∑ N j=1 P x (i) log 2 (P x (i) + ε), HY = − ∑ N j=1 P y (i) log 2 P y (j) + ε , and HXY2 = − ∑ N i=1 ∑ N j=1 P x (i)P y (j) log 2 P x (i)P y (j) + ε .