Assessing swidden land use in Myanmar by decision tree-based detection method using landsat imagery

Abstract

Swidden agriculture is a common land use found in the mountainous regions, especially in Southeast Asia. In Myanmar, the swidden agriculture has been practicing as an important livelihood strategy of millions of people, mainly by the ethnic groups. However, the extent of swidden agriculture in Myanmar is still in question. Therefore, we attempted to detect swidden patches and estimate the swidden extent in Myanmar using free available Landsat images on Google Earth Engine in combination with a decision tree-based plot detection method. We applied the commonly used indices such as dNBR, RdNBR, and dNDVI, statistically tested their threshold values to select the most appropriate combination of the indices and thresholds for the detection of swidden, and assessed the accuracy of each set of index and thresholds using ground truth data and visual interpretation of sample points outside the test site. The results showed that dNBR together with RdNBR, slope and elevation demonstrated higher accuracy (84.25%) compared to an all-index combination (dNBR, RdNBR, dNDVI, slope, and elevation). Using the best-fit pair, we estimated the extent of swidden at national level. The resulting map showed that the total extent of swidden in Myanmar was about 0.1 million ha in 2016, which is much smaller than other previously reported figures. Also, swidden patches were mostly observed in Shan State, followed by Chin State. In this way, this study primarily estimated the total extent of swidden area in Myanmar at national level and proved that the use of a decision tree-based detection method with appropriate vegetation indices and thresholds is highly applicable to the estimation of swidden extent on a regional basis. Also, as Myanmar is the largest country in mainland Southeast Asia in area with a great majority of the population living in rural areas, and many in the mountains, its land resources are of great relevance to the people’s livelihoods and thereby the nation’s progress. Therefore, this study will contribute to sustainable land management planning on both regional and national scale.

Background

Swidden, also known as slash-and-burn agriculture and shifting cultivation, is a common land use in mountainous regions, especially in Southeast Asian countries (Das et al. 2021; Schmidt-Vogt et al. 2009; Mertz et al. 2009a; Mertz et al. 2009b; Li and Feng 2016). Traditionally, the land is cleared for cultivation (normally by fire) on a rotational basis and then left for a few years for regeneration. This type of land use can be largely found on the forest-agriculture frontiers, and is still carried out in 40– 50 countries (Mertz 2009). According to Heinimann et al. (Heinimann et al. 2017), shifting cultivation landscapes, including fallows cover roughly 280 million hectares worldwide. In the Southeast Asia regions alone, about 14–34 million people engage in swidden cultivation (Mertz et al. 2009b). Therefore, it is self-evident that swidden agriculture plays an important role in mountainous communities and the livelihoods of the dwellers.

According to the reviews by Schmidt-Vogt et al. (Schmidt-Vogt et al. 2009) and van Vliet et al. (Vliet et al. 2012), the conclusive data on the extent of swidden is surprisingly limited, especially in Southeast Asia. The swidden system itself consists of temporarily cultivated land and associated fallows, and therefore it is a complex and dynamic land use. According to the review by Schmidt-Vogt et al. (Schmidt-Vogt et al. 2009), swidden fields are identified as agricultural lands, and fallows are grouped as either unclassified and/or barren wasteland or wood- and shrub-lands. The data gatherers find it difficult to define, find and measure the swidden system, and therefore often relegate it to a “residual category” of land use (Schmidt-Vogt et al. 2009; Padoch et al. 2007). Consequently, swidden lands often do not appear on land use maps or in statistical records (Schmidt-Vogt et al. 2009). Mapping swidden land use and its dynamics plays an important role in policy making and proper resource management for socio-economic and environmental benefits (Mertz et al. 2009a; Vliet et al. 2012). Yet, because of political control over swidden farming as well as the complexity and dynamics of swidden smallholdings, delineating swidden agriculture remains challenging (Padoch et al. 2007).

Several efforts to estimate the extent of swidden on a regional basis (Richards and Flint 1994; Uhlig et al. 1994) and on national scale in Southeast Asian countries have been made in the past, using different data sources such as censuses, forest inventories, aerial photographs and satellite images (Schmidt-Vogt et al. 2009). Attempts were also made to estimate solely on secondary forests including various phases of swidden fallows in Lower Mekong Region by Heinimann et al. (Heinimann et al. 2017) and Mittelman (Mittelman 2021), using various datasets. The first result came out of the regional analysis of swidden in southeast Asia by Spencer (1966) but with limited resources. Based on this result, a map of swidden distribution in the Asia–Pacific region was depicted in the Conservation Atlas of Tropical Forests: Asia and Pacific by Collines et al. (Collines et al. 1991).

