Spatio-temporal analysis of remotely sensed forest loss data in the Cordillera Administrative Region, Philippines

The Cordillera Administrative Region (CAR) in the Philippines is among the last forest frontiers in the country and is also home to 13 major watersheds in Northern Luzon that supply irrigation and hydroelectricity to other regions. However, it is faced with the deterioration of the quality of its watersheds due to forest loss driven mainly by agricultural expansion and illegal logging. Thus, this study was conducted to analyze the spatial and temporal patterns of forest loss that could serve as a basis for policy decisions. Also, this paper determined the strength of relationships using Pearson’s correlation coefficient (r) between forest loss and seven independent variables, which includes forest cover, agricultural areas, built-up, road network, and socio-economic data. This study utilized the Hansen Global Forest Change (HGFC), a Landsat-derived dataset from 2001 to 2019. Results revealed that 70,925 hectares (ha) of forest loss were detected with an annual deforestation rate of 3,744 ha/year across the region. Based on the validation, the accuracy of the HGFC data is 72%, but great caution should be observed when using the data with less than 0.2 ha due to very low accuracy. On a region-wide analysis, only the forest cover had a strong association with the forest loss with a computed r-value of 0.78. Conversely, on a provincial level, the explanatory variables had a strong to moderately strong correlation with deforestation. Hence, immediate, science-based, and sustained regional and multi-stakeholder efforts are necessary to conserve and protect the remaining forest cover in the region.


