Intra-Annual Identiﬁcation of Local Deforestation Hotspots in the Philippines Using Earth Observation Products

: Like many other tropical countries, the Philippines has suffered from decades of deforestation and forest degradation during and even after the logging era. Several open access Earth Observation (EO) products are increasingly being used for deforestation analysis in support of national and international initiatives and policymaking on forest conservation and management. Using a combination of annual forest loss and near-real time forest disturbance products, we provide a comprehensive analysis of the deforestation events in three forest frontiers of the Philippines. A space-time pattern mining approach was used to map quarterly deforestation hotspots at 1 km pixel size (100 hectares), where hotspots are classiﬁed according to the spatial and temporal variability of the 2000–2020 deforestation in the study area. Our results revealed that 79–81% of the hotspots overlap with primary forests and 27–29% are inside the state-declared protected areas. The intra-annual analysis of deforestation in 2020 revealed an alarming trend, where most deforestation occurred between the 1st and 2nd quarter (92–94% in hotspot forests; 87–89% in non-hotspot forests), highly overlapping within the slash-and-burn farming season. We also found “new” hotspots (2020) formed mostly from landslide scars and partly from selective logging, the latter is believed to be underesti-mated. Our study paves the way for rapid and regular assessment of the country’s deforestation, useful for the respective environmental institutions who convene several times a year. Moreover, our ﬁndings assert the imperative of alternative livelihoods to upland farmers, efﬁcient forest protection activities, and even the mitigation of landslide risks.


Introduction
Forests provide a broad range of ecosystem services including the mitigation of climate warming and its adverse consequences [1] and the regulation of water to sustain croplands and even mitigate potential floods [2]. Forests also provide aesthetic, cultural and recreational services for the wellness of humans and are home to diverse groups of flora and fauna. The tangible benefits from forests are obviously timber and also showed that alert systems could trigger immediate interventions in forest management and conservation by environmental stakeholders [21].
These EO products play a role as inputs to spatially identify hotspots of deforestation. One goal of identifying these hotspots is to prioritize forest areas that need immediate interventions i.e., forest protection. Mapping hotspots can either be through a subjective judgment by experts or using geospatial models [22,23]. Gandhi and Jones et al. (2019) [24] subjectively used the rate of mangrove changes from a set of country data to determine which among Southeast Asian countries are mangrove deforestation hotspots. Mapping hotspots using models may require spatial and temporal (spatio-temporal) data to identify hotspots based on the rate and extent of deforestation over time. Jardeleza et al. (2019) [25] used zonal statistics to identify 2010-2015 hotspots in Philippine forests based on tree cover from different forest covers including a geo-simulated version. Harris et al. (2017) [23] used a space-time pattern mining called Emerging Hotspot Analysis (EHA) to define statistically significant hotspots as a byproduct of the GFC forest loss data. Reddy et al. (2016) [22] did a similar hotspot mapping where each hotspot pixel has a confidence interval while considering the trends in deforestation since 1930, while Sanchez-Cuervo and Aide (2013) [26] applied EHA to assess hotspots in protected areas in relation to hotspot drivers.
No intra-annual deforestation hotspot products have yet been developed and in this sense, the forest disturbance alert products can be useful. The current available EO product for deforestation hotspots is the global EHA-based map [23], being updated every 1-2 years depending on the release of the annual forest loss data [27]. Hotspots can be seasonal, driven by illegal activities usually implemented when the weather is favorable in certain months. This seasonal activity cannot be easily revealed when looking at annual hotspot maps and forest loss per se. Regular hotspot mapping can also help evaluate any progress of interventions such as forest protection. This evaluation can take place during the regular assembly of environmental bodies several times a year.
Here we examine how annual and near-real time EO products can be used to regularly identify deforestation hotspots in the Philippine forests and understand their spatiotemporal trends in relation to the deforestation drivers towards an updated, efficient and practical forest protection and management.

