Protected areas slow down tropical rainforest disturbance in the Leuser Ecosystem, Indonesia

ABSTRACT

Tropical rainforest ecosystems that function as biodiversity pools had been undermined because of anthropogenic activities. Research has shown that protected areas (PAs) have become the first safeguard for biodiversity. However, how to measure the effectiveness of PAs remains unclear. We present spatiotemporal changes within the PAs and non-PAs in the Leuser Ecosystem, which is one of the significant global landscapes, using intensity analysis during two time periods and propensity score matching to investigate the effectiveness of PAs. We classified land cover using machine learning based on remotely sensed data. Our results revealed the effectiveness of PAs compared with non-PAs. The new conservation intervention after 2008 resulted in the deacceleration of deforestation from 2000–2010 to 2010–2020. In addition, PAs can reduce deforestation two times more effectively than non-PAs. Therefore, PAs and good governance within the Leuser Ecosystem are crucial in maintaining the natural ecosystem to address global conservation targets.

KEYWORDS:

• Intensity analysis
• land change
• Leuser Ecosystem
• protected areas
• statistical matching

Introduction

As a part of the Tropical Rainforest Heritage of Sumatra (Claudino-Sales, 2019), the Leuser Ecosystem significantly contributes to global conservation through its ecosystem services and biodiversity habitats (Le Saout et al., 2013). It is considered to be the last natural habitat for certain endangered mammals such as orangutans, elephants, tigers, and rhinoceros (Cochard, 2017). The Leuser Ecosystem notably contributes to carbon sequestration (Warren et al., 2017), provision of water (Cochard, 2017), and provides various products for community livelihoods in the surrounding landscape (Griffin et al., 2013). The Leuser Ecosystem falls under two general designations for conservation or economic developments through forest production and road infrastructure (Sloan et al., 2019, 2018). This landscape also functions as a National Strategic Area to preserve ecosystem services that restrict further road developments. However, a previous study highlighted some redundancies in the legal act related to the landscape, implying a shortfall in the conservation plan and strategies of the Leuser Ecosystem (Sloan et al., 2019).

Indonesia is rapidly expanding through infrastructure and accessibility development (Laurance et al., 2009; Sloan et al., 2018), across its major islands (e.g. Sumatra). This rapid expansion has caused a natural ecosystem disturbance. Recent reports indicate have indicated that anthropogenic activities have significantly changed the land cover (Margono et al., 2012). Agricultural incursions and logging activities are reportedly the significant drivers of deforestation and forest degradation in Sumatra (Margono et al., 2014, 2012), particularly in the Leuser Ecosystem (Gaveau, Wich, et al., 2009). Forest degradations are responsible for high carbon emissions that lead to climate change (Page et al., 2011). In addition, deforestation results in habitat losses for its threatened biodiversity (Gaveau et al., 2009; Kinnaird et al., 2003; Linkie et al., 2008; Weiskopf et al., 2019; Wich et al., 2016). Previous studies showed that protected areas could effectively function as habitat refugia and biodiversity hotspots (Condro et al., 2021a; Haight & Hammill, 2020). Thus, protected areas should be known as part of a comprehensive approach for sustainable management that can overcome the cumulative effects of anthropogenic activities such as biodiversity loss (Condro et al., 2021b; Magioli et al., 2021). It is essential to understand the role of the protected areas in the conservation process; however, elucidating the performance of protected areas has proved challenging (Geldmann et al., 2020).

Land use or land cover change (LULCC) analysis helps to evaluate the trade-offs between land categories to understand the coupled human–environmental mechanism (Lambin & Meyfroidt, 2011; Mendoza-Ponce et al., 2018). The spatiotemporal dynamics can represent the underlying anthropogenic activities that interact with the natural landscape over time (Song et al., 2018). The massive landscape alterations result in loss of biodiversity, climate disturbance, decreased land productivity, and tenurial conflicts in the community, affecting long-term environmental sustainability (Duveiller et al., 2020; Yohannes et al., 2020; Zhao et al., 2018). However, land transformation is investigated using data and techniques, and errors in LULC data can deliver through change analysis. Hence, robust classification methods and land cover change analyses are required to provide an in-depth understanding of the changes in the landscape over time (Xie et al., 2020).

