Extent and Area of Swidden in Montane Mainland Southeast Asia: Estimation by Multi-Step Thresholds with Landsat-8 OLI Data

Information on the distribution, area and extent of swidden agriculture landscape is necessary for implementing the program of Reducing Emissions from Deforestation and Forest Degradation (REDD), biodiversity conservation and local livelihood improvement. To our knowledge, explicit spatial maps and accurate area data on swidden agriculture remain surprisingly lacking. However, this traditional farming practice has been transforming into other profit-driven land use, like tree plantations and permanent cash agriculture. Swidden agriculture is characterized by a rotational and dynamic nature of agroforestry, with land cover changing from natural forests, newly-cleared swiddens to different-aged fallows. The Operational Land Imager (OLI) onboard the Landsat-8 satellite has visible, near-infrared and shortwave infrared bands, which are sensitive to the changes in vegetation cover, land surface moisture content and soil exposure, and therefore, four vegetation indices (VIs) were calculated, including the Normalized Difference Vegetation Index (NDVI), the Normalized Difference Moisture Index (NDMI), the Normalized Burn Ratio (NBR) and the Soil Adjusted Vegetation Index (SAVI). In this study, we developed a multi-step threshold approach that uses a combination of thresholds of four VIs and local elevation range (LER) and applied it to detect and map newly-opened swiddens and different-aged fallows using OLI imagery acquired between 2013 and 2015. The resultant Landsat-derived swidden agriculture maps have high accuracy with an overall accuracy of 86.9% and a Kappa coefficient of 0.864. The results of this study indicated that the Landsat-based multi-step threshold algorithms could potentially be applied to monitor the long-term change pattern of swidden agriculture in montane mainland Southeast Asia since the late 1980s and also in other tropical regions, like insular Southeast Asia, South Asia, Latin America and Central Africa, where swidden agriculture is still common.


Introduction
Swidden, a term invented by the Swedish anthropologist Karl Gustav Izikowitz, is resurrected to describe the age-old slash and burn agriculture in the tropics, in the context of the global collaborative program of the Reducing Emissions from Deforestation and Forest Degradation (REDD) in the developing countries [1][2][3]. In the uplands of mainland Southeast Asia (MSEA), swidden practice is typically household-based with a small size and randomly scattered [4,5], characterized by the agro-forestry's dynamics, complexity and diversity [6,7]. Being a unique and ubiquitous land use system, swidden agriculture provides livelihood necessities and maintains cultural identity for local ethnic groups [2,8,9]. Ever since the first call to overcome swiddening [10], governments in Southeast regions of MSEA. The awkward situation indicates that while the historical Landsat archive products are freely available to the public, on the one hand [23,27], the extent of ubiquitous swidden agriculture still remains ambiguous, on the other. Because of the rotational nature, the existing studies of swidden agriculture detection primarily utilize the single-date method, particularly by means of land use classification [28,29].
In this study, we developed algorithms for detecting and mapping swidden agriculture by comprehensively comparing the relationship of various VIs derived from Landsat OLI among the major land cover types, in combination with other environmental information in ArcGIS 10.2. Earlier studies have shown that other VIs beyond the Normalized Difference Vegetation Index (NDVI) [30] improve the delineation of burned surface [31,32]. Thus, the NDVI, the Normalized Difference Moisture Index (NDMI) [33] or the Land Surface Water Index (LSWI) [34], the Normalized Burn Ratio (NBR) [35] and the Soil Adjusted Vegetation Index (SAVI) [36] were obtained to discriminate the changes of land cover of the swidden system. Our objective is to establish a simple and effective method for delineating swidden in tropical highlands with limited optical imagery acquisition. The latest spatially-explicit information of swidden in 2014/2015 from this study will support forest management and REDD implementation.