Study Area
Our study sites are three forest frontiers in the country that encompass eleven provinces. The first is in the Sierra Madre, one of Southeast Asia's longest mountain ranges that extend within nine provinces. This mountain range is home to the last remaining primary forests in the Philippines, with very high biodiversity and the most number of protected areas. In these mountains, selective logging and agricultural expansion have persisted with alarming rates [28]. The second study site is Palawan province, an elongated island and the largest province in the country. Palawan is known to have primary forests with high endemism of flora and fauna, which are threatened by an increase in active mining operations, slash-and-burn farming, oil palm expansion and timber poaching [29]. The last study site is Apayao province, known to have high endemism of fauna including the rare Philippine eagle, roaming within three highly forested mountains in the province. See Figure 1 for the study area map.  [30]) and the state-declared protected areas are also shown.

Method Overview
We used Global Forest Watch (GFW) [27] as a starting point to gain an overview of deforestation in the Philippines over time. After gathering the needed deforestation data, we converted them to space-time data cubes as inputs to the Emerging Hotspot Analysis (EHA) to identify local hotspots and their classes, which were then interpreted spatially and temporally. See Figure 2 for the overview of the main steps of this study.

Deforestation Data Sources
We used EO products that are widely used to analyze deforestation in the local [11][12][13] and global context [31,32]. First, the GFC forest loss data from the years 2000 to 2019 were accessed from the GFC data portal. This dataset labels a 30 m pixel as "loss" if the pixel before the loss is a vegetated pixel greater than 5 m tall followed by a "removal", which can be caused by man-made and natural factors [18]. It is cautioned though by the data producers that forest loss is not always equivalent to deforestation. Second, we used GLAD forest disturbance alerts which were accessed using Google Earth Engine (GEE). The GLAD alert system is also a Landsat-derived product (30 m) and can detect disturbances with at least 50% tree canopy removal [20]. The main caveat for this product pertains to the confirmation of the alerts. Each alert pixel should be "confirmed" (detected again after the first detection) and the date of alert would be recorded as the date of said confirmation. Such confirmation can be delayed due to cloud interference. While the next alert product we used, the RADD system, is unaffected by the cloud-induced confirmation delay, the RADD system labels pixels slightly different than the GLAD system as RADD uses the first detection date as the alert date while retracting alert pixels if they are not confirmed within 90 days. Its producers have cautioned of possible undetected small-area disturbances (<0.1 ha). RADD alerts were also accessed via GEE. See Table 1 for additional details about these three deforestation data sources.

