Mapping Plantations in Myanmar by Fusing Landsat-8, Sentinel-2 and Sentinel-1 Data along with Systematic Error Quantification

Abstract

Forests in Southeast Asia are experiencing some of the highest rates of deforestation and degradation in the world, with natural forest species being replaced by cropland and plantation monoculture. In this work, we have developed an innovative method to accurately map rubber and palm oil plantations using fusion of Landsat-8, Sentinel 1 and 2. We applied cloud and shadow masking, bidirectional reflectance distribution function (BRDF), atmospheric and topographic corrections to the optical imagery and a speckle filter and harmonics for Synthetic Aperture Radar (SAR) data. In this workflow, we created yearly composites for all sensors and combined the data into a single composite. A series of covariates were calculated from optical bands and sampled using reference data of the land cover classes including surface water, forest, urban and built-up, cropland, rubber, palm oil and mangrove. This training dataset was used to create biophysical probability layers (primitives) for each class. These primitives were then used to create land cover and probability maps in a decision tree logic and Monte-Carlo simulations. Validation showed good overall accuracy (84%) for the years 2017 and 2018. Filtering for validation points with high error estimates improved the accuracy up to 91%. We demonstrated and concluded that error quantification is an essential step in land cover classification and land cover change detection. Our overall analysis supports and presents a path for improving present assessments for sustainable supply chain analyses and associated recommendations.

1. Introduction

Globally, up to 27% of forest loss has been attributed to permanent land use change for commodity production, and is currently estimated at 61% in southeast Asia [1]. Forests in the Greater Mekong region (Thailand, Vietnam, Lao PDR, Myanmar and Cambodia) are experiencing some of the highest rates of deforestation and degradation in the world [2,3,4]. The main drivers of land cover change include expansion of agriculture and plantation estates, extraction of natural resources [5], infrastructure development, and small and large-scale logging. The underlying drivers of land cover change include population and economic dynamics, often intensified by weak governance [6,7]. As a consequence, wildlife, hydrological and ecological functions, and local communities that rely on forest resources are experiencing significant negative impacts.

Natural rubber trees are the main source of hevea latex, used in a number of industrial products, most importantly automobile tires. Commercial rubber tree plantations are tapped from 5 years after planting until 25 or 30 years, after which their production declines and the plantation is no longer economically viable and the trees are removed and replaced with new seedlings, or converted to more profitable palm oil [8,9]. Mature rubber plantations are thus regarded as waste, although the wood has good properties for use in furniture manufacturing [10]. In Myanmar, there is potentially a vast amount of mature rubber wood which could thus be available for alternative use as a source of timber. However, efforts should be made to reduce the environmental and social impacts of this exploitation. While some capacity exists with regards to accessing these resources, land tenure, market access, policies and knowledge gaps are considered main obstacles for market entry and sustainable rubber wood use [11]. To act on opportunities and challenges for nature conservation, business operations, land tenure/rights, resource management and related policies within the country in relation to rubber wood, bamboo, rattan and potential other natural fiber products it is important to first understand the current resource availability. The inventory of available mature plantation wood, particularly relative to biodiversity assets, protected areas, and the extent, size, and stand age of rubber plantations is supporting improved characterization of the available rubber wood supply for sustainable management of these resources. This information is key for the development of responsible natural fiber and timber supply chains for these commodities which benefit local stakeholders. This information can help to identify high and low risk business opportunities for alternative timber sources in Myanmar in order to have reduced impacts on biodiversity and livelihoods.

A number of studies have presented various methods for discerning planted rubber from other natural and plantation forests such as palm oil in Asia, notably in China and Thailand [12,13,14,15,16,17]. It has been determined that an essential factor for accurate mapping of rubber plantations in SE Asian ecosystems is the availability of data at two crucial phenological periods: defoliation (leaf off) and new leaf emergence (leaf on), which distinguishes deciduous rubber trees from other vegetation. Many approaches consist of using optical satellite data [6,17,18,19,20,21] to detect rubber plantations. However, due to the persistent cloud cover in tropical regions, it is difficult to acquire sufficient high resolution satellite imagery which in turn compromises spatial resolution, as coarser satellites such as MODIS are frequently used. A combination of Synthetic Aperture Radar (SAR) [22,23,24] has shown promising results since SAR can penetrate clouds, haze, and dust. As such, satellites such as Sentinel-1 and Sentinel-2 are a new data rich source for land cover mapping purposes and enable land cover mapping on a higher spatial resolution.