Study Area
Montane mainland Southeast Asia (MMSEA), or the mountainous regions of MSEA, refers to the area over 300 m above sea-level (asl), covering the uplands of Cambodia, Laos, Myanmar, Thailand, Vietnam and Yunnan Province of China [15,37]. Five major river systems from west to east, namely the Irrawaddy, Salween, Chao Phraya, Mekong and Red rivers, run through MSEA from north to south. The region has a tropical monsoon climate with three distinct seasons, i.e., the cool-dry season (November-February), the hot-dry season (March-April) and the rainy season (May-October). The dry season (northeast monsoon) is characterized by low cloud cover with less than 20 mm of rainfall per month [38]. This greatly facilitates satellites (e.g., Landsat 8) to acquire cloud-free observations [25]. Typically, forests in MMSEA mainly comprise evergreen mountain forests (>1000 m¨asl), evergreen lowland forests (<1000 m¨asl), mixed deciduous forests and fragmented and degraded evergreen forest cover [38,39] and occupy the largest remaining tropical forests in MSEA [40]. Forest-dwelling people include various and diverse ethnic minority groups, and they have been extensively practicing swidden agriculture for subsistence for ages. Other land cover types include permanent farmland, waterbody, built-up land and swidden landscape (mosaics of newly-burned plots, swidden fields and fallows). Swidden agriculture used to be widely practiced in the uplands of Yunnan Province a long time ago. However, swidden agriculture in this province is currently distributed in the western, southwestern and southern border area (a long narrow belt) with Myanmar, Laos and Vietnam. In this study, as we focused on swidden agriculture, only the border area of Yunnan and the other five countries were viewed as the study area ( Figure 1).
Swidden agriculture generally starts with slashing and felling of natural or secondary forests in the early dry season, letting them dry by sun and air during the peak of the dry season and burning them to fertilize the sterile soils. Then, various food or cash crops are extensively grown before the rainy season comes [26]. After several years of crop cultivation, swidden fields experience a sharp decline of soil fertility and are left fallow for vegetation recovery. In the past few decades, under the pressure of population growth, economic development and national ecological conservation policies, swidden agriculture has undergone extensive transformations, including conversion into cash-crops (rubber, sugar cane, banana, etc.) cultivation and shortened fallow periods [3]. However, this traditional farming system remains a dominant land use category in MMSEA [6,41] and will maintain its viability into this century [5,42], although once ill-reputed for its deforestation and forest degradation [38,39]. This can be clearly demonstrated by the MODIS-derived fire counts (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)

Data Acquisition and Anomaly Values Processing
The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model Version 2 (ASTER GDEM V2), some of the latest publicly available topographical data, was jointly developed by the U.S. National Aeronautics and Space Administration and Japan's Ministry of Economy, Trade and Industry [43]. The ASTER GDEM covers the land surface between 83°N and 83°S with a finer spatial resolution (30 m). The DEM products were subjected to cloud masking, bad value and outlier removal and then partitioned into 1°-by-1° tiles, which contain at least 0.01% land area. Two-hundred seventy-four tiles were freely downloaded and then mosaicked to describe the topography of MSEA (Figure 2a).

Data Acquisition and Anomaly Values Processing
The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model Version 2 (ASTER GDEM V2), some of the latest publicly available topographical data, was jointly developed by the U.S. National Aeronautics and Space Administration and Japan's Ministry of Economy, Trade and Industry [43]. The ASTER GDEM covers the land surface between 83˝N and 83˝S with a finer spatial resolution (30 m). The DEM products were subjected to cloud masking, bad value and outlier removal and then partitioned into 1˝-by-1˝tiles, which contain at least 0.01% land area. Two-hundred seventy-four tiles were freely downloaded and then mosaicked to describe the topography of MSEA (Figure 2a). While the ASTER GDEM V2 has been improved substantially, particularly the removal of cloudy pixels, residual bad values and outliers [43], the products may still contain anomaly values at the site level. According to the maximum and minimum elevations of MSEA countries and Yunnan Province, China, the anomalies (negative values of DEM) were consistently assigned as zero (sea level) and the others (higher than the highest point) removed as no data, and then, values were allocated to the masked cells based on the Euclidean distance with the Nibble tool, a tool for value assignment with the value of their nearest neighbors in ArcGIS software.

Defining the Mountainous Area in MSEA
Vegetation fires in MSEA are caused by the upland swidden-related burning in most cases, yet they also include lowland agriculture residue burning [4,44]. With a view toward reducing the spectral noise effects of permanent farmland and the corresponding human settlements on the detection of swidden practice, MMSEA or the mountainous area of MSEA should be re-defined accordingly, as earlier studies primarily applied elevation (e.g., 300 m) to delineate mountain in MSEA [15,37]. Therefore, the mountain typology system proposed by the United Nations Environment Programme World Conservation Monitoring Center (UNEP-WCMC) was used to extract the mountainous part of MSEA, by considering the slope and terrain relief [45]. The Mountain Research Initiative typology was originally established at the global scale [45], which may fail to distinguish montane from non-montane pixels at the local scale or include scattered and isolated mountains [46,47]. Hence, for such an analysis at sub-global and smaller spatial scales, a parameterization of the appropriate local elevation range (LER) is needed [46].
In the context of MSEA, the mask comprises the land with different parametrization of the elevational and slope gradient and local relief (Table 1). Firstly, as a rule of thumb, slopes less than 2° and ranging from 2°-5° are termed as level or nearly level and very gentle slopes, respectively, and the corresponding lands are most appropriate for agriculture and living. Secondly, the minimum limit of elevation was set at 600 m to reduce the impacts of the upward permanent agriculture activities. Furthermore, some major cities and towns are located in areas below 600 m, such as While the ASTER GDEM V2 has been improved substantially, particularly the removal of cloudy pixels, residual bad values and outliers [43], the products may still contain anomaly values at the site level. According to the maximum and minimum elevations of MSEA countries and Yunnan Province, China, the anomalies (negative values of DEM) were consistently assigned as zero (sea level) and the others (higher than the highest point) removed as no data, and then, values were allocated to the masked cells based on the Euclidean distance with the Nibble tool, a tool for value assignment with the value of their nearest neighbors in ArcGIS software.