Introduction
Based on the Global Forest Resources Assessment (FRA) report in 2020, around 30% of the world's total land area are classified as forests [1], and it is home to over 80% of the Earth's terrestrial species of flora and fauna [2,3]. Forests play an important role in global ecosystem services and provide many benefits, including carbon sequestration, clean air and water, climate and water cycle regulation, and biodiversity conservation [4,5]. They also serve as a food source, medicine, and fuelwood for 2.6 billion people worldwide [6]. Moreover, the timber sector generates 54.2 million jobs and contributes $600 billion to the global economy [7], and around 1.6 billion people, the majority of whom belong to the world's poorest, rely on forests for their subsistence [8]. However, these forest areas, particularly in the tropical regions, are being depleted at an alarming rate [9,10], driven mainly by agricultural expansion [11,12], with an estimated annual deforestation rate of 10 million ha/year recorded globally between 2015 and 2020 [13].
Like many other developing countries in the tropics, the Philippines has undergone remarkable pressure from deforestation and forest degradation during and even after the logging era in the 1960s to '80s [14,15], due to forest product extraction, agricultural and infrastructure expansion [16]. Even a significant number of Protected Areas (PAs) in the country had a notable forest cover loss in terms of extent and rate over the last decades [17]. At present, forests cover around 7.014 million hectares (ha), or 23.38% of the country's total land area of about 30 million ha [18]. Most of the forested areas in the country are confined within three regions, which can be considered as the last forest frontiers. These are located in the Cagayan Valley or the Sierra Madre Mountains (Region 2) with 1.05 million ha, the MIMAROPA (Region IVB) with 0.95 million ha, and the Cordillera Administrative Region (CAR) with 0.81 million ha [18]. However, there has been no significant forest cover increase in these regions from 2001 to 2018 despite the country's reforestation efforts such as the National Greening Program (NGP) [19]. An alarming finding also reported that around 80% of the primary forests found in these regions are considered deforestation hotspots [20]. Thus, there is an urgent need to prioritize forest protection and conservation in these regions to halt further degradation and biodiversity loss in the country.
To address the impact of increasing forest loss on a national and global scale, the United Nation's (UN) Sustainable Development Goal (SDG) 15 (Life on Land) was formulated, which specifically aimed to ''protect, restore and promote sustainable use of terrestrial ecosystems, sustainably manage forests, combat desertification, halt and reverse land degradation, and halt biodiversity loss'' [21]. However, many nations are failing to meet their targets, as a result, these countries are falling short to halt deforestation which greatly contributes to the biodiversity loss crisis. The goal under the SDG 15 is to ensure the conservation, restoration, and sustainable use of terrestrial and inland freshwater ecosystems and their services, in particular forests, wetlands, mountains, and drylands, and to promote the implementation of sustainable management of all types of forests, halt deforestation, restore degraded forests and substantially increase afforestation and reforestation [22]. But despite the country's effort, the national review on the SDGs suggests that more work is still needed in protecting the forests and terrestrial sites in the country [23]. With this current scenario, supporting the government's monitoring effort of forest loss in the country, particularly the critical regions where most of the remaining forests are located, has never been so important. Aside from monitoring the spatial and temporal patterns of forest loss, identifying the drivers of deforestation are indispensable not only in supporting the SDG 15 but in planning sustainable forest management [24] and climate change mitigation and adaptation programs as well, such as Reducing Emissions from Deforestation and forest Degradation (REDD?) [25]. Though the ability to monitor forest loss is dependent on the availability of reliable data sets, as well as appropriate methodology and statistics [26].
With the advent of Geographic Information System (GIS) and remote sensing technology over the past years, regular forest loss monitoring on a national or global scale has become convenient [27]. Remotely sensed data, such as the Hansen Global Forest Change (HGFC) generated from high-resolution Landsat images, demonstrated its potential in monitoring forest cover changes [28,29]. The use of this data, through geospatial techniques and statistical analysis, will provide deforestation analysis that can be a basis for policy decisions towards protection and sustainable management of the remaining forests in the country, and in turn, halting biodiversity loss.
Here, the analysis of spatial and temporal patterns of forest loss was conducted in one of the critical regions which hold a significant portion of forests in the Philippines, the Cordillera Region. It is a non-coastal region situated in the northern island of Luzon in the Philippines and is home to 13 major watersheds in Luzon that supports four (4) of the biggest hydroelectric and multipurpose dams in the country. Also, these watersheds provide water for domestic, agricultural irrigation systems, and other ecosystem services to the region and its neighboring lowland regions. However, due to rapid agricultural expansion, uncontrolled timber poaching, including ''kaingin'' and forest fires, the region is currently under threat with the deterioration of the quality of its watersheds [30]. Hence, this paper aimed to provide local stakeholders and policy decision-makers with science-based spatio-temporal forest loss data and determined the annual deforestation rate on a provincial scale using remotely-sensed data. In addition, the present study evaluated the strength of relationships between forest loss and several spatially explicit thematic and socio-economic data variables using correlational coefficient analysis, which could enhance the public's understanding of the potential drivers of deforestation in the region.

Study area
The study sites are located in the six (6) provinces of CAR namely, Abra, Apayao, Benguet, Ifugao, Kalinga, and Mountain Provinces (Fig. 1). The region is characterized by mountainous areas with rolling to severely steep slopes and elevations that range from 3 m above sea level (masl) to 2,928 masl with its highest peak located in Mount Pulag located in Benguet, the third highest mountain in the country. The vegetative forest cover types of the region include dipterocarp forests typically located in the lowlying areas of Abra, Ifugao, Apayao, and Kalinga, molave forests, coniferous forests dominated by Pinus kesiya species, and mossy forests found in the higher altitudes of Benguet, Mountain Province, and Ifugao. At present, the region has four (4) national parks, the Mount Pulag National Park (MPNP), Balbalasang-Balbalan National Park (BBNP), Mount Data National Park (MDNP), and Cassamata Hill National Park (CHNP), the habitat of many vulnerable and endangered species.