Today, thanks to technological development, many kinds of remote sensors including airborne and space sensors provide remotely collected data for assessing burned areas like swidden patches (Lentile et al. 2006). At present, several kinds of remotely sensed data are available in different resolutions on different platforms including Landsat. With a growing availability of remotely sensed data (Potapov et al. 2019; Potapov et al. 2012), the extent, dynamics, and spatial characteristics of swidden agriculture are quantified on different scales in Asia–pacific countries (Castella et al. 2013; Hurni et al. 2013; Liao et al. 2015; Molinario et al. 2017; Messerli et al. 2009).

Using remote sensing, burned areas can be mapped with various image classification methods, inclusive of visual analysis, single channel density slicing, multitemporal thresholding of vegetation indices, principal component analysis, regression modelling, supervised and unsupervised classification, and spectral mixture analysis. Recent studies identified burned forest areas of shifting cultivation using a threshold value method with great overall accuracy (Das et al. 2021; Das et al. 2022; Swe 2020). A range of spectral indices such as the Normalized Burn Ratio (NBR), the difference in the Normalized Burn Ratio between pre- and post-fire images (dNBR), and the Normalized Difference Vegetation Index (NDVI) are commonly used for producing burned area maps (Lentile et al. 2006; Miller and Thode 2007; Rozario et al. 2018). The differences between pre- and post-fire indices are extensively used for characterizing burned areas because such values provide objective results in detecting changes (Landmann 2003). Previous studies showed that the dNBR- and RdNBR-based classification provided a higher degree of accuracy in research results (Miller and Thode 2007; Rozario et al. 2018; Miller et al. 2009) while Das et al. (2021) observed a great performance of NDVI and NBR in capturing differences in the burned patches in Northeast India. Even though further technological development has been made to identify and follow the trends of land use transition, we still cannot fully analyze the complex, dynamic and fragmented land use system of swidden (Mertz et al. 2009a).

Specifically in Myanmar, based on the field inventory data and satellite imagery by the Forest Department, the estimate of swidden area ranged from 0.29 million hectares within Permanent Forest Estate (PFE) to 10.18 million hectares (Schmidt-Vogt et al. 2009). The extent of swidden in Myanmar is extremely variable, even though the figures are provided by the same government agency (Schmidt-Vogt et al. 2009). A recent global scale study found that the shifting cultivation area in Myanmar has declined drastically since 2000 (Heinimann et al. 2017). Some studies researched on swidden cultivation in Myanmar in terms of traditional knowledge of swidden cultivators (Thet and Tokuchi 2020; Win 2004), biomass status in swidden fallows (Chan et al. 2013; Chan and Takeda 2016), and mapping of swidden land use at village level (Swe 2020; Chan et al. 2013; Chan and Takeda 2016; Swe and Nawata 2020; Chan and Takeda 2019). According to Chan and Takeda (2016) and Thet and Tokuchi (2020), the swidden fallow periods vary in different places in the 7 to 20-year range. Also, Swe et al. (Swe 2020) found that the fallow period decreased when permanent agriculture was introduced, and the availability of land for cultivation became limited as well. In the meantime, teak plantation was introduced in the Bago Mountains, Myanmar in the nineteenth century (Bryant 1996), and has since been well established. It is likely that the advent of plantation practices has made it difficult to assess the distribution of swidden at the national level. To fill this gap in swidden area mapping, we attempted to preliminarily estimate the extent of swidden in 2016 at the national level.

Objectives

The overall objective of this study is to produce a swidden land use/ land cover map of Myanmar by developing a decision tree-based swidden land use detection method. The national-level swidden land use data are scarce and thus urgently needed for the reform of the country’s land use policies pertaining to the land use rights of the mountain dwellers, as well as for the planning of sustainable development at local level. Therefore, this study is primarily intended to (i) delineate the swidden plots by developing a decision tree-based detection method, and with the use of remotely sensed data, (ii) test the accuracy of the method using the ground truth, and (iii) generate data on swidden land use areas of the years 2015‒2016.

Methods

Study area

Myanmar with an area of 676,557 sq. km has 15 administrative regions: one Union Territory (Nay Pyi Taw), seven States (Chin, Kachin, Kayin, Kayah, Mon, Rakhine, and Shan), and seven Regions (Ayeyarwady, Bago, Magway, Mandalay, Sagaing, Tanintharyi, and Yangon). The country lies between latitudes 9° and 29°N, and longitudes 92° and 102°E (Fig. 1), bordering Bangladesh and India in the west, Tibet Autonomous Region in the north, and Yunnan, Laos and Thailand in the east. Many mountain ranges such as the Rakhine Yoma, the Bago Yoma, the Shan Hills and the Tenasserim Hills run north-to-south. Tropical monsoon climate and geographical modification characterize various ecosystems of Myanmar, featured by ice-capped mountains in the north, desert-like dryland in the center, and tropical moist region in the south.

Fig. 1

Description of the study site with administrative boundaries