In recent years, several studies have assessed the spatiotemporal patterns of land changes (Chughtai et al., 2021; Naboureh et al., 2021). Most of these studies compared two-time land cover maps through a transition matrix to measure temporal change within a set of categories (Quan et al., 2019; Rogan & Chen, 2004). Some studies have also used a transition matrix to evaluate land cover changes over time (Hermon, 2016; Koh et al., 2011). However, studies have also shown the shortfalls of using this matrix in providing the underlying land cover dynamics and revealing gross loss/gain for each category (Huang et al., 2018; Pontius, 2019). Therefore, intensity analysis was introduced to investigate comprehensive land cover changes.

The intensity analysis is a framework to provide various signals of change through interval, category, and transition levels for land cover change investigation (Aldwaik & Pontius, 2012; Varga et al., 2019). This framework can better understand land cover change analysis than the transition matrix (Aldwaik & Pontius Jr., 2013). Nyamekye et al. (2020) assessed urban growth in Ghana using intensity analysis. Another study applied this framework to examine the land change dynamics in the African-emerging city (Akinyemi et al., 2017). Quan et al. (2018) implemented intensity analysis to communicate land cover dynamics in two regions of China during three time points. Xie et al. (2020) enhanced the visualization of the transition pattern in the intensity analysis framework. Other researchers used intensity analysis for evaluating land dynamics among categories (Cribari et al., 2021; Ekumah et al., 2020; Varga et al., 2019).

Measuring the effectiveness of protected areas is critical to determining future conservation strategies (Nelson & Chomitz, 2011; Reboredo Segovia et al., 2020). The most popular way to evaluate the performance of conservation policies within the protected areas is to estimate the deforestation rate in the landscape through land cover change analysis (Brun et al., 2015; Cuenca et al., 2016; Gu et al., 2020). However, the comparison of avoided deforestation between protected and non-protected areas would likely be flawed because it neglects the lack of randomness in the protected areas (Blackman et al., 2015; Graham et al., 2021). Recent studies have used statistical matching [e.g. propensity score matching (PSM) analysis] to account for the nonrandom location of protected areas within the landscape (Gaveau et al., 2009; Geldmann et al., 2020; Gu et al., 2020).

In this study, we present a comprehensive intensity analysis of the Leuser Ecosystem, Indonesia, which is one of the most significant landscapes in the world for biodiversity and ecosystem services (Le Saout et al., 2013), for two periods of time within two regions (i.e. protected and non-protected areas). We also evaluate the effectiveness of the protected areas within the landscape using PSM.

Methods

Study area

The Leuser Ecosystem is located on North Sumatra Island, Indonesia. It covers 26,495.89 km² and consists of protected areas (i.e. nature reserve areas, hunting parks, and nature conservation areas) and non-protected areas (i.e. production forests, limited production forests, and convertible forests) based on its forest management units. The size of the non-protected areas (16,810.62 km²) is almost two times larger than that of the protected areas (9,685.28 km²).

Land cover data

In this study, we captured deforestation dynamics and its drivers using six land cover categories defined by the Intergovernmental Panel on Climate Change (IPCC; Hirschmugl et al., 2018), based on machine learning classification of Landsat imageries for 2000, 2010, and 2020 (Gorelick et al., 2017). For our study, we modified the IPCC’s land cover classes for our study areas into the following categories: Forest, Wetland, Agriculture, Settlement, Shrubs, and Waterbodies. In this context, Settlement consists of built-up areas, roads, and other infrastructure objects that cover non-vegetated areas.