Dependent variable
The primary data used in the spatio-temporal analysis of forest loss in the study area was the HGFC version 1.7. This is a high-resolution forest cover change global map of the twenty-first century showing forest loss and forest gain datasets including the extent of tree cover in the year 2000 [27]. Tree cover data is defined as vegetation with a height of more than 5 m expressed as percent for each cell. Since the focused of the present paper is on the forest loss, the tree cover and forest gain datasets were not used. The main dataset used was the annual forest cover loss (lossyear) during the period 2001 to 2019 with a spatial resolution of 30 m per pixel or one arc-second. Forest cover loss is defined as a change from a forest to a non-forest state, from 2000-2019, coded as either 1 (loss) or 0 (no loss) [31]. These datasets were the results of time-series remote sensing analysis of various Landsat 7 Enhanced Thematic Mapper Plus (ETM?) and Landsat 8 Operational Land Imager (OLI) scenes performed using the Google Earth Engine (GEE), a cloud platform for earth observation. The global forest loss were characterized based on the calculated time-series spectral metrics from the per-pixel set of cloud-free image observations. The improved version can be visualized and downloaded for free from the earth engine partner's website. A worldwide set of 10 9 10 degree granules spanning the range 180 W-180 E and 80 N-60 S can be visualized and downloaded from the website. For the study area, the downloaded dataset was the lossyear with 20 N latitude and 120 E longitude granules.

Independent variables
The LULC data was derived from the National Mapping and Resources Information Authority (NAMRIA) 2015 Vegetative Land Cover Map (VLCM). This was requested from the DENR regional office in CAR. This dataset is also available on the NAMRIA geoportal website. The spatially Fig. 1 Location map of the study area explicit data extracted from 2015 VLCM and subset for each province and municipality were the forest cover, agricultural areas (both annual and perennial crops), and built-up areas (Fig. 2). The road network data was downloaded from the Open Street Map (OSM). It was then clipped to each municipality and its length was calculated in km. The population growth rate data were computed for each municipality using the census of population in the years 2010 and 2015 from the Philippine Statistics Authority (PSA), formerly the National Statistics Office (NSO). The projected 2020 populations were derived as a result of the 2015 census conducted by PSA in collaboration with the Interagency Committee on Population and Housing Statistics (ICPHS). The following formula was used in computing the population growth rate (PGR) in the study area: where PGR is the growth rate, NP is the new population or the 2020 data, the OP is the original population. The average income over the last 10 years per municipality was derived from PSA as well. Figure 3 shows the summary of PGR, the population in the year 2020, and the average income per province. The basis for selecting the independent variables were the suggestions and findings of various studies on the primary drivers of deforestation, either on a national or global scale. It was always observed that there is an inverse relationship between the density of forests and the depletion of forest [32]. Meanwhile, agricultural expansion remains to be the leading factor of deforestation [33,34]. Expansion of human settlements or built-up areas also leads to forest loss in highly-populated areas of the world [35]. Road density also has an impact on forest fragmentation and deforestation [15,36]. Apart from agricultural expansion identified as the forerunner of deforestation, population growth is also identified as one of the leading drivers based on numerous papers [33,35,37]. Lastly, some studies showed a positive relationship between deforestation and income level [38,39].

Geo-processing
The downloaded raster format of HGFC v1.7 lossyear data was used in quantifying the annual forest loss of the study area. The processes involved several geospatial and statistical techniques in a Geographic Information System (GIS) software, particularly the Quantum GIS (QGIS) 3. 16 Hannover version. The provincial and municipal politicaladministrative boundary of CAR in vector shapefile format, requested from the Department of Environment and Natural Resources (DENR) Cordillera Regional Office, were used in clipping the region-wide forest loss data (Fig. 4). After clipping, conversion from raster to vector (polygonize) data was performed. The conversion from raster file format to vector format was necessary for area computations in the attribute table of the map data. After generating the vector data, it was dissolved to the study area. After which, the areas in hectares (ha) of the forest loss per year (from 2001 to 2019) for each municipality were calculated. The data was split to get the individual forest loss data of each municipal boundary feature. After splitting each feature, the count (number) or patches (polygon), maximum, minimum, mean, standard deviation, and total area of forest loss in the study sites were determined. The number of polygons served as an indicator on the patchiness or fragmentation of forest loss in the study site. Finally, the attribute feature of the data along with all the necessary information was exported as a spreadsheet file for further statistical analysis.

