ASHFALL DISPERSAL MAPPING OF THE 2020 TAAL VOLCANO ERUPTION USING DIWATA-2 IMAGERY FOR DISASTER ASSESSMENT

Natural disasters incur many fatalities and economic losses for vulnerable and developing countries such as the Philippines. It is crucial that during calamities, on-ground surveillance is supplemented by low-cost and time-efficient methods such as satellite remote sensing. Diwata-2 is a Philippine microsatellite specifically equipped for disaster assessment. In this study, the capabilities of this satellite in ashfall detection were explored by closely examining the case of the Taal volcano eruption on January 12, 2020. Satellite images covering parts of CALABARZON and Metropolitan Manila before and after the phreatomagmatic eruption were compared. The presence and extent of heavy ash over the study area were identified after the image classification using the Support Vector Machine (SVM) algorithm. A decrease in vegetation cover and built-up areas was also observed. Upon validation, an overall accuracy of 91.4562 and Kappa coefficient of 0.8833 were achieved for the post-eruption ashfall extent map, exhibiting the potential of Diwata-2 imagery in monitoring volcanic eruptions and similar phenomena.


Background of the Study
For the first time in 43 years, the Taal volcanolocated in Talisay and San Nicolas, Batangas, Philippineserupted last January 12, 2020, at 2:30 pm PST. The explosive phreatomagmatic eruption yielded a giant plume of ash with a height of about 15km (Jing et al., 2020). A few hours later, plumes of hazardous particles and gases were released into the atmosphere. These were transported northward by high-altitude winds, causing an ashfall over the CALABARZON (Cavite, Laguna, Batangas, Rizal, Quezon) region and even several locations in Metropolitan Manila that are 42 to 90 km away from the Taal volcano (Leung et al., 2020). As monitored by the Philippine Institute of Volcanology and Seismology (PHIVOLCS), perilous gas emissions and earthquakes were brought about by the eruption, prompting the national government to impose a mandatory evacuation warning to the locals living within the 14-kilometer radius (Tripathy-Lang, 2020). The quick dispersal of the toxic gases made it quite hazardous for on-ground rescue teams and concerned agencies to identify the extent of the ashfall. It is therefore imperative that in the wake of catastrophic events like this, field surveys must be complemented by an efficient remote assessment of the affected area. One of the ways to do this is through the use of satellite imagery for remote sensing.
Last October 2018, the Philippine Scientific Earth Observation successfully launched Diwata-2, a microsatellite equipped with Near-Infrared (NIR) High Precision Telescope (HPT) and Spaceborne Multispectral Imager (SMI). In addition to the Diwata-2 satellite being supported by optical payloads that are built to determine the extent of disasters, its images have also shown potential for assessment of volcanic activity (Sabuito et al., 2020). This study aims to further explore the capabilities of Diwata-2 satellite imagery in ashfall dispersal mapping using the Support Vector Machine classification (SVM) algorithm. It is a supervised machine learning algorithm that employs the Kernel Trick technique and has been noted to achieve a higher level of accuracy compared to other established classification methods in remote sensing applications (Mountrakis et al., 2011).

Objectives
The objectives of this study are the following: 1. To perform corrections and classification on the satellite imagery to show ashfall extent; 2. To validate the presence of ashfall in the generated ashfall extent maps using ground reports; and 3. To perform an initial disaster assessment on the ashfall-affected areas in CALABARZON

Significance of the Study
Natural disasters incur many fatalities and economic losses for the ASEAN region, especially for developing countries such as the Philippines. The wide geographic stretch of incidences and increasing frequency of disasters require vulnerable nations to enhance their readiness and emergency response capacity. After catastrophic natural events such as floods, earthquakes, and volcanic eruptions, the main challenge is to assess the damage and destruction on-ground. Conducting fieldwork and sending survey teams to the damaged areas for disaster assessment tend to be costly, risky, and time-inefficient (Dashti et al., 2014).
Monitoring a high-speed disaster such as the dispersal of hazardous particles from a volcanic eruption using satellite remote sensing aids in more efficient disaster response and management. This will provide local government units a preliminary overview of which areas are more or less affected by ashfall so that appropriate relief plans are executed and the sites of evacuation centers are properly identified. The findings of this study also lead to a better understanding of the capabilities of the Diwata-2 microsatellite in observing calamities that strike the Philippines.