Defining the Mountainous Area in MSEA
Vegetation fires in MSEA are caused by the upland swidden-related burning in most cases, yet they also include lowland agriculture residue burning [4,44]. With a view toward reducing the spectral noise effects of permanent farmland and the corresponding human settlements on the detection of swidden practice, MMSEA or the mountainous area of MSEA should be re-defined accordingly, as earlier studies primarily applied elevation (e.g., 300 m) to delineate mountain in MSEA [15,37]. Therefore, the mountain typology system proposed by the United Nations Environment Programme World Conservation Monitoring Center (UNEP-WCMC) was used to extract the mountainous part of MSEA, by considering the slope and terrain relief [45]. The Mountain Research Initiative typology was originally established at the global scale [45], which may fail to distinguish montane from non-montane pixels at the local scale or include scattered and isolated mountains [46,47]. Hence, for such an analysis at sub-global and smaller spatial scales, a parameterization of the appropriate local elevation range (LER) is needed [46].
In the context of MSEA, the mask comprises the land with different parametrization of the elevational and slope gradient and local relief (Table 1). Firstly, as a rule of thumb, slopes less than 2˝and ranging from 2˝-5˝are termed as level or nearly level and very gentle slopes, respectively, and the corresponding lands are most appropriate for agriculture and living. Secondly, the minimum limit of elevation was set at 600 m to reduce the impacts of the upward permanent agriculture activities. Furthermore, some major cities and towns are located in areas below 600 m, such as Jinghong City. Thirdly, LER was set with three threshold values. Exploratory analysis shows a logarithmic relationship between the average of LER and neighborhood scales in Figure 3. Then, a 100-cell radius (or 3 km) was defined for the LER with the Focal Statistics tool under ArcGIS 10.2, hence using the LER of 300 m, 400 m and 500 m as thresholds. For the sake of brevity, we used LER300, LER400 and LER500 to represent the corresponding mountain area. MMSEA was defined ( Figure 2b) with an area of 1.06ˆ10 6 km 2 , 0.99ˆ10 6 km 2 and 0.90ˆ10 6 km 2 for LER300, LER400 and LER500, respectively, accounting for 45.34%, 42.51% and 38.73% of the total area of the five MSEA countries and Yunnan Province in China. It should be noted that the mountain area derived from the typology by UNEP-WCMC presents the theoretical distribution extent of mountainous area only from the aspect of physiography. However, the elevations lower than 600 m in the mountainous and hilly regions in China are preferentially considered for exclusive development purposes, such as city construction and permanent agriculture. As permanent farmlands and built-up land usually bring about noises for swidden detection, only the mountain area comprised of Classes 2, 3 and 4 is used to perform this study. 6 Jinghong City. Thirdly, LER was set with three threshold values. Exploratory analysis shows a logarithmic relationship between the average of LER and neighborhood scales in Figure 3. Then, a 100-cell radius (or 3 km) was defined for the LER with the Focal Statistics tool under ArcGIS 10.2, hence using the LER of 300 m, 400 m and 500 m as thresholds. For the sake of brevity, we used LER300, LER400 and LER500 to represent the corresponding mountain area. MMSEA was defined ( Figure 2b) with an area of 1.06 × 10 6 km 2 , 0.99 × 10 6 km 2 and 0.90 × 10 6 km 2 for LER300, LER400 and LER500, respectively, accounting for 45.34%, 42.51% and 38.73% of the total area of the five MSEA countries and Yunnan Province in China. It should be noted that the mountain area derived from the typology by UNEP-WCMC presents the theoretical distribution extent of mountainous area only from the aspect of physiography. However, the elevations lower than 600 m in the mountainous and hilly regions in China are preferentially considered for exclusive development purposes, such as city construction and permanent agriculture. As permanent farmlands and built-up land usually bring about noises for swidden detection, only the mountain area comprised of Classes 2, 3 and 4 is used to perform this study.