Correlation analysis
The present study evaluated the correlation of the forest loss and seven independent variables, which includes forest density, agricultural density, built-up density (units expressed in thousand ha), road network density (in hundred km) from OSM, population growth rate (percent), population data of 2020 (in thousands), and the average annual income (in millions) for ten years per municipality using Pearson's Correlation Coefficient (r) to measure the strength of the statistical association between these variables. A simple linear regression between the forest loss and each of the explanatory variables was performed. The formula for the r is presented below: The evaluation of the strength of relationships between forest loss and the different explanatory or independent variables, expressed in r, was conducted both on a provincial and regional scale to aid in the identification of potential drivers of deforestation in the region as a whole and in the six provinces.

Validation and accuracy assessment
To validate the accuracy of the HGFC data, visual observation using high-resolution historical satellite imageries in Google Earth Pro (GEP) was conducted. The forest loss polygons for validation were randomly selected for each province and overlaid to the GEP. Every province has 100 validation sampling points with a total of 600 validation points for the entire study area (Fig. 4). The sampling design considered was a binomial distribution with either disturbed or undisturbed data. The disturbed data are polygons with a non-forest state caused by either anthropogenic factors or natural hazards, whereas the undisturbed data are polygons with an existing forest cover as observed from the GEP. Each sampling polygon was visually observed on-screen through historical images within the GEP. This process was able to determine whether there were any disturbances (from forest to non-forest state) that occurred in the study area such as agricultural expansion, road and built-up establishment, landslides, and forest fires. Numerous authors conducted data validation using high-resolution images of GEP and produced promising results [31,40].
After completing the validation, a confusion matrix was created and the positive predictive value (PPV) or the precision rate was computed calculated using the Bayes' theorem of Positive Predictive Value (PPV) equation [41] as follows: where PPV is the Positive Predictive Value, TP is the True Positive or the forest loss data with observed disturbances, while the FP is the False Positive or the polygons with undisturbed forest area. Precision measures the portion of actual or true positives among those samples that are predicted as positive with a range value from 0 to 1. The larger the value the better predictive accuracy. Also, this paper evaluated the PPV of the dependent variable based on the different polygon sizes as follows: \ 0.2 ha, 0.2-0.5 ha, 0.5-1 ha, [ 1 ha. This was carried out to determine the applicability of the datasets as a whole in reporting forest Most studies only consider areas larger or equal to 0.5 ha [42], thus, they failed to recognize the small-scale forest disturbances. In this study, however, the considered minimum value for a forest disturbance is 0.2 ha or lesser. Smaller patches of forest loss may seem negligible compared to the larger areas of deforestation, however, the cumulative loss of these small pockets will add up and in turn will have a greater negative impact on the ecosystem as a whole.
3 Results and discussions 3