The agricultural sector is the main economy of Myanmar, and 75% of the total population lives in rural areas (Win 2004; Department of Population 2019). Over the mountainous regions, several ethnic groups (Kachin, Kayah, Karen, Chin, Mon and Shan) have been practicing swidden agriculture for generations. Each ethnic group has their own rules and practices on their swidden system. Agricultural crops such as upland rice, cotton, sesame, chili peppers, and groundnuts are grown for both subsistence and income generation in the first year of cultivation. Swidden cultivation plots are selected at the end of winter, then slashed and burned before the rains (i.e., before May). Then, swiddeners start growing crops. They weed during the rainy season, and harvest the crops in the winter (i.e., from the beginning of November) (Swe 2020; Thet and Tokuchi 2020; Win 2004; Chan et al. 2013; Chan and Takeda 2016; Chan and Takeda 2019). After crop cultivation, the cultivated areas are left as fallows for 7–20 years to recover for revegetation. The duration varies from location to location. (Swe 2020; Thet and Tokuchi 2020; Win 2004; Chan et al. 2013; Chan and Takeda 2016; Chan and Takeda 2019; Swe and Funakawa 2019). According to the previous studies (Swe 2020; Thet and Tokuchi 2020; Win 2004; Chan et al. 2013 Sep; Chan and Takeda 2016; Chan and Takeda 2019; Swe and Funakawa 2019) and field observation, the swidden cultivation in Myanmar could be found in the elevation above 150-m, and slope of > 10 percent.

In-situ data collection

Swe and Nawata (2020) tracked the boundaries of all existing swidden fields of SN and KC villages in the Bago Yoma using GPS (German Darkota 20) in May and June, 2016. Also, they collected the point information of lands for other uses such as forest lands in the studied villages. Furthermore, they conducted interview surveys with the locals to confirm the size of the swidden fields, their agricultural calendar and fallow period.

Using those ground truth data of Swe and Nawata (2020), we assessed the accuracy of the decision tree-based detection method developed in this study.

Landsat image selection and pre-processing

Previous studies recommend using the Landsat images to perform studies on burned areas due to its sufficiently high spatial resolution (Torralbo and Benito 2012) and free availability. In this study, we applied Landsat 8 Operational Land Imager (OLI)/ Thermal Infrared Sensors (TIRS) satellite images with a 30 m spatial resolution from the Collection 1 Tier 1 Surface Reflectance Data, as the prime source of data for the study. According to Foga et al. (2017), these images have the lowest geo-registration errors which include a pixel Quality Assessment (QA) band based on CF Mask (C code based on the Function of Mask or Fmask) algorithm, and is given preference to for operational cloud and cloud shadow detection. In addition, the USGS Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM) was applied to perform the landscape analyses of elevation and slope position.

Pre-processing CF Mask of Landsat 8 OLI/TIRS images requires application of quality bands to remove cloud and cloud shadow contamination by using respective QA band values. Pre-processing with the application of Landsat images and the landscape analyses with the SRTM DEM were performed on the Google Earth Engine (GEE) platform (Fig. 2).

Fig. 2

Flow chart with decision tree-based detection method

Image composites

In this study, several Landsat images acquired during April‒May in both 2015 and 2016 were applied for the accuracy and national level analyses (Fig. 2). When there are more than one image per region of interest, composites were created by using the median values of the Blue, Green, Red, Near Infrared, Short Wave Infrared-1 and Short Wave Infrared-2 bands of the images. In total, we applied 14 and 16 images from 2015 and 2016, respectively, for the accuracy test sites (Fig. 3). At the national level, we used 178 and 181 images from 2015 and 2016, respectively (Fig. 3).

Fig. 3

The number of Landsat 8 images used and Cloud Cover Percentage for the test site (upper panel) and for the whole country (lower panel)

Vegetation indices

To detect swidden plots, several vegetation indices were applied. In this study, two frequently used vegetation indices in swidden land use studies (Lentile et al. 2006; Miller and Thode 2007; Rozario et al. 2018; Landmann 2003; Miller et al. 2009) were computed to detect burnt areas: NBR and NDVI. NBR is one of the most widely used vegetation indices to highlight burnt areas in large fire zones. NDVI is also one of the most common vegetation indices to grasp the repetition and alternation pattern of swidden cultivation over different temporal intervals in test site and to characterize vegetation phenology.

Using the image composites, the pre- and post-NBR were calculated by Eq. (1). To know the changes in NBR values due to vegetation clearance or vegetation loss and burn, the dNBR value was computed by Eq. (2). Also, we recalculated the relative difference NBR (RdNBR) by Eq. (3).

$$NBR=(NIR-SWIR2)/(NIR+SWIR2)$$

(1)

$${\text{dNBR }} = {\text{ Pre}} - {\text{NBR }}{-}{\text{ Post}} - {\text{NBR}}$$

(2)