This study performed machine learning classification based on the Random Forest algorithm using cloud computing (Breiman, 2001). We used 25 variables derived from Landsat imageries (i.e. reflectance bands and spectral indices; Table S1). We performed visual interpretation to obtain training samples of the land cover classes within the study area (see Supplementary Materials for further details on land cover identification). Furthermore, we stratified the Leuser Ecosystem into protected and non-protected areas to capture the discrepancies in land dynamics within and outside protected areas.

Error analysis

We quantified the confusion matrix in 2000 and 2020 between the raster maps and sample references of land cover in the study area. We then performed intensity analysis to measure the error at the category level in terms of size and intensity (Pontius, 2019). In this analysis, the misses correspond to the omission error, false alarms correspond to commission error, and hits correspond to the agreement. We also calculated the component error of land cover maps in 2000 and 2020, i.e. Quantity, Exchange, and Shift (Pontius & Santacruz, 2014).

Intensity analysis

This study allocated the change during each period of time into three components: Quantity, Exchange, and Shift (Pontius & Santacruz, 2014). The Quantity component measures the absolute net change in the size of the classes. The Exchange component measures the transition of category i to category j in some areas or the reverse transition of category j to category i in other areas within a region during a period of time. The Shift component for category i occurs when both C_rtij < C_rtji and C_rtik > C_rtki for at least one combination of categories j and k. All equations used in this section are available in previous studies (Pontius, 2019; Quan et al., 2019).

Intensity analysis is a mathematical framework toward patterns of change in a hierarchical structure that compares uniform intensities to observed intensities of temporal difference among categories (Aldwaik & Pontius Jr., 2013). Intensity analysis includes interval, category, and transition levels. The interval level compares time intervals in terms of gross change. The category level measures the gains and losses of categories during each period of time in each region. The transition level compares the transitions in terms of intensity in each region during each period of time. We used transition pattern matrix introduced by Xie et al. (2020) to communicate the transition level. In addition, we visualized the land cover changes within protected and non-protected areas through the Sankey diagram (Muenchow et al., 2019) using ggplot2 package in R.

To analyse land change, we adopted the framework containing (1) data pre-processing, (2) land cover identification, (3) error analysis, and (4) change analysis from Quan et al. (2019). Figure 1 shows the conceptual framework for intensity analysis.

Figure 1. Conceptual framework for analyzing land change.

Effectiveness of the protected areas

The effectiveness of protected areas from the anthropogenic disturbances can be calculated by comparing the variables inside and outside protected areas (Geldmann et al., 2020). Slope and elevation, distance to forest, distance to roads, and distance to agricultural areas are the critical predictors that best represent human pressures within protected areas in Sumatran rainforests (Gaveau et al., 2009). However, statistical matching (Schleicher et al., 2020) for comparing protected and non-protected areas (e.g. propensity matching analysis) should be evaluated because the bias toward remote areas were having a high similarity to the non-randomly located protected areas (Joppa & Pfaff, 2009).

The optimal full matching method from the MatchIt package in R was used to perform PSM (Ho et al., 2011) based on the Random Forest algorithm (Breiman, 2001) with 500 number of trees. The binary values were created from protected (n = 2256) and non-protected areas (n = 12,080) for the treatment variables in PSM analysis. We performed exploratory analysis to assess statistical matching quality using love plots of standardized mean differences and distributional balance of distance between before and after matching (Oldekop et al., 2019). We also used Moran’s I and spatial distribution of post-matching residuals to investigate spatial autocorrelation and check for hidden bias within the model (Schleicher et al., 2017).

In this study, protected area boundaries were considered as dependent variables and were retrieved from the World Database of Protected Areas (WDPA; UNEP-WCMC and IUCN, 2020). Moreover, based on the previous studies regarding protected areas’ effectiveness (Gaveau et al., 2009; Geldmann et al., 2020), we used slope and elevation (Jarvis et al., 2006), distance to forest edge, distance to road, forest cover in 2000, climatological temperature from WorldClim v.2.0 (Fick & Hijmans, 2017), and human modification index (Kennedy et al., 2019) to capture the confounding variables. The present studies excluded buffers of 10-km around the protected area boundaries to deal with the neighborhood leakage effect (Cuenca et al., 2016; Gaveau et al., 2009).