.1 Spatial and temporal forest loss in CAR
The spatial and temporal analysis of forest loss using the HGFC v1.7 in the study area considered two attributes or fields; the polygon counts and the area in ha. The polygon counts determined the quantity or frequency of forest loss occurring in the site. The results of geo-processing showed that 204,861 polygon counts of forest loss were detected in CAR. The highest was recorded in Apayao province with 96,440 polygon counts, whereas the lowest was recorded in Benguet province with 14,833 patches. In terms of hectarage, the region recorded an accumulated area of 70,925 ha from 2001 to 2019. The province with the highest forest loss over the last two decades was Apayao, with 40,794 ha or 57% of the total loss in the region. This was followed by Kalinga, Ifugao, Mountain Province, Abra, and Benguet with 9467.55 ha (13%), 7126.47 ha (10%), 5770.37 ha (8%), 4574.61 ha (7%), and 3415.85 ha (5%), respectively. Based on the computation of the annual rate of forest loss (Table 1), Apayao province also recorded the highest rate, which is around 2,000 ha/year, while the rest of the provinces had below 500 annual rates of depletion with the lowest being recorded in the province of Benguet (179.78 ha/year). The provinces of Abra, Apayao, Ifugao, and Kalinga showed an increasing trend line or an upward direction of annual loss, while Mountain Province, to a certain extent, showed a sideways or horizontal trend. On a positive note, Benguet province showed a downward trend that keeps on decreasing over time (Fig. 5).
The forest loss statistics presented in the Global Forest Watch (GFW) are found to be inconsistent with the output of this paper. The total forest loss in the region as reported in GFW is around 65,900 ha whereas the total loss computed in this paper is around 70,900 ha. The most probable reason for this is because the political-administrative (provincial and regional) boundaries used in the dashboard of GFW are not consistent with the Cadastral administrative boundaries used by the government agencies in the region such as the DENR, LGUs, among others. Aside from the political boundary inconsistencies, the data from GFW cannot be further processed using geospatial analysis in GIS including statistical and correlation analysis per municipality, accuracy assessment, and data validation. Thus, great caution should be exercised when using the statistical data from GFW in reporting deforestation on a provincial or regional level because it could either overestimate or underestimate the data.
The minimum area or the smallest polygon of forest loss is computed at 0.073 ha or equivalent to 730 m 2 . This is smaller compared with the equivalent area of 30 m pixel (0.09 ha) because of the split process conducted in study. Some 30 m pixel polygons were cut following the boundary of each municipality. Thus, the smallest polygon area was computed at 0.073 ha. On the other hand, the maximum area (ha) detected by the HGFC data in the region is 35.68 ha located in the province of Apayao ( Table 1). The region-wide average mean of forest loss area is 0.30 with a standard deviation of 0.56. This paper suggests that around 60% of the polygons produced by HGFC in the study area had below 0.2 ha, on the other hand, the polygons with more than 1 ha comprise around  (Fig. 6). These findings implied that the majority of the detected deforestation in CAR are characterized by smaller and scattered patches of clearings, but the sum of its area constitutes only 16% of the overall deforested area. These small and less intensive deforested areas may be the result of slash-and-burn practices (locally known as kaingin) as reported by Araza et al. [20]. In line with this, the region's most pressing concern are the polygons with more than 0.5 ha which accounted for 65% of the total depleted area. Thus, this study recommends that rehabilitation efforts or forest landscape restorations should be prioritized in these polygons because, aside from landslides/soil-erosion and habitat loss, extreme warming usually originates in large patches of forest loss [42]. The region-wide temporal monitoring and analysis of forest loss revealed that the year 2014 and 2018 had the highest polygon count of loss with around 16,000. In contrast, the highest loss in terms of area (ha) was recorded from 2016 to 2018 consecutively with 7000 to 8000 ha, but by the end of 2019, a sudden decrease was recorded (Fig. 7). It was also observed on the temporal data per province, except for Apayao, that in the year 2017, forest loss started reflecting a downward trend over the succeeding years (Fig. 5), which is a good indicator as far as forest conservation is concerned. These findings may infer that the priority programs of the DENR which started in the year 2017 are effective in most of the provinces except in Apayao. This program includes an intensified forest protection and anti-illegal logging that aims to neutralize illegal logging hotspots; forest fire prevention; and continuous forest patrolling with the use of an advanced tracking technology which is the LAWIN system application that can be installed in smart phones or tablets [20,43]. The implementation of the Enhanced NGP that started in 2015 could also be seen as a possible factor in the decrease of forest loss during this period. However, an empirical investigation needs to be conducted to further support this inference. Also, it was observed that the highest peak of the trend line in two provinces (Benguet and Mt. Province) was recorded in the year 2009. It is worth mentioning that in the year 2009, two strong typhoons (Pepeng and Ondoy) struck the Northern Luzon Island, particularly in these two provinces, with heavy rainfall that caused a large number of landslides [44]. In this paper, it is plausible to consider that the HGFC data were able to determine the extent of forest loss caused by a natural hazard such as typhoon. Moreover, it has the potential to monitor the impact of government programs and efforts in forest conservation and protection.
From the results, it is clear that the province of Apayao has the most serious threat from deforestation among all the other provinces of the region since it had the highest annual rate of deforestation, as well as the highest area depleted. Generally, the present findings confirm that the study sites had undergone remarkable pressure from forest loss in the last two decades [45]. The major drivers of forest cover disturbances in one of the national park in the region were attributed to the expansion of agricultural areas with 48% of the total causes of forest disturbances; timber poaching or illegal logging (26%), disturbances caused by natural calamities such as landslides (12%); slash-and-burn Fig. 6 The percent distribution of the HGFC loss year data per province based on different area ranges. a Polygon count, b Polygon area (12%); and forest fires (2%) [46]. This may not be true in other parts of the region since site-specific empirical studies on the causative factors of forest loss are non-existent. However, a similar pattern of results was identified in other literature [16] with regard to the key drivers of deforestation in the country. This includes legal/illegal logging and poaching; agricultural expansion through ''kaingin'' making, conversion of forests and grazing; and biophysical factors such as landslides, floods, and typhoons. Overall, the result of the spatio-temporal analysis of deforestation performed in this paper may play an important role in the monitoring and evaluation of the UN's SDG 15 implementation in the region, which can be replicated in other parts of the country.