Landsat-8 OLI Imagery and Pre-Processing
In order to obtain the latest 30-m resolution spatial information of swidden agriculture, Landsat 8 OLI Level 1 terrain-corrected (L1T) products in MSEA during 2013-2015 are freely available from the USGS Earth Resources Observation and Science Center at the Global Visualization Viewer webpage [48]. The Landsat 8 L1T data have typically been processed with consistent radiometric calibration, systematic geometric correction and precision correction using the ground control data and a digital elevation model to collect parallax error caused by local topographic relief [23]. This pre-processing is essential for monitoring swidden agriculture in MMSEA. In addition, the Landsat 8 L1T data have improved geometric fidelity and geolocation accuracy [23]. A total of 120 OLI scenes acquired from mid-January-May (mainly March and April, 106 scenes) were applied to detect and map the latest swidden agriculture, as MODIS fire statistics showed that 75% of total Asian fire occurrences appeared between January and May, with the peak fire season in March [4]. For example, large fires are ubiquitous across the uplands of northern Laos in March, with yellow smoke rising from the slopes [11].
Landsat OLI Surface Reflectance High Level Data Products were ordered, processed and downloaded through the USGS Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) On Demand Interface at [49]. Products of the surface reflectance of OLI bands, cloud and cloud shadow (CS) mask, snow mask and VIs were obtained accordingly. Clouds and their shadows are other key issues to be considered. Among the 120 scenes, 101 of them have cloud coverage less than or equal to 5 percent, including 53 cloud-free ones. For the clouds and their shadow-contaminated pixels, the Fmask 3.2 standalone version was introduced to extract the corresponding mask [50]. The Fmask software deals with unmodified digital value imagery from the Landsat Archive and is widely applied in operational Landsat-related monitoring [51]. Four VIs, namely NDVI [30], NDMI [33], NBR [35] and SAVI [36], were calculated with the surface reflectance with the following equations: where ρ Red, ρ N IR, ρ SW IR1 and ρ SW IR2 refer to the surface reflectance values of the red band, NIR band, SWIR1 and SWIR2 bands in the Landsat-8 OLI sensor, respectively. For Equation (4), L is the soil brightness correction factor, and its value varies according to vegetation cover. For very high vegetation regions, L is assigned as zero, effectively turning SAVI into NDVI. For low to no vegetation areas, it is equal to 1. Typically, a value of 0.5 is the default setting for most situations, especially when vegetation cover is unknown.

Fieldwork on the Collection of Ground Truth Data in MSEA
Fieldwork for the collection of geotagged landscape photos on land use and cover types was carried out in MSEA from 2013-2015 ( Figure 1). During the field trips, the sampling strategy focused on deforestation and forest degradation induced by traditional swidden system (Figure 1a These photos observed and recorded specific attribute information about land use and cover types, such as primary or secondary forests, tree plantations (e.g., rubber, eucalyptus and teak), shrubland, grassland, swidden agriculture, croplands (e.g., paddy rice and cassava), villages and towns, and so on. In addition, differences in imagery between Landsat observation and local land cover in situ were visually cross-compared by means of the GPS toolbar under the ArcGIS Desktop. Both field photos and visual comparison greatly enrich our a priori knowledge about the remotely-sensed imaging characteristics of land use and cover types in MSEA.

Training Sampling and Swidden Landscape-Detecting Algorithms
Compared to other land cover types, such as natural or secondary forests and waterbody, swidden landscape simultaneously has features of low vegetation cover (biomass) and low moisture content. The features can be quantified by NDVI, NDMI, NBR and SAVI, respectively. To quantitatively reveal the differences of VIs among various land cover types, 2388 points of interest (POIs) for land cover types, including built-up land (238)  8 teak), shrubland, grassland, swidden agriculture, croplands (e.g., paddy rice and cassava), villages and towns, and so on. In addition, differences in imagery between Landsat observation and local land cover in situ were visually cross-compared by means of the GPS toolbar under the ArcGIS Desktop. Both field photos and visual comparison greatly enrich our a priori knowledge about the remotelysensed imaging characteristics of land use and cover types in MSEA.