This study used percentage deforestation as the response variable with two periods (2000–2010 and 2010–2020), following Gaveau et al. (2009). Then, we compared deforestation estimated from matching analysis within the protected and non-protected areas (Cuenca et al., 2016). We defined ‘deforestation’ as the alteration of Forest and Wetland categories into the other land cover types (i.e. non-Forest and non-Wetland) with the values ranging from 0% to 100% (Laurance et al., 2009).

Results

Land cover data and error analysis

Figures 2 and 3 show the maps and areas of the six land categories at the three time-points having the same interval of years (i.e. 2000, 2010, 2020). Within the Leuser Ecosystem, the Forest category dominated the protected areas from 2000 to 2020. Shrubs was found to be the largest category under non-protected areas at all time points (Figure 3). Figure 2 presents the net changes in land cover for each category. From 2000 to 2020, Forest experienced a net loss of 19% of the protected areas, while Shrubs experienced a net gain of 12% of the protected areas. With respect to non-protected areas, Forest experienced a net loss of 18%, while Agriculture experienced a net gain of 16%.

Figure 2. Maps of the Leuser Ecosystem, Indonesia over three time points within protected and non-protected areas.

Figure 3. Sankey diagram of each land category in the Leuser Ecosystem within protected and non-protected areas at three-time points.

Figure 4 describes the category level errors for (A) size and (B) intensity of the Leuser Ecosystem in 2000 and 2020. In Figure 4(a), the union of hits and misses indicates the number of observations in the reference information for each category. For a particular category, if the false alarms were more extensive than the misses, it indicates that the map overestimated the size of that category. In the case of smaller false alarms than the misses, the map underestimated the size of that category. In 2020, the size of the Forest, Agriculture, and Settlement categories was insignificantly overestimated compared with the size of the rest of the categories, which was insignificantly underestimated. The underestimation of Shrubs resulted in the most significant component of error in 2020. In 2000, the size of the Forest, Wetland, and Agriculture categories was insignificantly overestimated compared with the size of the other categories, which was insignificantly underestimated. In addition, the underestimation of Shrubs and the overestimation of Forest resulted in the most significant error in 2000. In Figure 4(b), the uniform dashed lines represent that the overall error in 2020 and 2000, which was 3% and 5% of the Leuser Ecosystem extent, respectively.

Figure 4. Summary of category level of errors in 2000 and 2020 for (a) size and (b) intensity. The dashed lines indicate the overall error as a percentage of category size.

The maps overestimated the size of the Forest, Wetland, and Agriculture and underestimated the size of the Shrubs and Settlement; this resulted in the quantity error in 2000. Similarly, the quantity error in 2020 occurred because of the overestimation of the size of the Forest, Agriculture, and Settlement categories and the underestimation of the size of the Wetland and Shrubs categories by the map. The exchange occurred due to confusion regarding the pairs of categories. No shift component of error was observed in 2000. The exchange component accounted for most of the errors within all categories, except for the Settlement category in 2000. For the period of 2020, the quantity component accounted for most of the errors within all categories, except for the Agriculture and Settlement categories, which had the most significant intensity of shift component error in 2020 (Figure S1).

Change analysis

The results showed an extensive land cover change in the non-protected areas compared with the land cover change in the protected areas (Figure 5). As shown in Figure 6(a), the intensity of interval level change produced a consistent pattern because the protected areas and non-protected areas occupy a relatively similar spatial extent. Quantity was found to be the largest component during the two periods, particularly from 2010 to 2020. Our results also showed that the land cover change decelerated across the two-time intervals. The changing intensity with respect to a temporal extent during 2000–2010 was faster than the uniform speed, while the change during 2010–2020 was slow.

Figure 5. Interval level change components (a) and category level losses and gains (b) in terms of intensity at two time periods within protected and non-protected areas of Leuser Ecosystem.