Data and Software Used
One of the objectives supported by the optical payloads carried by the Diwata-2 satellite is the determination of the extent of damages from disasters ("PHL-Microsat," n.d.). The satellite images used for this study were obtained from the Ground Receiving, Archiving, Science Product Development and Distribution (GRASPED) project under the Space Technology and Applications Mastery, Innovation and Advancement (STAMINA4Space) Program of the Advanced Science and Technology Institute (DOST-ASTI).
Images taken before (January 06, 2020) and after (January 27, 2020) the Taal volcano eruption were used for satellite data processing. A Diwata-2 Spaceborne Multispectral Image (SMI) with Top-of-Atmosphere (TOA) Reflectance was used as input for this workflow. The software utilized in this study to process these satellite images are Environment for Visualizing Images (ENVI) and Quantum GIS (QGIS).

Atmospheric Correction
The process of removing atmospheric effects to produce surface reflectance values, which can significantly improve image interpretability and use, is known as atmospheric correction (Yu et al., 2002). Because the Earth's atmosphere is constantly changing, which affects several surface variables such as land surface temperature, vegetation health, and change detection, it is critical that a satellite image be subjected to an atmospheric correction before processing.
The Diwata-2 satellite images were processed for this section using the Dark Object Subtraction (DOS) method in QGIS. In the DOS method, a specific value is subtracted in the top-ofatmosphere (TOA) reflectance. It is primarily assumed in this technique that there is a pixel in the satellite image that corresponds to the surface reflectance value in Equation 1, such that the TOA reflectance above that certain pixel agrees with Equation 2.
where ⍴surf = surface reflectance value where ⍴TOA = top-of-atmosphere reflectance ⍴atm = atmospheric reflectance The histogram of the Diwata-2 image was computed. For each of the 9 bands, the value of the TOA reflectance where the histogram starts to peak was recorded. In the Raster Calculator of QGIS, the values recorded were subtracted from the respective bands of the image.

Training Data Collection
The training data used were based on the preprocessed satellite images, applied with visual context, and are in reference to reports and a priori information about the volcanic eruption. Land features such as cloud, vegetation, and ash that were visible were assigned to corresponding classes for each image.
A new shapefile layer was created and training polygons were added to it while classifying the various land features on the atmospherically corrected images. The classes that were used to categorize the features on each satellite image are listed in Table  1. Samples of the training polygons that were included in the Regions of Interest (ROIs) for each class are also shown. The 2:30 pm color image of DIWATA-2 was used as the base layer in creating the polygons because it exhibits more clearly the vegetation, heavy ash, and cloud shadow on the land. The shapefile layer of the newly-created training polygons was used as the input for Training Areas for the Support Vector Machine classification in QGIS.

Support Vector Machine Classification
The Support Vector Machine classification was used on the satellite images using ENVI, specifically with its Radial Basis Function (RBF). Figure 2 depicts the parameters used for the SVM classification of the DIWATA-2 satellite image taken on January 27, 2020. The RBF Support Vector Machine classification took into account two hyperparameters: cost (C) and gamma. The former is the trade-off between training error and classifier generalization capability, whereas the latter determines the amount of influence that each training data has on the decision function.
The cost parameter is indirectly proportional to the margin of the decision boundary that separates different classes. The lower the value of C is, the higher the number of possible misclassifications. Meanwhile, a smaller value of gamma yields to more generalized regions that separate the classes. These two must be simultaneously optimized because a very large value of one parameter will make the effect of the other negligible (Yildirim, 2020).
The training data that were collected previously were used for the SVM classification in ENVI. Apart from the gamma and C values, the rest of the parameters were set to default. The most ideal C and gamma parameters that will not result in overfitting and underfitting are 100 and 0.0058 -0.0283, respectively (Sothe et al., 2020). C = 100 and gamma = 0.01 were used for this study. As a result of the Support Vector Machine classification, the satellite images were classified according to the set categories based on the regions of interest. The symbology of the classes was then changed intuitively according to the specified class.