Training Sampling and Swidden Landscape-Detecting Algorithms
Compared to other land cover types, such as natural or secondary forests and waterbody, swidden landscape simultaneously has features of low vegetation cover (biomass) and low moisture content. The features can be quantified by NDVI, NDMI, NBR and SAVI, respectively. To quantitatively reveal the differences of VIs among various land cover types, 2388 points of interest  Considering the dynamics of swidden landscape, especially the fallows with distinct recovering processes, remote sensing monitoring of swidden landscape via a single threshold method may not be appropriate. Therefore, a multi-step threshold method (MST) was proposed to delineate the fallows of the swidden landscape. The procedures for monitoring swidden landscape are shown in Figure 5. To begin with, a combination of VI criteria (Table 2), including NDVI greater than 0.10, Considering the dynamics of swidden landscape, especially the fallows with distinct recovering processes, remote sensing monitoring of swidden landscape via a single threshold method may not be appropriate. Therefore, a multi-step threshold method (MST) was proposed to delineate the fallows of the swidden landscape. The procedures for monitoring swidden landscape are shown in Figure 5. To begin with, a combination of VI criteria (Table 2), including NDVI greater than 0.10, NDMI, NBR and SAVI simultaneously less than or equal to five threshold values (0.20, 0.25, 0.30, 0.35 and 0.40), respectively, were applied to discriminate non-vegetation (e.g., newly-burned fields, swidden fields, unplanted farmlands and built-up) and little vegetation cover (i.e., fallow fields at various rehabilitating stages) from well-grown vegetation (e.g., natural or secondary forest, tree plantations and growing crops). It is noted that swidden landscape and non-swidden, including agricultural residue burning, unplanted farmlands and built-up land, were extracted in this step. Particularly, the first level of VI criteria (NDVI ě 0. NDMI, NBR and SAVI simultaneously less than or equal to five threshold values (0.20, 0.25, 0.30, 0.35 and 0.40), respectively, were applied to discriminate non-vegetation (e.g., newly-burned fields, swidden fields, unplanted farmlands and built-up) and little vegetation cover (i.e., fallow fields at various rehabilitating stages) from well-grown vegetation (e.g., natural or secondary forest, tree plantations and growing crops). It is noted that swidden landscape and non-swidden, including agricultural residue burning, unplanted farmlands and built-up land, were extracted in this step. Particularly, the first level of VI criteria (NDVI ≥ 0.10, NDMI, NBR and SAVI ≤ 0.20) was used to extract the distribution of newly-cleared swiddens. The other four levels of VI criteria (NDVI ≥ 0.10, NDMI, NBR and SAVI ≤ 0.25, 0.30, 0.35 or 0.40, respectively) represented the relevant different-aged fallows.   Next, three masks of mountainous area (LER300, LER400 and LER500) in MSEA were utilized to map the swidden landscape in the mountainous region and remove the crop residue burning pixels of farmlands. In addition, the topographical masks obtained from elevation, slope and LER contribute to excluding the unplanted farmlands and built-up land. In general, local inhabitants are accustomed to settle down in the regions with low elevation (below 600 m) and low slopes (e.g., 5˝). Meanwhile, any topographical limit factor of high elevations, large LER values and steep slopes could sharply reduce the agricultural productivity and enhance the production cost. However, there is an exception for the rubber plantations in MSEA. Rubber trees are usually cultivated at an elevation of 600-1000 m with large slopes [52]. For the mature rubber plantations, they put forth new leaves in March with high NDVI [53]. For the young rubber plantations (less than three years), the canopy is open and may be misclassified as swidden landscape. However, it is should be noted that new cultivated rubber plantations were typically transformed from natural or secondary forests through slash-and-burn practice several years ago.

Validation
Accuracy assessment of Landsat-derived maps of swidden agriculture in the tropics is a challenging task, as swidden agriculture is normally a small-sized and household-based farming practice, typically scattered in the uplands with certain slope gradients, which leads to access restriction and exhibits rotational and dynamic features. In spite of this, we attempted to validate and evaluate the resultant maps based on GE and ground surveys. Due to access limitation of swidden plots, field photos were often not taken in situ within the plots, but at a short distance. In the course of algorithm development, we noticed that unplanted permanent farmland could produce noise with respect to swidden identification. Therefore, these field photos were used as ground truth data to select regions of interest (ROIs). Specifically, we linked the OLI-derived maps of swidden landscape with GE and randomly selected 373 ROIs (150,972 pixels) for swidden agriculture and 373 ROIs (52,639 pixels) for permanent farmland as the ground reference for accuracy assessment via a confusion matrix mainly based on the finer spatial resolution imagery provided by Digital Globe and CNES/Astrium acquired from 2012-2014.

Results
Latest Maps of Swidden and Different-Aged Fallow Landscape in 2014/2015. The multi-step threshold (MST) method with various criteria of VIs was applied to detect and delineate swidden landscape consisting of newly-opened swiddens and multi-year fallows, typically less than seven years, in combination with a mountainous area mask. According to some recent studies, swidden fields with fallow periods up to 7-10 years usually regenerate to a normal level as natural primary forests in the aspects of soil fertility and vegetation coverage [54][55][56]. Swidden agriculture is an agroforestry system characterized by rotational and dynamic features in land use and land cover. It is even more complex for fallows at different regeneration stages. This brings about many challenges for monitoring swidden agriculture landscape with traditional one-size-fits-all approaches. In this study, the MST method involved two facets for setting multi-step threshold values. One is to set multi-stages for swidden and fallow field based on the VIs (NDMI, NBR and SAVI) with a step of 0.05 from 0.20-0.40 and; the other is to define various mountainous regions with the LER values (300, 400 and 500 m).
To reduce the effect of agricultural residue burning, topographical factors, including elevation, slope and LER, are foremost for the upland swidden monitoring. According to the UNEP-WCMC's mountain typology scheme, elevations above 300 m and LER over 300 m are two of three basic indicators (including slope) for defining mountains globally. However, areas with elevations lower than 600 m are preferentially used for living and permanent agriculture [53]. In addition, LER was used to discriminate the swiddening system in regions with larger relief, even though elevations may not be particularly high from permanent agriculture in regions with conversely high elevation but lower relief [45]. Therefore, we used the LER300-and LER500-derived mountainous masks to optimistically and conservatively estimate the area of swidden agriculture, respectively. A reasonable estimate was obtained based on the LER400-derived mountainous mask, as the error ranges between 9% and 16%. Table 3 shows the estimated area of swidden agriculture and the corresponding proportion in MMSEA. It should be noted that the area of swidden agriculture was accumulated from newly-opened swidden to the fourth stage of swidden fallow. Spatially, Figure 6 displays the distribution patterns of swidden landscape under different scenarios of LER derived from the available Landsat OLI imagery in MMSEA during 2013-2015. Figure 7 displays the spatial differences between the reasonable estimation and the optimistic and conservative estimates. Table 3. Land area and its proportion of swidden agriculture landscape with the MST method in MMSEA.

Area (km 2 ) Proportion (%) of the Total Land Area in MSEA
Range LER300 LER400 LER500 LER300 LER400 LER500 11 SAVI) with a step of 0.05 from 0.20-0.40 and; the other is to define various mountainous regions with the LER values (300, 400 and 500 m).
To reduce the effect of agricultural residue burning, topographical factors, including elevation, slope and LER, are foremost for the upland swidden monitoring. According to the UNEP-WCMC's mountain typology scheme, elevations above 300 m and LER over 300 m are two of three basic indicators (including slope) for defining mountains globally. However, areas with elevations lower than 600 m are preferentially used for living and permanent agriculture [53]. In addition, LER was used to discriminate the swiddening system in regions with larger relief, even though elevations may not be particularly high from permanent agriculture in regions with conversely high elevation but lower relief [45]. Therefore, we used the LER300-and LER500-derived mountainous masks to optimistically and conservatively estimate the area of swidden agriculture, respectively. A reasonable estimate was obtained based on the LER400-derived mountainous mask, as the error ranges between 9% and 16%. Table 3 shows the estimated area of swidden agriculture and the corresponding proportion in MMSEA. It should be noted that the area of swidden agriculture was accumulated from newly-opened swidden to the fourth stage of swidden fallow. Spatially, Figure 6 displays the distribution patterns of swidden landscape under different scenarios of LER derived from the available Landsat OLI imagery in MMSEA during 2013-2015. Figure 7 displays the spatial differences between the reasonable estimation and the optimistic and conservative estimates. Table 3. Land area and its proportion of swidden agriculture landscape with the MST method in MMSEA.