$${\text{RdNBR}}\, = \,{\text{dNBR}}/ \, \left[ {\sqrt {\left\{ {{\text{ABS}}\left( {{\text{pre}} - {\text{fireNBR}}/{1}000} \right)} \right\}} } \right]$$

(3)

Similarly, to assess the reduction in vegetation presence, the pre- and post-NDVIs were computed by Eq. (4), and the difference NDVI (dNDVI), by Eq. (5).

$$NDVI=(NIR-Red)/(NIR+Red)$$

(4)

$${\text{dNDVI }} = {\text{ Pre}} - {\text{NDVI }} - {\text{ Post}} - {\text{NDVI}}$$

(5)

The differences in NBR and NDVI and the relative difference in NBR were adjudged to account for any change pixels due to vegetation clearance or vegetation loss and burn for swidden cultivation. Such changes in index value ranges were optimized based solely on their occurrences in various images acquired during the April‒May of 2015 and 2016 periods (Fig. 2). After computing all the indices, we overlaid on them the elevations and the slopes derived from SRTM DEM images.

Threshold value extraction and selection

Within the ground truth polygons of the swidden plots collected in 2016, sample points were randomly taken, especially at the centers of the ground truth polygons. In addition, outside the ground truth site (test site), we took sample points by way of visual interpretation, comparing the vegetation statuses in 2015 and 2016 (Fig. 4). A total of 108 sample points in and outside the test sites were applied to extract the dNBR, RdNBR, and dNDVI values from the stacked layers (Fig. 2). The extracted values from the sample points were then calculated to obtain minimum, 1st Quartile, 2nd Quartile, and 3rd Quartile, and maximum values for each index (except elevation and slope). Therefore, a total of 15 threshold values of dNDVI, dNBR, and RdNBR indices with their corresponding minimum, 1st Quartile, 2nd Quartile, and 3rd Quartile and maximum values were assessed to select the most suitable threshold for each index for further accuracy assessment.

Fig. 4

Taking sample points using Landsat Image composites (April–May, 2016) in (upper panel) and outside (lower panel) ground truth sites

Before accuracy assessment, we considered another factor, i.e. the size of the swidden plot. Due to different patterns of swidden practice (such as random shift by the Karen people, and collective shift by some ethnic groups of the Chin), the extent of swidden varies in the number of plots and in area. Also, at the time of plot-burning for swidden cultivation, site preparation for tree plantation takes place in Myanmar, especially in mountainous areas, which may have been misdetected (wrongly detected) as swidden plots. However, the size of tree plantation sites is often bigger than that of swidden plots. According to the Forest Department (unpublished), the smallest size of tree plantation is 25 acres (about 10.12 ha). The smallest size of swidden plot is 0.5 acres (0.2 ha) based on the ground truth data of (Swe 2020). Therefore, we excluded the areas, which are smaller than 0.5 acres and larger than 25 acres, to extract the actual area of swidden plots.

When assessing the five thresholds for each of the three indices (in total 15 threshold values), we adopted the decision tree-based detection (DTD) method to select the most accurate threshold for each index. The DTD method is a commonly used data mining method for developing classification systems based on multiple covariates or for establishing prediction algorithms for a target variable (Song and Lu 2015). We conducted an accuracy assessment using the DTD method based on the ground truth data collected in 2016. The accuracy was assessed for two parameters (number of plots, and the extent of area) by using the following Eqs. 6 and 7 (Swe 2020):

$$Accuracy\left( \% \right),\,based\,on\,the\,number\,of\,plots = \frac{{100 \times total\,number\,of\,accurated\,etection\,plots\,identified\,by\,ground\,truth}}{{total\,detected\,area\,(Number\,of\,plots)}}$$

(6)

$$Accuracy\left( \% \right),\,based\,on\,the\,area = \frac{{100 \times total\,accurate\,detection\,are\,a\,identified\,by\,ground\,truth}}{{total\,detected\,area\,(ha)}}$$

(7)

From this accuracy assessment, we selected three criteria with the most appropriate corresponding thresholds.

Threshold value pairing and accuracy assessment

We combined the three indices (dNDVI, dNBR, and RdNBR) and their corresponding thresholds of the highest accuracy, with both elevation and slope percent derived from SRTM DEM to make the threshold value pairs. Seven pairs of indices (Pair 1–7) were obtained. Again, we performed a similar accuracy assessment on these seven pairs of indices, using the ground truth polygons collected in the test site, to select the most accurate pair for further assessment of swidden plots at national level.

Estimation of swidden extent at national level

Using the best-fit model, we estimated the extent of swidden for the whole of Myanmar. It should be noted that cloud cover was not taken into account in this estimation. Additionally, we calculated the swidden area by the administrative boundaries to know which regions need to be paid attention to in terms of land use planning and sustainable development of swidden system at the regional level. Due to cloud cover, we computed the area percentage of each region by the total swidden area detected in this study.

Results