Hotspot Mapping and Assessment
To assess the 2020 deforestation scenario spatially and temporally, we generated quarterly hotspot maps that account for the potential deforestation events from 2000 to 2020 using the annual forest loss data (2000-2019) in combination with each alert product (2020), resulting in two hotspots map versions. We included areas with >30% of the GFC's tree cover density (TCD) product. This threshold is slightly above the forest mask of the global hotspot data (25% TCD, [23]) and the maximum threshold stipulated in the United Nations Framework Convention on Climate Change [34]. Such forest mask would likely capture disturbances from upland agriculture encroachment. We then resampled the RADD pixels using nearest neighbour method from 10 m to 30 m for consistency in spatial scale i.e., to avoid having three times more pixels than the other products because of the pixel size difference. We also used the same forest mask for the deforestation alerts so we masked out GLAD alerts outside the primary forest baseline year 2019 [30] included in the RADD product. Herein, we call the pixels from the annual forest loss and disturbance alerts as deforestation pixels.
For hotspot mapping, we used EHA (via ArcGIS) being the widely used method from local to global hotspot mapping using EO products. The spatial scale of the EHA matters and we aimed for 1 km × 1 km hotspot bins to preserve small potential hotspots (despite a long computation time), deemed useful for forest protection at local jurisdictional scales. For each bin, we assigned a 3 × 3 window as neighboring bins. The consistency of these hotspots is routinely tested at regular time intervals and we aimed for a 1-year period while updating the 2020 deforestation pixels quarterly (as sufficed by the alert products), hence producing four hotspot maps in 2020.
The technical steps of EHA first involve a clustering method based on the density of deforestation pixels inside the 1-km bins and among neighboring bins, followed by a comparison of such density for each bin relative to all deforestation pixels in the study area. If a bin is found to be relatively dense, a high standard deviation (Z-score) is expected for that bin as evaluated by the Getis-Ord Gi statistic [35]. The second EHA step is a temporal assessment of the deforestation trend for each bin and is performed using the Mann-Kendall trend test to determine whether the trend among bins is increasing or decreasing. Statistical significance (p < 0.05) from these spatial and temporal evaluations would assign a 1-km bin as a hotspot with corresponding classifications described in Table 2. Hotspots are hereby defined as local areas with "statistically significant deforestation". Table 2. Hotspot classes and their description. The time step represents 2001-2020 years (2001-2019 = GFC loss data; 2020 = alert products). The x rows depict "statistical significance of deforestation" per year or simply those classified by EHA as a "hotspot". New 1-km hotspot bins that are statistically significant (>95%) only in 2020 x Consecutive 1-km hotspot bin that are statistically significant (>95%) only in recent years x x x Sporadic 1-km hotspot bins with discontinuity in statistical significance over the years (on-and-off hotspot) Intensifying Same with persistent but recent years have increasing deforestation Non-hotspot 1-km hotspot bin without statistical significance in the past 2 decades Hotspot classes were mapped, counted and analyzed for the year 2020. Hotspots that overlap with the primary forest in 2019 and the government-declared protected areas were also tallied per hotspot class and province. Including these three boundaries would enable joint assessment and interventions by both the national and local governments. The provincial boundaries were obtained from Global Administrative Areas (GADM) while the protected areas were provided by the Biodiversity Management Bureau of the Department of Environment and Natural Resources (BMB-DENR). The primary forest boundary used in the RADD alerts was adapted.
Temporal assessment of the hotspot and deforestation areas was also initiated to investigate deforestation peaks in 2020. Using the 30 m alert pixels (GLAD and resampled RADD), we graphed the accumulated and per month deforested areas. This analysis was extended until April 2021 to assess whether there is a seasonal trend every 1st quarter. We then assessed the deforested areas per hotspot class and outside hotspots (forest nonhotspots), from one quarter to another in 2020.
We finally related the hotspot classes to their main drivers using Sentinel 2 and Planet satellite imageries viewed at GFW, and Google Earth Pro particularly its terrain functionality. Both platforms are useful to easily navigate the landscape and distinguish between two drivers based on spatial patterns i.e., landslide (upslopes) and slash-and-burn (foot slopes). Further validity of the drivers was confirmed by multiple local experts based on their familiarization and access to on-site reports. Lastly, local and international news were collated and summarized in Table A1 as proof that deforestation existed in 2020.

Hotspot Maps and Analysis
Hotspots that account for the deforestation from 2000 to 2020 were found to exist in the whole study area for both maps that use the RADD and GLAD alerts. These hotspots shown in Figure 3 are a mix of the three hotspot classes located within the remaining primary forests of the country. The detailed breakdown of these hotspots is shown in Table 3 where 1978-2246 hotspots within nine provinces were classified as hotspot areas. Alarmingly, 79-81% of the hotspots overlap with the primary forests in 2020; while 27-29% of them overlap with protected areas. A higher number of hotspots was identified when using RADD alerts (12% higher than GLAD alerts). Results from the two hotspot versions were particularly close for hotspot classes with on-and-off deforestation (sporadic hotspot) and those hotspots formed in recent years (consecutive hotspot). For the hotspots in 2020 (new hotspot), the hotspot map with RADD alerts was found to be higher by 69% than the map with its counterpart alert product. We found one persistent and one historical hotspot and we just merged them to the sporadic class.
The hotspot maps show the regions with the highest number of hotspots such as Southern Palawan, the mountains of Apayao province, and foot slopes of Sierra Madre mountains. The very dense red pixels in Southern Palawan are a combination of sporadic and consecutive hotspots. Notable also are new hotspots forming in the mid-part of the Sierra Madre mountain range and in the northernmost part of Apayao. Visible also in most provinces are scattered hotspots that are either sporadic or consecutive.   GLAD  RADD  GLAD  RADD  GLAD  RADD  New  57  189  100  100  72  40  Consecutive  765  802  72  73  25  26  Sporadic  1156  1255  82  83  27  29  Total  1978  2246  79  81 27 29 The highest number of hotspots that overlap with PAs were found in the provinces of Palawan followed by Isabela and Quirino. These PAs are Mt. Matalingahan and Malamplaya Protected Landscapes in Palawan, and Quirino Protected Landscape and Northern Sierra Madre Natural Park in Sierra Madre. Relatively, the provinces with the highest % hotspot in PAs (almost proportional hotspot count in both PA and non-PA) were ranked and enumerated as follows: Quirino, Isabela, Rizal and Palawan. See Figure 4 for the provincial-level results.