Correlation analysis
Based on the region-wide correlation analysis as presented on scatter plots (Fig. 8), only the forest cover density has a strong relationship (r-value = 0.78), all other variables were not correlated with the forest loss. The result of other empirical studies highlighted the role of agriculture as a potential driver of deforestation [47,48] but this may not always be the case on a regional scale analysis as shown in this study. There is, however, a strong association between the loss of forest and agricultural density recorded in some of the provinces such as in Mt. Province, Kalinga, and Apayao (with r values of 0.92, 0.76, and -0.59, respectively). Thus, it is more appropriate to do the correlation analysis on a provincial level or even at the municipal level and not on the national or regional level to fully understand the potential drivers of deforestation.
Compared with the regional analysis which gave a strong linear relationship between the dependent variable and the forest cover density, the output of the provincial analysis, as shown in Table 2, suggests that the presence of forest does not always correlate with deforestation as in the case of the provinces of Ifugao (r = 0.47), Kalinga (r = 0.41), and Mt. Province (r = 0.44). This is contrary to the common notion that there is always an inverse 123 relationship between forest density and deforestation [32]. The probable reason for this might be attributed to the strong indigenous forest management systems and practices in these three provinces such as the Muyung system in Ifugao [49]. It could also be attributed to the conservation of natural forests through the Community-Based Forest Management (CBFM) program of the government since the region has 84 CBFM areas with a corresponding area of 61,949 ha [18].
For the built-up density and road network density variables, only two out of the six provinces in the study area recorded a positive relationship. These are the provinces of Kalinga (r values = 0.78 and 0.74) and Mt. Province (r values = 0.58 and 0.85). This result indicates that forest loss in the region does not appear to have been affected by built-up or road network variables but rather it is influenced by a more complex pattern of anthropogenic encroachment [36]. This is contrary to the popular belief that the closer a forest was to roads, the higher the rate of deforestation [15]. The forest loss is not correlated with road network because the country is already in a post-logging era which means that logging roads are nonexistent. In the case of CAR, there are no longer logging concessions. Further, the no correlation r-values may also be the result of the implementation of a total log ban on natural and residual forest in the country through the Executive Order (EO) No. 23 series of 2011.
The income per municipality could also be a potential driver of forest loss as observed in the provinces of Apayao, Kalinga, and Mt. Province. It is interesting to note that these three provinces had the lowest income among the other provinces. This conveys the idea that the lower the income the higher rate of deforestation. Other studies claimed that when there is a low human development index (HDI), which includes the income, with a high population growth rate, there was a high rate of deforestation, but when this explanatory variable is high, the rate of deforestation was low, regardless of the rate of population growth [37]. In this paper, however, it was proven that forest loss in a particular province can be linked to its low income irrespective of the population growth rate.
In terms of correlation coefficient analysis per province, Kalinga and Ifugao exhibited a moderately strong relationship between the dependent and independent variables (population growth rate) with r values of 0.54 and 0.52, respectively. Furthermore, the province of Apayao showed a strong (-) association between forest loss and population growth rate (r = -0.75) while the provinces of Abra and Mt. Province both showed a weak correlation. It is worth noting that only Benguet has no correlation (r = 0.11) among the two variables despite the fact that this province had the highest recorded population growth rate (0.17) and the most populous province with 46% of the total population in the region. Moreover, this province is where the City of Baguio (Hill Station) is geographically located, an urban area and the economic powerhouse of the region where the central business district, regional government offices, popular tourist attractions, and major universities in the region are located [50]. Contrariwise, the province of Apayao, with the least population (only 7% of the total population in the region) and a negative population growth rate (-0.03), a province devoted to agricultural production, furniture industry, crafts, and housewares making [51], recorded the highest forest loss rate. These findings indicate that population growth in an urban area is not always associated with the reduction of forest cover on a provincial scale as compared to rural areas where most of the population depends on upland agriculture or slash-andburn for subsistence. As Yameogo [52] puts it, population density, urbanization, and economic growth have a significant and positive impact on forest loss. Also, Liu et al. [53] claimed that the relationship between urbanization and vegetation/forest degradation appears to be nonlinear, which means that urbanization will not necessarily result in forest and/or vegetation degradation. In relation to this, Basnyat [54] also claimed that the relationship between population dynamics and forest degradation is not linear. It is, in a sense, a logical conclusion that the migration of peoples from rural to urban areas will decelerate the loss rate in a particular province. Counter-intuitively, this is not always the case as shown in other studies [33,47] that forest loss is positively correlated with two major driving forces: population growth and agricultural expansion. Either way, this study highlights the fact that population growth in a rural or an urban area, does not always associate with forest loss in the study site.