Figure 6. Transition patterns for two time periods. Each square’s size depicts the transition intensity. The gradient of colors indicates the difference between transition and uniform intensities. Red colour means targeting and blue means avoiding.

During each period within the protected areas, Forest exhibited the most significant losses (ranging from 633.90 km² to 1,218.62 km²), while Shrubs exhibited the most significant gains (ranging from 564.40 km² to 1,202.00 km²). The Settlement category did not show any losses during all periods. In the non-protected areas, Forest showed the most extensive loss (ranging from 391.95 km² to 1,534.87 km²), while Shrubs showed the most significant growth from 2010 to 2020 (416.11 km²). In comparison, Agriculture showed the most significant increase from 2000 to 2010 (1,561.68 km²). Figure 6(b) demonstrates the uniform change intensity as a black line for each region during each period. For both protected and non-protected areas during all periods, Forest showed active losses (i.e. more intensive than uniform), while Settlement and Agriculture showed active gains. Moreover, Wetland decreased actively in the non-protected areas during all periods, while Agriculture decreased actively in the protected areas from 2000 to 2010. Meanwhile, Wetland and Shrubs increased in the non-protected areas during 2010–2020 and the protected areas during 2000–2010. Water experienced dormant gains and losses in both areas during all periods.

Figure 6 represents the area and intensity of transition patterns during the two periods. The largest transitions were found in the protected areas from Agriculture to Shrubs during 2000–2010 and Forest to Shrubs during 2010–2020. A relatively large transition was also observed from Forest to Shrubs within the protected areas from 2000 to 2010. With respect to the non-protected areas, the most significant transitions occurred from Wetland to Agriculture during 2000–2010 and Wetland to Shrubs during 2010–2020. In addition, intensive transitions from Forest to Shrubs were observed within the non-protected areas during 2000–2010. In the non-protected areas, a substantial portion of the gain in Shrubs was derived from Forest during 2000–2010 and Wetland during 2010–2020. Besides, the gain in Agriculture targeted Wetland and Shrubs during all periods. The transition from Wetland to Shrubs was relatively large during the first period. In the protected areas, the gain in Shrubs was derived mainly from Forest during all periods and from Agriculture during 2010–2020. The gain in Settlement was derived from Agriculture during the first period in all regions. This study found a more intensive transition during the 2000–2010 period than the 2010–2020 period for both regions (i.e. protected and non-protected areas). Our results also showed more significant transitions in the non-protected areas than in the protected areas during all periods.

Effectiveness of the protected areas on avoided deforestation

Significant differences in means and statistical distributions of each control variable that represents anthropogenic and physical drivers were observed between treatment (protected) and control (non-protected) groups. This is because all the p-values from the t-test from matching groups were found to be highly significant (Table S2). Insignificant p-values toward control variables in the post-matching groups indicated covariate balance after the matching (Table S2). Thus, PSM has adequate balanced measures of anthropogenic disturbance across treatment and confounding variables (see Supplementary Materials for further details of post-matching quality assessment; Figure S2, Figure S3). The discrepancy in deforestation between protected and non-protected areas could be connected to the key predictor after matching.

The mean difference in the percentage deforestation rate between the non-protected vs. protected areas during the 2000–2010 period was highly significant (t = – 6.3701; d.f. = 12,277; p-values < 0.001), and positive (Figure 7(a), Table 1). In addition, this result showed consistency of robustness in means comparison when measured using a Wilcoxon test (z = – 5.83; p-value < 0.001). The mean difference in the percentage deforestation rate between the non-protected vs. protected areas during the 2010–2020 period was also significant (t = – 6.5542; d.f. = 12,087; p-values < 0.001), and positive (Figure 7(b), Table 1). This result showed the consistency of robustness in means comparison as confirmed using a Wilcoxon test (z = – 9.88; p-value < 0.001). The deforestation rate in the protected areas was significantly lower than that in the non-protected areas. In addition, PSM analysis indicated that protected areas can reduce the impacts of deforestation within the landscape.