Area (km 2 ) Proportion (%) of the Total Land Area in MSEA Range
LER300 LER400 LER500 LER300 LER400 LER500 Newly opened swidden 32  The resultant map (LER300 derived) of swidden agriculture landscape was validated according to the confusion matrix based on the ground reference points ( Table 4). The classification results are encouraging with the overall accuracy up to 86.9% and the Kappa coefficient of 0.864. By contrast, the accuracies of LER400-and LER500-derived maps were higher. The estimated producer's accuracies (±95% confidence interval) for swidden agriculture landscape and permanent farmland are 0.839 ± 0.001 and 0.955 ± 0.002, respectively. The estimated user's accuracies for swidden agriculture landscape and permanent farmland are 0.982 ± 0.002 and 0.674 ± 0.002, respectively. The estimated The resultant map (LER300 derived) of swidden agriculture landscape was validated according to the confusion matrix based on the ground reference points ( Table 4). The classification results are encouraging with the overall accuracy up to 86.9% and the Kappa coefficient of 0.864. By contrast, the accuracies of LER400-and LER500-derived maps were higher. The estimated producer's accuracies (˘95% confidence interval) for swidden agriculture landscape and permanent farmland are 0.839˘0.001 and 0.955˘0.002, respectively. The estimated user's accuracies for swidden agriculture landscape and permanent farmland are 0.982˘0.002 and 0.674˘0.002, respectively. The estimated overall accuracy is 0.869˘0.001. The area estimate with a 95% confidence interval of the swidden agriculture landscape class is 11,610.4˘13.0 ha and 6741.6˘28.0 ha for permanent farmland.
Remote Sens. 2016, 8,44 overall accuracy is 0.869 ± 0.001. The area estimate with a 95% confidence interval of the swidden agriculture landscape class is 11,610.4 ± 13.0 ha and 6741.6 ± 28.0 ha for permanent farmland. Maps showing the differences of between the reasonable estimation (LER400 derived) and the optimistic and conservative estimates, namely LER300 derived and LER500 derived, respectively.