Potential Ashfall Extent Map
For the post-eruption satellite image, the heavy ash class that was derived from the previous SVM classification was then retained while the rest of the classes were set to zero opacity. It was then overlaid on different base maps such as Google Satellite, Google Maps, and OpenStreetMap using the QuickMap Services tool in QGIS. The cities and towns that were covered by the classified ash were observed.
Afterward, image templating was done to produce map layouts of the potential ashfall extent on the satellite images. The resulting pre-eruption and post-eruption ashfall extent maps were compared and analyzed.

Accuracy Assessment
To validate the coverage of the classified ashfall over certain areas in CALABARZON, ground reports were collected and their coordinates were retrieved. The reported locations were matched with the cities and municipalities with heavy ash as identified on the potential ashfall extent maps.
Metrics were computed to quantify the thematic accuracy of the image classification. A confusion matrix was generated using ground reference data from Google Earth. The points used for the classes in each image were different due to the dissimilarity of the acquisition dates and area coverage.
The overall accuracy is the measurement of the ratio of the correctly classified pixel values to the total number of values and expressed as a percentage. On the other hand, the kappa coefficient (κ) measures the agreement between classification outcome and ground reference values. The kappa coefficient is a value from 0 to 1, where 1 is the representation of a perfect agreement. Their computations are shown in Equations 3 and 4.

=
(3) The user's accuracy and producer's accuracy for each class in the pre-eruption and post-eruption images were also computed. The user's accuracy is a quantitative measurement of how reliable the image classification is and how well it represents what is really on the ground; the producer's accuracy defines the probability that a reference pixel is correctly classified.

Preprocessing
To remove the effect of the atmosphere on the radiation measured by the satellite sensors, an atmospheric correction was applied to the images. The Dark Object Subtraction (DOS) method was used and a visual comparison of the satellite images before and after the preprocessing is shown in Figure 2. The atmospheric correction on the satellite data brought about significant changes in both the spectral signatures and visual quality of the images. The DIWATA-2 images were clearer and the haze was greatly reduced.

Image Classification
The classification outputs for the images before and after the Taal eruption are shown in Figure 2. It can be observed that no land cover in the pre-eruption image was classified as heavy ash. This was validated by referring to news reports indicating the timeline of the events that occurred around the occurrence of the eruption of the Taal volcano in 2020. According to Rappler (2020), the volcano started spewing ash only on the afternoon of January 12, 2020. Its alert status was raised to Level 2 at 2:30 pm, then to Level 3 at 4:00 pm, and Level 4 at 7:30 pm on the same day. By then, the ash column has already reached up to 15 kilometers, coupled with volcanic lightning and lava mud (Rappler, 2020). For the post-eruption image, it can be seen that the heavy ash was concentrated on and around the Taal volcano. It can also be observed north of Batangas, over the provinces of Cavite and Laguna.

Post-processing
Some of the pixels classified as heavy ash were located over cloud shadows, on the edges of the image, or along the shoreline of Batangas. Post-processing was applied to further smoothen the image. The result is shown in Figure 4.