Deforestation in the Year 2020
The deforestation analysis in 2020 is shown in Figure 5a as cumulative deforested areas per week (default temporal resolution of the alert products) together with its monthly aggregates in Figure 5b. Overall, the highest deforestation was recorded in dry months from February to May (peak in April) resulting in a sharp increase in cumulative deforestation. This trend gradually decreased afterward and hit the lowest alerts around October only to increase again until March 2021, hence suggesting a seasonal trend.
There was a slight difference in the trend between the two alert products. While both had their highest detection in February and April, the RADD alerts were relatively higher at the beginning of the year. These alert pixels can be the forest loss pixels in the late months of 2019 that were not accounted for when the forest mask used for RADD was created. The GLAD alerts fluctuated from August to October, which can be late confirmations of earlier disturbances i.e., in June as the start of the wet (and cloudy) season. In May and July, the GLAD alerts were slightly higher than the RADD alerts.
For the quarterly summary of deforested areas (Figure 6), aside from the higher detection by RADD, we found two key observations. First is the highest deforestation between the 1st and 2nd quarter of 2020 both in forest areas identified as hotspots (92-94%) and forests that are non-hotspots (87-89%), with the exception of the GLAD alerts inside new hotspots. Second is the relatively higher alerts outside hotspots (53-76% of all alerts), which are likely those small deforested areas scattered all over the study area.  Spatial patterns were also observed in the development of hotspots within the year. The example shown in Figure 7 shows that sporadic hotspots increased by 400 ha from 1st to 2nd quarter in 2020, going uplands and within the borders of the primary forests, see Figure 7a,b. This is also depicted in Figure 7c,d where the actual conversion indicates that these hotspots are likely to be driven by slash-and-burn farming judging by the soil color. Another example of hotspot development (500-ha increase from 1st to 3rd quarter) is shown in Figure 8 where landslides are the main driver of deforestation in the mountainous and forested parts of northern Sierra Madre. Strong typhoons also passed through the area in the 2nd and 3rd quarters of 2020. The example hotspots were found to overlap in different municipalities.

Hotspot Mapping and Deforestation Detection
Our hotspot maps revealed that 79-81% of the hotspots overlap with the primary forests and around 27-29% are inside the declared protected areas. These hotspots are formed in clusters while some are scattered within the provinces' forests. These numbers are alarming and relatively higher than the global hotspot map mainly due to the smaller spatial scale (hotspot pixel size) in this study, which is suited for local to regional analysis. Nevertheless, there is still a high overlap between our hotspot maps and the global hotspot particularly in the dense hotspots within Southern Palawan and Apayao mountains mainly driven by the same forest loss input (2000-2019) and hotspot mapping approach (EHA). The map producers of the global hotspot also advised adjusting the spatial scale of hotspot bins to preserve small hotspot patches for local-scale forest protection i.e., from 10 km to 1 km spatial bins [23]. Another study also used the same annual forest loss data and they revealed that protected areas inside Sierra Madre and Palawan experienced the highest forest loss rate from 2010 to 2012 [11]. Using again the loss data and a time series vegetation index, the continuity and severity of deforestation in Sierra Madre were reported to cause a net forest loss even if the potential forest gains from the national reforestation program in 2011 are accounted [13].
Because of the alert products, 57-189 new hotspots (in 2020) and 765-802 consecutive (recent years) hotspots were revealed; while the rest are sporadic hotspots (on-and-off from 2000 to 2020), which would be unidentified without the annual loss data. Though forest loss was reported to be increasing from 2017 to 2019 [13], the relatively high number of new (and consecutive) hotspots could also be linked to the higher deforestation pixels provided by alert products compared to the annual loss pixels. It was reported that the 2000-2019 GFC forest loss may underestimate the actual forest loss in a given year, particularly in the tropics [36]. This potential omission error (missed detection) is mainly caused by the limitation of optical satellites to secure cloud-free images over large scales within the year as well as the inability of the product to detect small-scale disturbances [36]. This limitation can also be the reason for the higher number of hotspots detected when using RADD (radar) compared to GLAD (optical). Such comparison was demonstrated in Congo where higher detection was reported using the radar-based system in comparison with GLAD and the GFC loss data [33] and in Peru when also considering forest edges to detect deforestation [37]. Though it was reported also that the early 2020 RADD alerts in Asia may be overestimated due to the forest baselining of the RADD product [33]. On the other hand, there can still be potential commission errors (false positives) for an alert product since they detect losses based on single-date images instead of image composites [21]. Nevertheless, both (high-confidence) alert products would be useful for future hotspot mapping i.e., GLAD alerts in areas outside primary forests, where radar signals decrease in sensitivity e.g., flooded forest, mangroves and peatlands. The complementation of the two alerts can be further explored using their average for a conservative estimate of deforestation and hotspots.
Those deforested areas unclassified as hotspots (53-76% in 2020) are likely deforestation events that are small and less intensive within a 100 ha spatial scale (the hotspot bin size). A prominent example is selective logging in selected small forest areas to supply domestic needs of lumber. Moreover, upland farmers routinely practice slash-and-burn (locally known as Kaingin) and charcoal making in several small areas, often practiced as a household effort and have been increasing due to upland migration [38]. Technically, the classification of the deforestation pixels as non-hotspot is also due to the EHA itself. The hotspot mapping algorithm only looks for statistically significant hotspots relative to the whole forest area, where most are still naturally protected because of inaccessibility at the very least. The local statistics in a highly disturbed area with a size of 100 ha would be very different from the global statistic (study area), e.g., dense deforested pixels from large landslides were classified as new hotspots while Kaingin areas by households belonged to non-hotspots. Many deforested pixels would be unclassified as hotspots once the hotspot bin size increases. See Figure A1 for omitted hotspots when increasing the hotspot size to 1000 ha instead of 100 ha. The hotspot algorithm also considers the longevity of forest disturbance and those areas with dense disturbances in limited time, e.g., 1-3 years (except recent years) would likely be statistically insignificant and hence non-hotspot. Potential examples of this case would be the conversion of forests to oil palms and surface mining in Palawan.