Validation and accuracy assessment
The validation of the accuracy of the HGFC data based on randomly selected points for each province revealed that only 54% of the validated polygons (600 sampling points) regardless of their area, had an accurate detection of forest loss in the study area. The low validation accuracy could be attributed to the geolocation error of smaller polygons (30 m) in the reference points such as the GEP. However, a higher percent accuracy was observed based on the equivalent area of the validated polygons (1928.32 ha) with an accuracy rate of 72% (Table 3). The accuracy assessment of the HGFC data based on different polygon sizes pointed out the highest and lowest accuracy relative to the area. In this case, polygons with more than 1 ha demonstrated a higher accuracy (with 66%) compared with the polygons with smaller areas. It is also important to note that HGFC polygons with less than 0.2 ha gave the lowest accuracy with only 34% ( Table 4). The validation conducted in the study strongly suggests that great caution should be observed when using the forest loss data with polygons lower than 0.2 ha, particularly in reporting deforestation rates as well as in planning management interventions.
Overall, the HGFC data demonstrated its applicability in spatial and temporal monitoring of the annual forest loss in the study area based on the validation which produced a 72% accuracy rate. In the Philippines, the HGFC data can be used with moderate to high reliability as suggested by Fallacurna and Perez [31], but the forest loss in an actual scenario may even be worst due to underestimation [19]. Moreover, its applicability was also observed in the other parts of the world as shown in the studies of Austin et al. [28] in Indonesia, and other countries of Southeast Asia [29]. However, some concerns were raised by other authors in terms of its accuracy. According to Tropec et al. [55], the HGFC datasets failed to distinguish the change in vegetation from tropical forests to herbaceous crops such as pineapple, soybeans, and tea plantations, which resulted in an underestimation of forest loss and compromises its reliability for policy decisions. This claim was further supported by Curtis et al. [40], stating that the data does not distinguish permanent forest loss from temporary disturbances that include shifting cultivation, forest fires, and forestry plantations. Also, there was an underestimation of forest loss in countries with dry forests like in Ghana [56,57]. On the contrary, Tracewski et al. [58] claimed that the HGFC represents the best high-resolution forest cover change data available worldwide. A promising result was also reported in detecting forest loss in the Protected Areas (PA) worldwide using the same data [59]. Based on the study of Vieilledent et al. [60] in Madagascar, they combined the HGFC data and the historical national forestcover maps to monitor deforestation and forest fragmentation over the last 60 years, which proved its applicability in monitoring forest loss. Furthermore, a high overall accuracy (0.88) was derived from the accuracy assessment conducted by Bos et al. [61] in four tropical countries (Indonesia, Peru, Tanzania, and Vietnam). This is congruent with what has been found by Hansen et al. [27], which recorded a 13% false positive, or 87% true positive (accuracy) globally. Nonetheless, due to the consistent methodology in generating the data, the HGFC is so far the best remotely sensed data, with a high-resolution (30 m), applicable for forest cover change monitoring on a local or a global scale analysis that is freely accessible [62]. This paper also highlights the applicability of GEP in the validation of the HGFC data owing to its high-resolution images and historical data. As shown in Fig. 9, one of the validation samples in the province of Kalinga reveals the status or condition of its forest area across different periods. The sampled forest loss polygon was detected in 2014 by HGFC data which is very consistent with the GEP on the year the forest was depleted.