Threshold determination and accuracy

After pre-processing, compositing and analyzing vegetation indices for the pre- and post-burn images, we extracted the threshold values for NBR and NDVI indices with data from 108 sample points. These threshold values for each index were statistically extracted by quartiles, and the accuracies for each pair of index and threshold assessed by two parameters _ the number of plots and the area of the detected plots.

The threshold values for dNBR are 0.058 at minimum threshold (dNBR_1), 0.297 at 1st quartile (dNBR_2), 0.371 at 2nd quartile (dNBR_3), 0.440 at 3rd quartile (dNBR_4), and 0.588 at maximum (dNBR_5) (Table 1). The threshold values at minimum (RdNBR_1), at 1st quartile (RdNBR_2), at 2nd quartile (RdNBR_3), at 3rd quartile (RdNBR_4), and at maximum (RdNBR_5) are 8, 18, 21, 23, and 29, respectively (Table 1). The threshold values for dNDVI are 0.032 at minimum (dNDVI_1), 0.158 at 1st quartile (dNDVI_2), 0.235 at 2nd quartile (dNDVI_3), 0.303 at 3rd quartile (dNDVI_4), and 0.515 at maximum (dNDVI_5) (Table 1).

Table 1 Threshold values and accuracy assessment

The dNBR models produced 752, 94, 58, 36 and 1 plots by each threshold. According to the assessment based on the number of plots, dNBR_5 threshold produced one plot which overlapped with the ground truth plot (100% accuracy), followed by dNBR_3 threshold (79.31% accuracy, 46 overlapped plots), dNBR_2 threshold (76.60% accuracy, 72 overlapped plots), dNBR_4 threshold (75.00% accuracy, 27 overlapped plots), and dNBR_1 threshold (9.27% accuracy, 53 overlapped plots) (Table 1). The assessment based on the area of the ground truth provided 100% accuracy for dNBR_5 threshold, followed by dNBR_2 threshold (77.67%), dNBR_4 threshold (77.53%), dNBR_3 threshold (77.20%), and dNBR_1 threshold (65.44%) (Table 1). Although the dNBR_5 threshold was found to be 100% accurate in term of the number of plots, it could detect only one accurate plot with the area of 0.51 ha which did not cover the ground truth area of 129.42 ha. Also, in terms of the number of plots detected and overlapped with the ground truth plots, dNBR_3 showed a higher percentage (79.31%) than dNBR_2 (76.60%) (Table 1). However, considering the area-based accuracy with a higher number of plots detected by the model and overlapped with the ground truth plots, we selected the dNBR_2 threshold for further analysis.

In case of RdNBR index, the models gave 310, 112, 98, 91 and 26 plots for each threshold (Table 1). The assessment based on the number of plots showed the highest accuracy (71.43%, 65 overlapped plots) with RdNBR_4 threshold, followed by 66.33% (65 overlapped plots) with RdNBR_3, 65.38% (Potapov et al. 2019) with RdNBR_5, 62.50% (70 overlapped plots) with RdNBR_2, and 22.58% (70 overlapped plots) with RdNBR_1 (Table 1). The area-based assessment also gave the highest accuracy (81.33%) with RdNBR_4 threshold, followed by 80.07% with RdNBR_3 threshold, 77.59% with RdNBR_2 threshold, 65.27% with RdNBR_5 threshold, and 52.03% with RdNBR_1 threshold (Table 1). Therefore, we applied the RdNBR_4 threshold for the later analysis, based on the accuracy assessment in terms of the number of plots detected and the area extent of the overlapped plots.

The dNDVI models detected 682, 124, 56, and 16 plots for each threshold, except dNDVI_5 (Table 1). The threshold dNDVI_5 could not detect any accurate plots. According to the assessment based on the number of plots, it showed the highest accuracy (73.21%, 41 overlapped plots) with dNDVI_3 threshold, followed by dNDVI_2 threshold (62.10%, 77 overlapped plots), dNDVI_4 threshold (56.25%, 9 overlapped plots), and dNDVI_1 threshold (6.16%, 42 overlapped plots) (Table 1). The area-based assessment for the dNDVI index showed the highest accuracy (73.55%) with dNDVI_3 threshold, followed by dNDVI_2 threshold (71.70%), dNDVI_1 threshold (68.77%), and dNDVI_4 (63.26%) (Table 1). The results showed that the dNDVI_3 produced the highest accuracy based on the number of plots and area assessments.

Therefore, we performed further analyses using the threshold values of 0.297 for dNBR, of 23 for RdNBR, and of 0.235 for dNDVI (Table 1).

Index pair analysis and accuracy assessment

Using the threshold values of dNBR, RdNBR and dNDVI, indices which have the highest accuracy, and together with elevation and slope, we performed analyses for seven index pairs (Table 2). Like the accuracy assessment in threshold determination, we calculated the accuracy, in terms of the number of plots and area, for the seven threshold pairs (Table 2).