Deforestation in the Year 2020
A clear seasonal trend of deforestation was observed in 2020, highest from February to May of the year (peak in April). That period is also the time of slash-and-burn season as verified by the key informants from environmental task force members at the sub-municipal (barangay) level. They further shared that aside from Kaingin and charcoal making, timber poaching still exists. Although the informants also mentioned that loss of jobs during the COVID-19 lockdown has driven the increase of deforestation i.e., halt of eco-tourism in the case of Palawan, and since most people (including forest rangers) were prohibited outside their houses, attributing the increase of deforestation due to the lockdown would need a separate analysis. Once GLAD alerts from 2015 are made publicly available, we would assess whether the spike from 1st to 2nd quarter in 2020 is statistically significant. Nevertheless, we summarized online reports that deforestation existed in 2020 in Table A1. Six published news online reported illegal activities from different parts of the study area, consistent with the information from local experts.
The deforestation in 2020 hints that previously deforested areas in recent years have been expanding as evident in the high number of consecutive hotspots. A recent study used GLAD alerts to conduct an inter-country comparison of deforestation in 2020 and it was concluded that the continuation of deforestation from 2019 to 2020 is a more reasonable root cause of the 2020 deforestation statistics than the pandemic lockdown effect [39].

The Usefulness of the Hotspot Maps
The quarterly maps of hotspots can be helpful to guide the necessary interventions, depending on the kind of hotspot and their drivers since each driver can have different spatial and temporal patterns [40]. While the sporadic hotspots reasonably consist of mostly Kaingin areas due to the cyclic nature of such practice, the consecutive hotspots, as recently mentioned, can be signs that these cultivated lands are encroaching, likely a consequence of upland migration and high crop yield demands [38]. Similarly, those new hotspots that would exist as a single 1-km pixel (without neighboring pixels) would likely appear in places nearby the first single pixel as an indication of selective logging. This kind of illegal logging is an organized mission by poachers that extends from weeks to months until the target quota of logs is reached and would result in forest degradation [28]. The new hotspots found mostly in Sierra Madre and Apayao (see Figure 3) are formed mostly from landslide scars, reasonably because these regions are mountainous and within typhoon belts. The hotspots in Southern Palawan are mainly sporadic and consecutive, suggesting that deforestation is mainly Kaingin-related as such practice has been a persistent livelihood source in the province for decades [41].
Consider these hotspots an early warning so authorities would know where and when to intervene-which saves time and resources-instead of considering all deforestation pixels. First, the new hotspots mostly from landslides during typhoon months are valuable to the government in the context of disaster risk reduction and management (DRRM) and landslide rehabilitation. Here we provide the first provincial-scale information about the severity of landslides after typhoons. Secondly, potential new hotspots attributed to illegal selective logging should be given priority in forest protection. Several teams can be deployed to monitor strategic checkpoints and other potential transport routes of the logs with regard to the hotspot location. Hotspots sometimes extend over two municipalities so local coordination is of utmost importance for mobility and security reasons. The use of high-technology forest protection tools such as the Lawin app and flying of Unmanned Aerial Vehicles (UAV) for surveillance are additional options given the inaccessibility of these timber poaching zones. Thirdly, the slash-and-burn practice has been a challenge for the past three decades and is still very far from its ideal transition to sustainable farming. Those kaingeros responsible for any identified large and expanding kaingin areas should be the priority recepients of any upland livelihood programs. Lastly, any new hotspots related to mining should be closely monitored now more than ever because of the recent Executive Order 130 of the government that would allow new mining explorations and permits. The impact of any interventions can be evaluated during the quarterly assembly of multi-stakeholder environmental bodies. Evaluations should be in line with the anticipated forest gains from the existing national reforestation program to determine whether such program would indeed result into a net forest gain [13].
The number of hotspots inside protected areas and primary forests are key inputs to Forest Land Use Plans of municipalities. A similar spatial analysis of hotspots with and without tenurial instruments can also be initiated, particularly those under the management of upland communities.