Conclusions
The validation results of the forest loss detected by the HGFC lossyear data are a reminder that immediate, science-based, and sustained regional and multi-stakeholder efforts are necessary to conserve and protect the remaining forest cover of one of the last forest frontiers in the country. This is to minimize and/or eradicate the negative impacts of forest loss, both in the upstream and downstream communities, such as water shortage, river siltation, flooding, landslides, erosions, and biodiversity loss, thereby achieving the SDG 15 (Life on Land) target. Furthermore, the result of this paper could be utilized by concerned government and non-government organizations to efficiently monitor the forest loss in the region, as well as to plan out present and future interventions through the local and provincial-level development plans such as the Comprehensive Land Use Plan (CLUP), Forest Land Use Plan (FLUP), Integrated Watershed Management Plan (IWMP), among others. The HGFC data could potentially determine the extent of forest loss caused by a natural disaster such as typhoon. Also, it can be useful in monitoring the impact of government programs and efforts in forest conservation and protection. However, great caution should be observed when using the polygons with below 0.2 ha due to their very low accuracy. The Radar for Detecting Deforestation (RADD) and GLAD-S2, developed by the GFW, also offer deforestation alert systems, which can be used as an alternative to the HGFC dataset. A deforestation hotspot analysis could also be conducted in the region to identify priority areas for conservation. Correlation analysis of provincial or municipal level will likely produce a more reliable result compared with the national or regional scale analysis. Finally, the identified potential drivers of deforestation in this paper should be considered by various stakeholders to ensure a holistic implementation of forest landscape restoration as well as to facilitate decisionmaking processes.
completed without the full support of Ms. Sarah Jane and Mr. Paul Isaac.
Funding No funding was received for conducting this study.

Declarations
Conflict of interest The authors have conflicts of interest to declare that are relevant to the content of this article.
Data availability The Hansen Global Forest Change version 1.7 datasets generated during and/or analysed during the current study are available in the earth engine partner's website repository http://earth enginepartners.appspot.com/science-2013-global-forest. The software QGIS version 3.16 Hannover used in the geospatial analysis of this study can be downloaded from the QGIS website https://www.qgis. org/en/site/. The census of population in the Philippines, including the projected populations and income level, used in this study can be retrieved from the Philippine Statistics Authority (PSA) website https://psa.gov.ph/statistics/census/projected-population. All data generated or analysed during this study are included in this published article.