Table 2 Description and accuracy assessment for seven index pairs

When assessing the index pairs, different models detected different numbers of plots. There were 97 plots produced by Pair_1 model, 88 by Pair_2, 76 by Pair_4, 55 by Pair_3 and Pair_5, and 48 by Pair_6 and Pair_7 (Table 2). Figure 5 also depicted the generated polygons by seven index pairs (blue) and the ground truth polygons (red) which visualize the accuracy assessment.

Fig. 5

Maps showing the accuracies for [1] dNBR + Slop + Elevation (Pair_1), [2] RdNBR + Slope + Elevation (Pair_2), [3] dNDVI + Slope + Elevation (Pair_3), [4] dNBR + RdNBR + Slope + Elevation (Pair_4), [5] dNBR + dNDVI + Slope + Elevation (Pair_5), [6] RdNBR + dNDVI + Slope + Elevation (Pair_6), and [7] dNBR + RdNBR + dNDVI + Slope + Elevation (Pair_7)

According to the accuracy assessment by the number of model-produced plots which overlapped with the ground truth plots, the Pair_4 model showed the highest accuracy (82.89% with 63 overlapped plots), followed by Pair_2 (77.27% with 68 overlapped plots), Pair_6 and Pair_ 7 (77.08% with 37 overlapped plots each), Pair_3 and Pair_5 (74.55% with 41 overlapped plots each), and Pair_1 (74.23% with 72 overlapped plots) (Table 2).

The area-based accuracy assessment showed that the Pair_4 model produced the largest accuracy (84.25%), followed by Pair_2 (82.04%), Pair_1 (78.98%), Pair_6 and Pair_7 (75.20% each), Pair_3 (72.72%), and Pair_5 (72.69%) (Table 2).

According to both accuracy assessments, the Pair_4 model showed the highest accuracy, and we will use this model to estimate the extent of swidden for the whole country.

Estimation of the extent of swidden at national level

Using the Pair_4 model, the spatial extent of swidden plots was estimated at national level for the year 2016. This, however, does not include some regions due to an extensive cloud cover during the burning and leaf initiation periods of the study year. The cloud cover was 8379623.25 ha, especially in the Northern parts of Myanmar such as Kachin and Sagaing. The spatial distribution of swidden plots and cloud cover is depicted in Fig. 6. The average proportions of cloud cover for the test site in 2015 and 2016 are 32.45% and 21.75%, respectively (Fig. 3). The average proportions of cloud cover for the whole site in 2015 and 2016 are 33.20% and 30.61% (Fig. 3).

Fig. 6

Generated Swidden Land Cover Map Classified by Decision Tree Detection Method for Myanmar

According to the estimate, the swidden area covered about 100311.10 ha in 2016. It comprises a very large proportion of the total land area of the country (about 0.15% of the total country area). If we exclude the cloud cover area from the total country area, the estimated swidden plots encompass about 0.17% of the total national area.

According to the estimate of the extent of swidden by the administrative boundaries in Myanmar (Table 3), swidden is estimated to be practiced in almost every state and region, except the Yangon Region. A disproportionately large area of swidden cultivation was observed in Shan State with an area of 57733.90 ha (57.56% of the total detected swidden area in 2016), followed by Chin State (13409.73 ha, 13.37%).

Table 3 Swidden Land Cover Areas of Myanmar

Discussion

In Myanmar, due to extremely variable figures on swidden extent, even from the same governmental department, we attempted to estimate the swidden at national level using free available Landsat images in combination with decision tree-based plot detection method in GEE platform.

Remotely sensed data, indices, thresholds, and models for detecting swidden system

To assess the overall system of swidden, several attempts were made in a number of studies (Das et al. 2021; Lentile et al. 2006; Swe 2020; Landmann 2003). Due to the constraints of monitoring techniques and data at macroscopic scale, previous geographical studies of swidden agriculture used to be mainly descriptive (Inoue 2000). The advancement of remote sensing techniques since 1970s, however, has made it possible to monitor environmental changes ever more closely, in that remotely sensed data provide invaluable information on fire events and burned area with its synoptic, multi-temporal, multi-spectral and repetitive coverage capabilities (Vadrevu and Justice 2011). With a growing availability of remotely sensed data (Potapov et al. 2019; Potapov et al. 2012), the extent, dynamics, and spatial characteristics of swidden agriculture are quantified at different scales in Southeast Asian countries (Heinimann et al. 2007; Castella et al. 2013; Hurni et al. 2013; Liao et al. 2015; Molinario et al. 2017). Singh and Dubey (Singh and Dubey 2012) suggested that remote sensing provides land resource data in both digital form and different bands of the electromagnetic spectrum. Satellite data that can be used for delineating swidden cultivation are Medium spatial resolution Landsat images such as Landsat TM, Enhanced Thematic Mapper (ETM), and Operational Land Imager (OLI), coarse spatial resolution by MODIS and Synthetic Aperture Radar (SAR). Medium spatial resolution Landsat images were the most commonly used satellite data in many previous studies (Li et al. 2014). They are particularly suitable for mapping burned areas (Li et al. 2014; Petropoulos et al. 2011). Therefore, we utilized free available Landsat imageries in this study to detect swidden patches.