Figure 7. Comparison of the deforestation rate during the 2000–2010 period (a) and the 2010–2020 period (b) within protected areas and non-protected areas of the Leuser Ecosystem.

Table 1. Comparison of mean differences in percentage deforestation rates per year and 95% confidence intervals across between non-protected areas (non-PA) and protected areas (PA) during 2000–2010 and 2010–2020 periods in the Leuser Ecosystem.

non-PA vs. PA	Deforestation 2000–2010	Deforestation 2010–2020
Before matching
Mean difference (%)	22.55	17.30
95% CI	(19.72–25.38)	(16.14–18.46)
Post-matching
Mean difference (%)	15.23	15.64
95% CI	(13.49–16.97)	(14.29–16.99)

Discussion

This study used intensity analysis to measure the map errors and assess land cover changes in two different land management areas (i.e. protected and non-protected areas) of the Leuser Ecosystem, Indonesia. The study period covered 20 years, comprising two periods of time: 2000–2010 and 2010–2020. Some researchers believed that maps effectively analyzed land cover change when its accuracy (based on the hits) at each time point was greater than 85% (Quan et al., 2015; Xie et al., 2020). However, the metric that depicted the percentage of hits did not reveal to bring out the various types of errors (Pontius et al., 2011). Thus, some investigators commonly used the confusion or error matrix to summarize the various types of errors. They assumed the remarked alterations to be the correct changes in order to explain their results. Based on the Random Forest algorithm, our results revealed a relatively low percentage of error for 2000 (5%) and 2020 (3%) in the Leuser Ecosystem. Therefore, we believe that data error has no significant influence on land changes for both periods of 2000–2010 and 2010–2020. Previous studies showed that the Random Forest algorithm performs well in identifying land cover based on remotely sensed data (Buchhorn et al., 2020; Condro et al., 2020; Oliphant et al., 2019). Moreover, the present study considered the neighbourhood leakage effect across the buffer zone of protected areas and the bias associated with controlling the non-random location of the designated protected areas using statistical matching analysis. Thus, this study could significantly improve assessments of the performance of the protected area within the Leuser Ecosystem.

The intensity analyses at interval levels showed a more significant land cover change in the non-protected areas than in the protected areas. This phenomenon indicated the high cultivation activities (e.g. oil palm and timber plantations) outside the protected areas within the Leuser Ecosystem (Margono et al., 2012). We also observed more intense land cover changes in the first period (2000–2010), which implies the development pressures before the government legally enacted the Leuser Ecosystem as environmental planning priorities based on Law Number 26/2007 and Government Regulation Number 26/2008 (Sloan, 2014). Despite their ecological and social importance, Forest and Wetland in the non-protected areas were found to decrease during all periods. This decrease could be linked to anthropogenic activities due to Settlement and Agriculture expansions. The changes in peat-swamp forests (i.e. Wetland) were driven mainly by oil palm conversion. In comparision, the changes in forests were considerably driven by rubber plantations, mixed agroforestry, etc. (Yarrow Robertson & Van Schaik, 2001). Moreover, Wetland within the protected areas was more stable with minimal alterations than that within the non-protected areas. However, Forest losses were also occurred in the protected areas within the Leuser Ecosystem. In addition, the Leuser Ecosystem is threatened by the legal or illegal development of infrastructure (e.g. roads and encroachments) in the area (Sloan et al., 2018).

Our findings provide information on the magnitude of transitions of the various categories (Varga et al., 2019). In the non-protected areas, significant transitions from Wetland to Agriculture, from Forest to Shrubs, and from Shrubs to Agriculture were observed in the first interval. Besides, Agriculture and Wetland categories avoided Shrubs in the first interval within the non-protected areas. In the non-protected areas, the transition from Forest to Shrubs and Wetland to Agriculture decreased in terms of intensity in the second interval. A significant transition was observed from Agriculture to Shrubs in the first interval within the protected areas, whereas more uniform transitions were observed in the second interval within the protected areas. This result revealed that faster change processes in land cover were occurred during the first interval in both areas than in the second interval. A previous study also showed that natural ecosystem loss (i.e. Forest and Wetland) was mainly driven by Agriculture and Shrubs, particularly in non-protected areas due to cultivation and logging activities (Margono et al., 2012).