Discussion
Our swidden agriculture landscape detecting approach may have a common feature that differentiates low to no vegetation/moisture from medium to high vegetation/moisture with the landscape mosaic method by Hurni and his colleagues [20,57]. However, the algorithms proposed in this study progresses beyond other approaches that just provide an overall pattern of swidden agriculture landscape, but delineating swidden agriculture at different fallow stages in detail. One innovative point of the MST method should be the consideration of both newly-burned plots and different stage fallows. Our research work is an exploratory study, but proved to be an important attempt. In this study, however, a preliminary conceptive classification system of different-phased fallow was qualitatively applied to describe the fallow vegetation recovery. Further research needs to compare and determine the value range of VIs for different year-old fallow with year series data.
Several factors could potentially affect the detection and mapping of swidden agriculture landscape using the MST method, which involves the multi-step of VI thresholds and topographical relief thresholds. The first factor is related to the difference of acquisition date of the Landsat OLI imagery spanning from mid-January-May. The Landsat OLI sensor is scheduled for acquisition in the morning (nearly 10:00 a.m.) every 16 days, while the swiddening practice does not have a fixed date, although comparatively concentrated from 12:30 p.m.-2:30 p.m. during March and April [19]. As we used a simple dated method, swidden practice detecting depends on images acquired at a time point that could potentially omit some newly-opened swidden after the acquisition. Acquisition Maps showing the differences of between the reasonable estimation (LER400 derived) and the optimistic and conservative estimates, namely LER300 derived and LER500 derived, respectively.