The new generation of satellites such as the Sentinel constellation offer exciting opportunities to monitor land cover (changes) on a high spatial resolution. Fusion products typically include more observations within a given period [25]. However, different sensors vary in spectral, spatial and temporal resolution [26]. With earth observations available from a variety of satellite sensors, data fusion techniques are being increasingly developed and utilized for various applications. The merging of data achieved by such data fusion approaches overcomes limitations from a single sensor, which is focused on either passive or active remote sensing. Hence, both single sensor spatial and temporal drawbacks can be overcome or partially overcome by data fusion techniques. With the development and testing of data fusion techniques, this is a rapidly developing area in remote sensing theory with applications ranging from water quality management [27] to forest mapping [28].

There is a large variety of data fusion techniques in which different types of sensors are combined. The Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) [29] and the Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model (ESTARFM) [30] are classic examples. These techniques involve combining lower spatial resolution optical imagery such as MODIS with higher spatial imagery such as Landsat in order to overcome the coarser temporal resolution in the latter sensor. This is achieved via the assumption that MODIS and Landsat observations are correlated with each other. Data fusion techniques involving optical and SAR data are also gaining acceptance with recent applications to better forest mapping [31] and surface water detection [32].

We present the following assessment in the context of the Regional Land Cover Monitoring System (RLCMS) developed by SERVIR-Mekong, which is a joint USAID and NASA collaborative project. It is aimed at providing support to dedicated development and sustainable landscape projects in the Mekong region. The RLCMS leverages state-of-the-art cloud computing technologies and includes an explicit quantification of accuracy and uncertainty. We describe these methods in detail in the following sections. Innovations of this study include the data fusion of different sensors and systematic error quantification to reduce errors. All calculations are done using a cloud-based remote sensing approach.

2. Materials and Methods

2.1. Study Region

The region of study is the Dawna Tenasserim Landscape (DTL) (Figure 1) which straddles the border between Myanmar and Thailand and covers the Kayin, Mon, and Tanintharyi states along the southeast coast of Myanamar. The area stretches from 8.82 to 19.5 degrees northern latitude. It has a tropical climate with significant rainfall most months and a short dry season. The area includes the Tenasserim Hills, where the Myinmoletkat Taung is the highest point with 2072 m. The area harbours one of the last large intact forest landscapes in the region, which hosts a large number of ethnic groups, and a wide variety of endangered wildlife such as elephants and tigers [33]. The region is experiencing numerous threats including rapid increase in agricultural land use and associated deforestation, notably for rubber and palm oil, road and infrastructure development, and logging resulting in habitat degradation and fragmentation [34].

2.2. Typology

This study focuses on mapping rubber and palm oil plantations, but also includes other existing dominant land cover types in the area. Classification thus requires the definition of class boundaries, which is clear, precise, possibly quantitative, and based on objective criteria (FAO; [35]). The land cover typologies used in this study are described below and visualized in Figure 2:

• Cropland: land covered by crops and then followed by harvest and a bare soil period [36]. Examples include cereals, oils seeds, rice, vegetables, root crops and forages. This excludes orchards, forest croplands, and forest plantations.
• Forest: Land spanning more than 0.5 hectares with trees higher than 5 m and a canopy cover of more than 10 percent, or trees able to reach these thresholds in situ. It does not include land that is predominantly under agricultural or urban land use [37].
• Mangrove: coastal sediment habitats with more than 10% woody vegetation canopy cover and the majority of cover is higher than 2 m [36].
• Palm oil: The land area designated primarily for production of palm oil. Palm oil grows in tropical climates within 10 degrees of the equator and high rainfall (minimum 1600 mm/year) [38]. Palm oil is a productive crop, it is planted as mono-culture and its expansion comes at the expense of tropical forests [39].
• Rubber plantation: Forest area with rubber tree plantations. Rubber plantation can be inter-crops with other plants, and has about 30–40 years of rotation.
• Surface water: Open water larger than 10 m by 10 m and open to the sky, including fresh and saltwater
• Urban & built up: Cultural lands covered by buildings, roads, and other built structures.