Updating of the Hotspot Maps
Since our mapping workflow is mostly automated and reproducible, we can efficiently update and modify hotspot maps if the need arises, i.e., adjust time intervals per wet and dry season or modify the 2000 forest baseline. The same routine where alert products are used only in the latest year (e.g., 2021) will be applied for future hotspot maps. We can also produce mangrove deforestation hotspots given a time series mangrove extent input. Similarly, forest degradation hotspots would need to be further integrated into our current hotspot mapping by preserving the RADD alerts' spatial resolution while resampling the GFC and GLAD data to 10 m, and further adjusting the hotspot size from 1 km to 250-500 m. This, however, would still not account for the disturbances undetected at 10 m pixel size (e.g., small trees for charcoal making) not unless an EO product of degradation alerts comes out or we use time series above-ground biomass as input to EHA instead.
Lastly, the validation of deforestation drivers by local experts especially for timber poaching zones can also be complemented by platforms such as GFW and Google Earth Pro to be more certain on the hotspot drivers (see Figure A2). This method of validation is even more promising now that Planet images with very high spatial and temporal resolution are reported to be available openly.

Conclusions
Quarterly deforestation hotspot maps in 2020 were created at 1 km pixel size using the Emerging Hotspot Analysis (EHA) and a time series input of deforestation data from Earth Observation products. The annual forest loss time series is useful to identify sporadic hotspots, likely slash-and-burn farming areas, while the forest disturbance alerts are helpful to reveal consecutive hotspots that indicate encroachment of these cultivated lands, and new hotspots formed mostly from landslide scars. Relating the hotspot classes to the deforestation drivers is best attained with the help of local experts e.g., validation of timber poaching zones. Hotspots solely driven by timber poaching (forest degradation) would be integrated into the next hotspot versions. Succeeding work for this may include the use of the EHA inputs with 10 m pixel size (instead of 30 m), creating smaller hotspot bins, and substituting the forest loss inputs with time series above-ground biomass. Regular updating and assessment of hotspot maps are expected to happen prior to or during the quarterly assembly of environmental institutions, with close attention to any potential hotspot caused by the anticipated new mining areas. Once local hotspots are identified, interventions can be more proactive in order to prohibit hotspot expansions, both preserving forest ecosystem services and government resources. to all Filipino forest rangers who died in service, and to our friend, the late Jun Sasuman, Municipal Environment and Natural Resources Officer of San Vicente, Palawan.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Figure A1. The comparison between 1 km (left) and 10 km (right) spatial sizes of hotspots. Spatial detail and local hotspots are preserved in the former which is suited for local enforcement of forest protection activities.