Regarding the index selection, Swe and Nawata (2020) applied the NDVI break value to detect swidden plots in the Bago Yoma, Myanmar. They developed remote sensing-based detection methods for swidden land use with the application of the Landsat images. Das et. al (2021) also used, dNBR, RdNBR and dNDVI indices to develop ways to observe long-term mapping of shifting cultivation. With the application Landsat data. they developed a decision tree-based multi-step threshold (DTMT) method in Northeast India. Thereafter, other studies applied the same vegetation indices to assess swidden agriculture and its trend (Lentile et al. 2006; Miller and Thode 2007; Rozario et al. 2018). In this study, we also applied these three common indices (dNBR, RdNBR and dNDVI), and the most appropriate thresholds figured out by a decision tree-based detection method.

Our accuracy assessment on three indices showed that dNBR_2 has the highest accuracy (76.60%), followed by dNDVI_3 (73.21%) and by RdNBR_4 (71.43%) in terms of the number of overlapped plots. In the aspect of swidden area, RdNBR_4 also showed the highest accuracy (81.33%), followed by dNBR_2 (77.67%) and by dNDVI_3 (73.55%). In previous studies, dNBR- and RdNBR-based classifications showed a higher degree of accuracy in their research results (Miller and Thode 2007; Rozario et al. 2018; Miller et al. 2009) and NDVI and NDVI difference also provided a great performance in capturing burned areas in Northeast India (Das et al. 2022).

In addition, Swe and Nawata (2020) suggested that the maximum likelihood method was better than the NDVI method for swidden land use detection. In a case study by Das et al. (2021), the combination of dNDVI, dNBR, and RdNBR produced higher overall accuracy. However, in this study, NBR showed the better accuracy than NDVI. Of seven pairs, dNBR in combination with RdNBR produced the highest accuracy (84.25%) compared to any of these indices alone with slope and elevation, and among all pairs of indices. Therefore, the applicability of each index needs to be further assessed.

Swidden plot detection and accuracy assessment

Previous studies stated that the complex nature of swidden made it difficult to achieve an accuracy higher than 70% by using Landsat images (Castella et al. 2013; Liao et al. 2015; Müller et al. 2013; Thatheva and Yasuyuki 2009). However, recent studies in Northeast India by Das et al. (2021), in which the spatiotemporal dynamics of swidden cultivation was assessed using Landsat images and applying the decision tree multi-threshold classification framework based on vegetation indices such as dNDVI, dNBR, and RdNBR, showed the higher overall accuracy (above 85%). This study also produced the highest accuracy of 84.25% and the rest are also higher than 72%. These degrees of accuracy could have resulted from the application of image composites during the study period. Similarly, Das et al. (2021) and Swe and Nawata (2020) pointed out that multi-image classifications can increase the accuracy of swidden detection methods.

Regarding the accuracy assessment on these selected indices with corresponding most appropriate thresholds, Swe and Nawata (2020) presented the area-based accuracy of 62.8% in SN village, and 42.2% in KC village.

In this study, we applied a similar accuracy assessment method to Swe and Nawata (2020), and our decision tree-based detection method produced a relatively higher accuracy result. Similarly, Shimizu et al. (2017) conducted a swidden plot level detection test in the Bago Mountains, Myanmar, with Landsat images using a area-based accuracy assessment method. In their findings, the producer’s and user’s accuracy levels for shifting cultivation during the period from 2000 to 2014 were 45.1% and 74.4%, respectively (Shimizu et al. 2017). Although Sakai (2002) conducted a similar analysis to accurately detect “slash and burn” cultivation plots by considering the elements ˗ (i) the decrease in vegetation between two periods (ii) sizes of plots (iii) shapes of plots, the accuracy assessment results were not reported. In our study, the area-based accuracy was relatively higher compared to the accuracy result of Shimizu et. al. (2017) and Swe and Nawata (2020). Even though the accuracy levels vary, Das et al. (2022) demonstrated an overall accuracy of above 85% and Kappa values of above 0.69, for all four periods (a) 1975–1976, (b) 2000–2001 (c) 2014–2015 (d) 2017–2018. The present study also produced the highest accuracy of 84.25% and the rest are also higher than 72%.

Therefore, the selection of indices, an appropriate threshold value for each index by DTD method provided better accuracy compared to other similar studies, although similar accuracy was given when assessed with different methods.