2.3. Methods Overview

Optical imagery from Landsat 8, Sentinel 2, and SAR images from Sentinel-1 were combined into a single stack of images. This image stack was used as predictors in the classification process. The training sample was created by assigning the available reference data the coincident image values. A random forest algorithm was applied to the training sample and then used to calculate biophysical probability layers (called primitives from hereafter). Primitives were used in a decision tree to create the final assemblage with plantations. During this process we used a Monte-Carlo simulation to create associated uncertainty layers. We discuss the specifics in more detail in the following sections. The general methodology is outlined in Figure 3. The links to the source code can be found in the Appendix A.

This study was conducted in Google Earth Engine (GEE), a cloud-based computing environment that includes access to the full archive of Landsat and Sentinel imagery [40]. GEE combines a large data archive of satellite imagery with a computational platform. The platform enables scientists to conduct research on environmental issues on a variety of spatial and temporal scales [41,42,43,44]).

For this study, we used the USGS Landsat 8 surface reflectance product. This product contains atmospherically corrected, orthorectified surface reflectance data. The images have been atmospherically corrected using the Landsat Surface Reflectance Code (LaSRC) [45,46,47] and also contains the data produced by CFMASK [48]. Landat 8 has a spatial resolution of 30 m. Images with more than 40% cloud cover were excluded from the analysis due to issues with haze.

The Copernicus Sentinel-2 mission provides high spatial resolution images. We used the Sentinel-2 image collection in Google Earth Engine which contains spectral bands representing Top Of Atmosphere (TOA) reflectance from the Sentinel 2a and 2b satellites. The spatial resolution for Sentinel-2 varies for the different bands. The blue, green, red and near-infrared bands have a resolution of 10 m, the red-edge and shortwave-infrared bands 20 m, and all others 60 m.

The Copernicus Sentinel-1 mission provides synthetic aperture radar images with high temporal and spatial resolution. We also used the Sentinel-1 data from the GEE archive. Google Earth Engine makes scenes available after their ingestion team has applied thermal noise removal, radiometric calibration, and terrain correction using the Sentinel-1 Toolbox processing algorithms. Pre-processing includes thermal noise removal, radiometric calibration, and terrain correction. We used the available polarization bands, which were the VV and VH dual bands. A filter was applied to de-speckle the image [49]. Sentinel-1 has a spatial resolution of 10 m.

For the optical imagery, composites were created for the rainy and dry season. Seasons were defined based on timing of rubber tree leaf senescence and emergence. The rainy season was defined as the period from May until December, the other months were used in the composites of the dry season. A yearly composite was used for the SAR imagery. Composites included a medoid, 20th and 80th percentile from the annual image collection.

Processing steps for Sentinel-2 included shadow and cloud removal and atmospheric correction. BRDF correction and topographic correction were applied to Landsat 8. The atmospheric correction was done by applying the 6S radiative transfer model. The original Second Simulation of the Satellite Signal in the Solar Spectrum (6S) model was developed by [50] and translated into Python by [51]. We used the Sentinel-2 6S implementation for the Google Earth Engine to atmospherically correct the images.

2.4. Cloud Shadow Masking

Cloud shadow removal is an essential step because of the negative influence cloud shadow can have on data analysis [52]. We used the Temporal Dark Outlier Mask (TDOM) algorithm [53]. This algorithm identifies pixels that are dark in the infrared bands but are found to not always be dark in past and/or future observations. It detects statistical outliers with respect to the sum of the infrared bands. Next, dark pixels are identified by the sum of the infrared bands (NIR, SWIR1, and SWIR2). The pixel quality attributes generated from the CFMASK algorithm (pixel-qa band) was also used for Landsat-8 shadow masking [48].