Discussion
Our swidden agriculture landscape detecting approach may have a common feature that differentiates low to no vegetation/moisture from medium to high vegetation/moisture with the landscape mosaic method by Hurni and his colleagues [20,57]. However, the algorithms proposed in this study progresses beyond other approaches that just provide an overall pattern of swidden agriculture landscape, but delineating swidden agriculture at different fallow stages in detail. One innovative point of the MST method should be the consideration of both newly-burned plots and different stage fallows. Our research work is an exploratory study, but proved to be an important attempt. In this study, however, a preliminary conceptive classification system of different-phased fallow was qualitatively applied to describe the fallow vegetation recovery. Further research needs to compare and determine the value range of VIs for different year-old fallow with year series data.
Several factors could potentially affect the detection and mapping of swidden agriculture landscape using the MST method, which involves the multi-step of VI thresholds and topographical relief thresholds. The first factor is related to the difference of acquisition date of the Landsat OLI imagery spanning from mid-January-May. The Landsat OLI sensor is scheduled for acquisition in the morning (nearly 10:00 a.m.) every 16 days, while the swiddening practice does not have a fixed date, although comparatively concentrated from 12:30 p.m.-2:30 p.m. during March and April [19]. As we used a simple dated method, swidden practice detecting depends on images acquired at a time point that could potentially omit some newly-opened swidden after the acquisition. Acquisition before March or after April could weaken the spectral differences between swidden agriculture landscape and other land cover types. It should be noted that the time variations of acquisition only influence the accuracy of newly-opened swidden identification.
A second factor is the lower limit of elevation set for defining mountain regions. According to UNEP-WCMC, elevation above 300 m¨asl is the fundamental requirement for mountain definition. In this study, the parts of swidden agriculture located between 300 m and 600 m were not considered. In the past, this area was given highest priority for swidden agriculture due to the easier accessibility. However, cash crop (e.g., sugarcane and tea gardens) farming supported by the Alternative Development policy is increasingly occupying the frontiers of the uplands in MSEA. This phenomenon was relatively common in the border area of Myanmar, Laos and Yunnan Province, as we noticed during the field trip earlier this year. This treatment could fail to monitor the swidden farming in this area. The reason why we handled it in this manner was attributed to reducing the accuracy.
A third factor is the agriculture residue burning of upland permanent farming. Intuitively, it is very difficult to discriminate agriculture residue burning from swidden-related burning based on the Landsat spectral differences. Topographical factors, as we considered in this study, greatly reduce the classification error. However, the permanent croplands in the upland area could still bring disturbances to swidden system monitoring. Actually, permanent agriculture and built-up land are two factors that may possibly increase the difficulties for swidden agriculture extraction.
A fourth factor is the cloud and cloud shadow contamination. The issue of cloud and its shadow is a serious problem for remote sensing in the tropics. The dry season with less and occasional rain in MSEA helps with acquiring cloud-free imagery, which facilitates the formation of this study [25]. However, 55.8% of the total imagery was still partially contaminated by clouds and cloud shadows, with the mountainous area of 0.14% for cloud shadow and 0.53% for cloud, respectively. The pixels affected by cloud and cloud shadow were excluded for further analysis, which also may omit a small part of swidden agriculture.
In summary, the study has demonstrated the potential of the newly-launched Landsat 8 OLI imagery and the swidden agriculture landscape identifying algorithms for the large area mapping of swiddens and different-aged fallows. This is made possible due to the availability of Landsat-derived VIs covering the elements of vegetation, moisture and soil, which enable us to develop comprehensive and diversified algorithms. Our swidden agriculture landscape-detecting algorithms are developed based on the understanding of land cover changes, i.e., conversions from well-grown forests to newly-burned plots and gradually rejuvenating back to secondary forests. On the one hand, we focused on the detection of sharp declines in the four VIs during the burning or post-burn stage for newly-burned pixels and the delineation of the continuous restoration of different-aged fallows, on the other hand. The application of the swidden agriculture landscape algorithms in other tropical regions in insular Southeast Asia, South Asia, Latin America and Central Africa, where swidden agriculture (or a local term) dominates, especially in the REDD+ implementation countries, could provide a consistent dataset of the distribution, area and extent of swidden agriculture in the tropics. In addition, the algorithms can also be used to monitor the long-term change patterns of swidden agriculture landscape with Landsat 5/7/8 historical archive data in MMSEA ever the late 1980s. The long-term dataset of the distribution and extent of swidden agriculture landscape in the pan-tropical regions serves as a database for evaluating the relationships with infrastructure accessibility (typically road construction and settlements) and geographical (topography and drainage) factors, investigating the swidden usage intensity by analyzing the swidden cycle, fallow length and vegetation restoration process. Swidden agriculture has been ecologically and economically sustainable with low levels of population density and market formation. As population grows and regional economic cooperation deepens, heavily driven by market profits, swidden agriculture has been experiencing rapid transitions to industrial agriculture and plantations with swidden cycle and fallow periods shortened. The ensuing transformations are widely recognized challenges to cultural identity, rural livelihoods, biological diversity and carbon sequestration. SEA is currently one of the most vigorous economic zones globally and also the forefront of dramatic land use and cover changes heavily driven by infrastructure construction and economic development; therefore, there is an urgent need for deepening and comprehensive investigations of the environmental changes and the impacts supported by multidisciplinary research work.

Couclusions
In this study, an up-to-date distribution and extent of swidden agriculture landscape at 30-m resolution was identified and delineated using a multi-step threshold method in MMSEA. To our knowledge, the current estimation and mapping of swidden agriculture landscape were the first attempt to monitor the centuries-old agroforestry practice itself at the regional scale. Besides, the 30-m map of the swidden agriculture landscape provides the most comprehensive and detailed information about its distribution, shape, size and extent. This is an obvious improvement compared to previous studies, as they either belittled swidden agriculture or just investigated it at a more local scale [57]. Swidden agriculture is a widespread and dominant land use category across the mountainous and hilly regions of South Asia and SEA. Large-scale analysis of swidden agriculture will contribute to understanding the related differences among various countries, which helps with implementing the REDD program. An effective and robust approach for detecting the scattered distribution of swidden and dynamic fallows becomes very necessary, especially for its feasibility and suitability across the whole MSEA. The MST method highlighted the unique rotational and dynamic features of the swidden system, namely opening new swiddens from secondary forests and leaving soil-depleted swidden fields to fallow for soil and vegetation regeneration. Consequently, this method introduced a combination of VIs (NDVI, NDMI, NBR and SAVI) involving the aspects of vegetation coverage, land surface moisture content and soil exposure to highlight the differences between the swidden agriculture landscape and other land cover categories, such as forests, waterbody and permanent crops. As there are always newly-cleared swiddens each year, we set the base criteria of VIs to extract these pixels. Then, an increment of 0.05 for NBR, NDMI and SAVI simultaneously were added to delineate the swidden fallow at four successive stages.