Moreover, to fully detect all ground truth swidden fields quite challenging in this study. Similar finding was reported by Swe and Nawata (2020), even the combination of the NDVI break value method and the Maximum Likelihood method allowed just over 70% of the actual swidden plots to be detected. Such detection failures were possibly due to incomplete burned fields and small-sized swidden plots as well as the topography (Swe 2020). It is also possible that the most appropriate threshold was not selected for each of the applied methods. However, our study suggested that the use of ground truth data in combination with a decision tree-based detection method produces better accuracy in detecting swidden plots on the study sites and proved the applicability of our methodology to detect highly complex swidden plots.

Swidden extent at national level

This study is the first attempt of its kind to estimate the swidden extent at a national level in Myanmar, using free available Landsat images and the GEE platform. The results showed that swidden cultivation is mostly practiced in areas where ethnic minority groups live. According to this study, swidden was most extensive in Shan State, followed by Chin State. Win (2004) reported similar findings that swidden cultivation practice mostly prevails in Chin and Shan states.

Without cloud cover areas, we estimated that there was a total swidden area of 0.1 million ha nationwide in 2016. Compared with other previously reported Figs. (0.29 – 10.18 million ha), our estimate is much smaller likely due to an actual reduction in swidden land use and its demise in some places due to several factors. Mertz et al. (2009b) and Schmidt-Vogt et al. (2009) reported some farmers downscaled swidden farming and others ceased it altogether. Swe and Nawata (2020; 2020) also presented a decreasing trend in swidden agricultural practice in two villages and the demise in one village that they studied. Whether or when such changes occurred in Myanmar is yet to be known and we need to make further time-series analyses using the same method for more details.

Limitations

Limitations in this study were heavy clouds and other atmospheric constituents mostly in regions such as Sagaing and Kachin. CF mask products were applied to exclude cloud pixels prior to image classification, and these pixels were treated as no data values. Even though the image classification was done and change detection was modelled based on the best classification result on the GEE platform, we could not avoid some degree of cloud cover when we detected the swidden patches at national level (Fig. 7).

Fig. 7

Land 8 Composite images (April–May 2016, SWIR-NIR-RED Band combination) with the swidden plots detected by the best-fit model

Also, higher mapping accuracy would likely be achieved with larger ground truth data (Anders et al. 2021; Maglogiannis et al. 2007). The ground truth data utilized during the testing phase in this study were limited in quantity and involved positional uncertainty, even though we collected them form 108 sample points in different parts of the country. A challenge we encountered was discrepancies in time between the available CF mask of Landsat 8 OLI/TIRS images and the ground truth data. Therefore, the accuracy needs to be further assessed with more ground truth data from other parts of the country. Despite all these limitations, the proposed method and the accuracy assessment provided a better performance in detecting swidden plots in Myanmar than ever before.

Conclusions

This study attempted to fill the gap of swidden land cover map of Myanmar at a national level. Our primary goal was to develop a swidden land cover map for the whole of Myanmar by applying free available Landsat images, and using GEE cloud computing in combination with the DTD method of several indices. When we applied the DTD method to select the most appropriate indices, and to find the best suited threshold values for each index, our study suggested that the dNBR model combined with RdNBR, slope and elevation produced the highest accuracy (83% and 84% in terms of the number of plots and the area of the overlapped plots, respectively), while the all index pair (dNBR, RdNBR, dNDVI, slope, and elevation) showed lower accuracy even with the ground truth data. Therefore, NBR is likely to outweigh NDVI in detecting burned areas like swidden patches, though we may need to further assess the applicability of different indices. Using the best-fit combinations of the indices and thresholds (Pair_4 which is the combination of dNBR > 0.297 and RdNBR > 23 and slope > 10 and elevation > 150) with high accuracy, we could produce an overall national swidden land use map of Myanmar for 2016.

In addition, our estimate on the extent of swidden showed that about 100,311.10 hectares of land was used for swidden in Myanmar in 2016 with the highest extents in Shan State and Chin State. It showed a smaller Fig. (0.1 million ha), compared to other previously reported data (0.29 – 10.18 million ha). This is likely due to transition away from, and in some areas, demise of swidden. The other factor is cloud cover. There was a very large proportion of cloud, especially in the Sagaing Region and Chin State. Without the cloud cover area, the estimated swidden plots comprised about 0.17% of the total national land area.

In conclusion, this study contributed to fill the gap between the actual swidden area and the area depicted on an existing available map of Myanmar with the use of free available images and the GEE platform. As Myanmar is the largest country in mainland Southeast Asia in area with a great majority of the population living in rural areas, and many in the mountains, its land resources are of great relevance to the people’s livelihoods and thereby the nation’s progress. This study using the DTD method with the appropriate vegetation indices and thresholds, also provided the potential to estimate swidden extents both at regional and national levels. This will largely contribute to land management planning, especially for the farmers in the mountainous areas.