Cloud removal is another essential step in optical remote sensing. Clouds were removed using the quality bands for Landsat-8 and Sentinel-2. Furthermore, the Google Earth Engine cloudScore algorithm was used. The algorithm uses the spectral and thermal properties of clouds to identify and remove pixels with cloud cover from the imagery. The algorithm finds bright and cold pixels and compares them to the spectral properties of snow. The Normalized Difference Snow Index (NDSI) is calculated to prevent snow from being masked. The algorithm uses the visible, near-infrared, and shortwave infrared for a scaled cloud-score and then takes the minimum.

2.5. BRDF Correction

In the case of anisotropic surfaces, surface luminance depends on the properties of the surface, the location and type of light and the conditions of observation [54,55,56]. Nicodemus [57] described the characteristics of the light reflection for a specific surface by the BRDF (sr${}^{-1}$). The nadir view angles cause directional reflection on the surface which were described by [58,59,60]. We applied the BRDF correction to all images in the image collection to account for these effects.

2.6. Topographic Correction

The slope, aspect and elevation can cause variations in reflectance for similar features with different terrain positions [61,62,63]. Topographic correction is the process to account for these effects. The Modified Sun-Canopy-Sensor Topographic Correction method as described by Soenen et al. [64] was applied to account for these effects. The method uses the sun-canopy-sensor (SCS) [65] with a semi-empirical moderator (C) to account for diffuse radiation [66,67,68]. The ALOS Global Digital Surface Model was used [69,70].

2.7. Representing Vegetation Phenology Using Sentinel-1

The harmonic model represents the characterization of the temporal variability in a data-series using harmonic components [71]. Harmonic trend analysis [72] was applied to the Sentinel-1 time series as Sentinel-1 offers a complete time-series which is not affected by clouds. We applied one harmonic term (shown in Equation (1)), which represents one seasonal cycle per year [73]. Equation (1) has a linear component (${\beta }_0$ and ${\beta }_1$) and two harmonic coefficients (with ${\beta }_2$ and ${\beta }_3$). The phase and amplitude estimates at the pixel level were added as bands to the image stack.

$$\begin{array}{c}\hfill SAR={\beta }_0+{\beta }_{1}t+{\beta }_{2}cos\left(2\pi t\right)+{\beta }_{3}sin\left(2\pi t\right)\end{array}$$

(1)

2.8. Processing

The seasonal Landsat-8, Sentinel-1 and Sentinel-2 composites were combined into a single image, where the different images are represented as bands. A series of indices were calculated from the optical satellites. We applied a random forest classifier using the probability mode on these images to create probability maps for each land cover class. These maps were then assembled into a final land cover map using a decision tree logic and Monte-Carlo simulations to estimate uncertainty.

Covariates

We refer to covariates as the composite band values and derived indices that are used as predictors in the random forest algorithm.

Table 1. Normalized difference metrics (ND; Equation (2)) were calculated for every Landsat-8 and Sentinel-2 image. This table provides an overview of the band combinations.

Blue	Green	Red	nir	swir1
green	red	swir1	red	swir2
red	nir	swir2	swir1
nir	swir1		swir2
swir1	swir2
swir2

provides an overview of the band combinations used in the normalized difference (Equation (2)) calculations. For some combinations there are more common names such as Normalized Difference Water Index (NDWI) [74]), Normalized Burn Ratio (NBR) [75]), Normalized Difference Snow Index (NDSI) [76]) and Normalized Difference Vegetation Index (NDVI) [77].

$$ND(band1,band2)=\frac{band1-band2}{band1+band2}$$

(2)