Visualization and Validation
The potential ashfall extent map of the provinces circling the eruption site is shown in Figure 5. It can be observed that the dispersal is prominent over areas to the north and west of the volcano. Because of the boundaries of the image captured by the Diwata-2 satellite over the CALABARZON region, the presence of heavy ash over the province of Quezon and some parts of Cavite were not recorded.
The dispersal of ashes and other hazardous particles after a volcanic eruption is rapid and while this potential ash extent map is not representative of their location at any given point of time, it shows the areas where there is a high agglomeration of heavy ash. If complemented with meteorological findings such as wind direction and other weather conditions, this map is likely to provide insights on where the ash is headed next. The corresponding areas for each land cover as measured from the classified images are shown in Table 2. Although the preeruption and post-eruption satellite images do not cover the same extent of the CALABARZON region, there is a significant change in the area that can be observed for the vegetation and heavy ash classes. More insights were derived from the user's and producer's accuracies that were computed for each class, as shown in Table  4. Most of the accuracies of both classified images are highly satisfactory, suggesting that the classifications are quite close to the ground truth. Some outliers include the user's accuracy for the built-up areas and the producer's accuracy for the heavy ash class. The latter yielded the lowest accuracy, which can be attributed to the low spatial resolution of the Diwata-2 satellite which is at 127 m. Another method that was done to validate the map was to collect reported statements and photos about the presence of heavy ash in specific locations in CALABARZON. These were then matched with the potential ashfall extent maps derived from the satellite images. An example is shown in Figure 6, showing an ash-covered Ferris wheel at the Sky Ranch amusement park in Tagaytay, which matches with the presence of ash as indicated by the map. Figure 6. Validating the potential ashfall in Tagaytay City with ground reports The caption on the photo says "A Ferris wheel is covered with volcanic ash in a park in the city of Tagaytay, on January 14, 2020 (Taylor, 2020)," in an online report by The Atlantic Media, a Washington-based multi-platform publisher.
Another match with a ground report is shown in Figure 7. The photo is included by CNN Philippines in one of their online news reports on the aftermath of the Taal volcano eruption in nearby affected areas. The aerial photo shows ash-covered roads, trees, and rooftops in a residential area located in Barangay Buso Buso, Laurel, Batangas (Regan, 2020). No specific road or establishment was mentioned in the report, so the centroid of Barangay Buso Buso was considered as the area of heavy ash. The ash extent map derived from the DIWATA-2 image shows that potential heavy ash covers a portion of the Municipality of Laurel, including its centroid. With the accuracy and validity of the generated potential ash extent map, there are initial findings that can be derived to aid in disaster response such as relief efforts, evacuation, and rescue operations.

CONCLUSION
This study detected the presence and dispersal of ashfall over parts of the CALABARZON region and Metropolitan Manila, a week after the explosive eruption of the Taal volcano. Changes in land cover were also observed in the decrease of vegetation and built-up areas in a span of 21 days. A total area of approximately 302.53 km 2 was covered with heavy ash, which was concentrated in locations north of Talisay, Batangas. The direction of this dispersal can be confirmed by the aforementioned reports of high-altitude northward winds in the area on the same date. Furthermore, the high overall accuracy and Kappa coefficient of 91.4562 and 0.8833, respectively, signify that the classification of the post-eruption ash extent is sufficiently validated. It can be concluded that the method used in this study for mapping the dispersal of ashfall is adequate for remote disaster monitoring and assessment purposes. This study can be applied and utilized by national agencies, local offices, and stakeholders for rescue operations, policymaking, and rehabilitation purposes.

RECOMMENDATIONS
For the satellite data pre-processing, the researchers recommend using more complex atmospheric correction tools for the Diwata-2 satellite image. For the image classification of the satellite data, it is recommended to input other values of gamma and cost parameters in the Support Vector Machine classification and to also explore other classification methods and test their capabilities in classifying ash. Due to the conduct of this research during the peak of the pandemic, it was not plausible to coordinate with local government agencies for more field reports and data on the 2020 Taal volcano eruption, so it is suggested to add more sources for the ground truth data. The inaccuracy of the classification of some features could be accounted for by the low resolution of the Diwata-2 satellite data. High-resolution satellite images could be explored to identify heavy ash with the methodology utilized in this study.