Large-scale industrial logging concessions and land clearance for agricultural expansion were responsible for the transition of Forest to Shrubs (Abood et al., 2015; Gaveau et al., 2009). In addition, the category level of intensity analysis captured the road-infrastructure development within the Leuser Ecosystem. This phenomenon reflected the proliferation of anthropogenic activities through natural resource accessibility (Poor et al., 2019; Prasetyo et al., 1995; Sloan et al., 2018). Nevertheless, with the help of good conservation governance and practices, protected areas can safeguard biodiversity and reduce deforestation and ecosystem disturbance in the Leuser Ecosystem. Within this ecosystem, Forest and Wetland play an essential role as biodiversity habitats and provide various ecosystem services (Cochard, 2017; Olson et al., 2001). In agreement, previous works have also shown the resilience of protected areas toward global climate change (Condro et al., 2021a; Struebig et al., 2015). Therefore, rapid natural ecosystem disturbance can be overcome by increasing the network of protected areas and by protecting the current protected areas from degazettement (Le Saout et al., 2013).

A shown in previous studies, the protected areas were generally more effective in avoiding deforestation and human disturbance compared with the unprotected areas (Gaveau et al., 2009; Geldmann et al., 2020; Schleicher et al., 2017). In addition, the present study found that the protected area network suffered two times less deforestation than the non-protected areas during 2000–2020 despite the high levels of anthropogenic activities in the tropical rainforest landscape (Graham et al., 2021; Figure 7). These results also concur with prior studies that revealed that the protected areas avoided deforestation over earlier periods (Gaveau et al., 2009).

Conclusions

This study examines the application of intensity analysis for two different land-use management areas (i.e. protected and non-protected areas) across two periods in the Leuser Ecosystem, Indonesia. This study aimed to evaluate land-use change and interpret the obtained results concerning the conservation policy within the landscape. Besides, we also evaluated the land cover maps through error analysis to avoid inaccuracies of the commonly used metric (e.g. Kappa). We also measured the effectiveness of protected areas within the Leuser Ecosystem using statistical matching analysis to reduce bias toward several effects. Our results show that the new regulation of environmental planning priorities resulted in the deacceleration of land change during the period of 2010–2020 compared with the period of 2000–2010. In addition, we found that the land change was more intensive in the non-protected areas than in the protected areas. Our findings suggest that protected areas within the Leuser Ecosystem play an essential role in reducing deforestation and forest degradation compared with non-protected areas. Forests and Wetlands categories showed decreased activity, while Settlement and Agriculture categories accounted for significant net expansions during all time intervals in both regions. In addition, the gain in Agriculture and Shrubs targeted Forests or Wetlands, indicating that the incursions from cultivation and logging activities were still going on, mainly in the non-protected areas within the forests of the Leuser Ecosystem. It is worth noting that the protected areas are our last hope for reaching optimistic global conservation targets.

Acknowledgments

The authors would like to express gratitude toward The Laboratory of Environmental Analysis and Geospatial Modelling, Faculty of Forestry and Environment, IPB University, for supporting this research by providing computational devices and Irene Carolina Fatmayanti for her administrative assistance.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available on request from the corresponding author, A.A.C. The data are not publicly available due to their containing information that could compromise the privacy of research participants.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/1747423X.2022.2115571

Additional information

Funding

This research was funded by the Ministry of Research and Technology of the Republic of Indonesia through PMDSU Scholarship, grant number 2127/IT3.L1/PN/2021; Kementerian Riset Teknologi Dan Pendidikan Tinggi Republik Indonesia.