Two ratio (R) bands were included, calculated by dividing one band by another. This was done for the SWIR1 and NIR bands and the red and SWIR1. The Enhanced Vegetation Index (EVI, Equation (3)) [78] was also included, as well as the soil-adjusted vegetation index (SAVI, Equation (4) using L = 0.5) [79]. The Index-based Built-Up Index (IBI) [80] was calculated from Equation (5).

$$EVI=2.5\ast \frac{NIR-red}{NIR+6\ast red-7.5\ast blue+1}$$

(3)

$$SAVI=\frac{(1+L)(NIR-red}{NIR+red+L}$$

(4)

$$IBI=\frac{NIR}{(NIR+red)}+\frac{green}{(green+swir1)}$$

(5)

The total number of bands in the final image stack was 234. We sampled the complete stack of covariates for each land cover class and evaluated the importance of the covariates in R [81,82,83]. A smaller number of covariates reduces the computational expense and eliminates noise. Whereas the bands have different spatial resolutions, sampling was done on a 10 m spatial resolution, the smallest spatial resolution of all bands.

2.9. Reference Data

Reference data was collected from both field data and through photo-interpretation of satellite imagery. High resolution data from RapidEye and PlanetScope have a spatial resolution of 5 and 3 m respectively and were ingested into the GEE platform. Opportunistic sampling was done using the high resolution data and Landsat-8 and Sentinel-2 composites. Reference data for forest and cropland were obtained by applying a stratified sample on the SERVIR-Mekong regional land cover product [84].

Table 2. Tally of the number of reference data points collected for each class.

Land Cover Class	No. Training Points
Cropland	1758
Forest	1422
Mangroves	873
Palm oil	1467
Rubber	1802
Water	499
Urban and built up	993

provides an overview of the total number of points that were collected in each class.

2.9.1. Training Data

The reference data was used to create a training data set for the random forest classifier. The data of each class was combined with a random sample of all classes to create a sample representing both class presence and absence (other land covers).

2.9.2. Primitives

The table with training data was used in a random forest model in R [81,82,83] to select the most important covariates. Metrics with a correlation in excess of 0.90 with other variables were removed to reduce problems associated with highly collinear predictor variables. Various methods are available to select a parsimonious set of metrics to use as predictors in a model, including dimensional reduction of the data, such as a principle component analysis or canonical correlation analysis. We used the information on the variable importance measures to select the covariate list. The random forest classifier was then applied in GEE using the selected most important covariates and the training data. The classifier was trained with 100 trees in probability mode. The classifier was then applied to each class (class versus other), resulting in the probability map for the class.

2.9.3. Assembly Logic and Monte-Carlo Simulations

A decision tree and Monte-Carlo simulation were used to create the final land cover assemblage. The decision tree was used to set the order and thresholds used to combine each of the primitives together into one final land cover map. The decision tree was run 100 times with a Monte-Carlo simulation process [85]. The Monte-Carlo process entailed adding a grid with random numbers to each of the primitives. It produces a final land cover map, which is the mode of the 100 simulations and a probability map which is the count of the mode divided by the total number of model runs.

2.10. Validation

We placed a stratified random sample in the final land cover product. The strata used were pixels classified as rubber plantation, palm oil, and an aggregate other class. 1500 plots were located in each of the three strata. The land cover at each plot was analyzed in Collect Earth Online using available imagery from Digital Globe, Planet, and Bing maps [86]. From the 1500 plots, 47 were flagged as bad for various reasons such as no or bad imagery. We then calculated a confusion matrix, overall accuracy, and disagreement [87,88,89]. We also plotted how overall accuracy changes when pixels with a low level of uncertainty are iteratively filtered out. We also plot the distribution of the probabilities as boxplots by class type.

3. Results

We created seven primitive maps, including cropland, forest, mangrove, palm oil, rubber plantation, urban, and water. It is important to have a clear separation in the probability distribution between the presence and absence of each class primitive layer. The boxplots of probability values by land cover type for each of the seven primitive presence/absence maps is shown in Figure 4. There is distinct separation for most class primitive layers. However, there is notable overlap between urban, water, cropland and forest. There is also overlap between mangrove and water.

To specify the probability distribution for the Monte Carlo simulation we aggregated all absence classes into one sample and calculated the 2.5 and 97.5 percentiles. The same process was repeated for the presence sample. These are represented in Figure 4 as dashed lines. We chose these thresholds rather than the minimum and maximum to exclude noisy data points. These thresholds are far apart for the cropland, forest, palm oil and rubber plantation primitives. They are much closer for the mangrove, urban, and water primitives. In the Monte-Carlo simulation step, we used the mean of the two percentiles as the threshold for the decision tree while the distance from the mean was used as the maximum variation in the simulations.

The data from Figure 4 were used in a decision tree classifier. Mangrove was placed on top followed by cropland, palm oil, forest, rubber, water and urban. The resulting land cover and uncertainty map from the Monte-Carlo simulation are shown in Figure 5. It can be seen that the hilly inland areas are classified as forests with extensive cropland areas in the north. Mangroves are mainly found in the Tanintharyi Region. Palm oil plantations are found throughout the region, but more concentrated along the coastal areas. The uncertainty map shows the spatial distribution of uncertainties throughout the region.

Validation was done using a confusion matrix. The confusion matrix of the 2017 and 2018 land cover maps can be found in

Table 3. Confusion matrix for 2017 and 2018

	Surface Water	Forest	Urban and Built Up	Cropland	Rubber	Palm Oil	Mangrove
Surface water	40	0	6	0	0	0	1
Forest	2	1531	14	18	66	17	0
Urban and Built up	0	0	17	1	0	0	0
Cropland	9	16	2	476	14	6	0
Rubber	4	68	2	84	749	26	0
Palm oil	2	89	0	2	27	262	1
Mangrove	11	79	0	22	6	4	163

. The overall accuracy was 0.84, with an overall quantity disagreement of 0.07 and an overall allocation disagreement of 0.09, following the definitions of Pontius Jr and Millones [88]. Highest accuracy values were found for mangroves, rubber, and forest types. Lower accuracy values for water and Cropland. Surface water was mostly confused with mangroves; cropland was confused with rubber. An uncertainty map was used to filter the validation data. Validation points with lower accuracy values were masked out by varying the threshold. The results are shown in Figure 6. With this filtering, it can be seen that the accuracy increases from 0.84 up to 0.91. The increase in accuracy starts at around 50%.

The probabilities of the land cover map are shown in Figure 7a. It can be seen that the overall confidence of the map is generally high (above 75%). Figure 7b shows the probabilities for the different land cover classes. It was found that other, palm oil and rubber plantations have the largest variability. The probabilities of classes that change from 2017 to 2018 were also calculated. The probability distribution of the combined classes is shown in Figure 7c and specified according to class in Figure 7d. It can be seen that the variability increases dramatically when zooming in on those classes. Also, it is noted that cropland and urban have high probabilities in both situations.

Figure 8 illustrates a close up example of the land cover class transitions from 2017–2018. For each transition, we selected all points that undergo change and map to which class they change (as a percentage) and the probability distribution of that specific class change. For example, cropland sees most changes to rubber plantations, forest and other. It was found that changes to other land-use is generally a large portion of the changes, whereas probabilities are variable. For mangroves, most transitions are to forest.

4. Discussion

We have presented a method to fuse optical and SAR imagery in the context of nature conservation in Myanmar. We have demonstrated that a combination of Sentinel-2, Landsat-8 and Sentinel-1 provides a sufficiently dense spatio-temporal data series, even in regions with persistent cloud cover, to separate between different land cover and plantations types. Due to the nascence of best pre-processing workflow and algorithms for Sentinel 2 data, the processing steps of the two optical satellites–sentinel 2 and Landat 8–differed. BRDF and topographic corrections were applied to the Landsat 8 time series, but were not available for the sentinel 2 images. Work of [90] and [91] present a method to apply the BRDF correction to Sentinel 2, but they used a different model for atmospheric correction. As such, values presented in their work could not be applied to our study. Similarly, we used a conservative approach with regards to applying topographic correction. Initiatives such as the Harmonized Landsat and Sentinel 2 surface reflectance data set [92] are very exciting in this regard, as they provide a consistent integration of Landsat and Sentinel 2 images which has been extensively tested and include consistent methods to account for BRDF effects [93].

Error quantification is a key quality check that provides accountability to all stakeholders. We have demonstrated that the use of a Monte-Carlo simulation with a decision tree logic is a very useful tool to map uncertainty. The accuracy values reported in this study range from 84% up to 91% when taking error quantification into account. These numbers are in the same range as e.g., Dibs et al. [94] and Dong et al. [95]. Information on the potential error of a pixel can be used for additional reference data collection or might provide valuable information on the spatial distribution of land cover classes that overlap. Land cover type specific information on uncertainty and overlap between classes is also useful in the decision tree logic, where classes with higher uncertainty or classes with notable overlap can be re-arranged according the purpose and priority of the final land cover map.

Time-series analysis along with error quantification and stand age estimation [96] might provide other valuable information. An example is give in Figure 8 where land cover change is mapped with the associated uncertainty of the change. However, the presented method uses a single composite for every year whereas more dedicated time-series analysis methods such as Landtrendr [97] and Breaks For Additive Season and Trend (BFAST) [98] use complete time-series to detect changes. Other methods such as intensity analysis can be powerful to analyze changes in land categories over time [89,99]. Future work can cross-walk these different methods and improve future land cover and land cover change products.

The greatest amount of rubber plantations was found in Tanintharyi state (approximately 468,000 ha), followed Kayin state (approximately 344,000 ha) and then Mon (approximately 246,000 ha) for 2018. These estimates differ from the numbers reported by the Land Management and Statistics Department (MOALI) for the years 2015–2016. They reported 138,828 ha, 107,978 ha, and 198,741 ha respectively. Similarly, the numbers vary from a study by Connette et al. [4] who estimated rubber occupied 127,500 ha in the Tanintharyi state in 2016. The same study reported 136,500 ha of palm oil, where we found 125,253 ha for the year 2018. In contrast, another study from Torbick et al. [100] reported an area of 750,822 ha for palm oil and rubber plantations combined in Tanintharyi state, where we found 594,000 ha. However, more importantly, nearly half of the plantations in Dawei and Myeik districts in Tanintharyi state are located within key biodiversity areas, protected areas or wildlife corridors, highlighting the risk associated with extraction of rubber wood in these areas. This analysis provides improved assessments for sustainable supply chain analyses and recommendations.

In this study, we have demonstrated how sensor fusion and explicit error quantification can be beneficial for earth science applications using a cloud-based platform. Additional sensors or imagery from higher resolution resources might improve the results by better spatial and/or temporal representation to discern forest types [101]. Including additional statistical information such as elevation, slope, aspect, and orientation of slope. We also found ancillary remote sensing products such as the Joint Research Center (JRC) Global Surface Water product [102], global forest cover products [103] and data products on cropland extent [104] can be used to further improve results. Combining various classifiers based on a priori knowledge could further enhance the method [105]. Machine learning algorithms using neural network approaches are also evolving and are promising for integration into the methods presented here to potentially improve the results by mapping canopy species [106,107].

5. Conclusions

The study provides valuable insights into discriminating between natural and anthropogenic forest types. This is particularly applicable in regions where natural forest is quickly being replaced with plantations, and these conversions may not be detected with traditional remote sensing methods. Given this situation, we fused data from Landsat 8, Sentinel 2, and Sentinel 1 data to overcome the various limitations present in each individual sensor. In such land cover classification efforts, systematic error quantification is a key quality check and provides accountability to all stakeholders. Starting with fused data and primitives as biophysical probability layers, we achieved both final classification and system error quantification in the RLCMS workflow. This process resulted in an overall good accuracy which can be further improved if validation data with higher error estimates are used. We finally contend that these results are essential inputs to guide sustainable management of wood resources in Myanmar.