Assessment of forest restoration with multitemporal remote sensing imagery

Climate variability and man-made impacts have severely damaged forests around the world in recent years, which calls for an urgent need of restoration aiming toward long-term sustainability for the forest environment. This paper proposes a new three-level decision tree (TLDT) approach to map forest, shadowy, bare and low-vegetated lands sequentially by integrating three spectral indices. TLDT requires neither image normalization nor atmospheric correction, and improves on the other methods by introducing more levels of decision tree classification with inputs from the same multispectral imagery. This approach is validated by comparing the results obtained from aerial orthophotos (25 cm) that were acquired at approximately the same time in which the Formosa-2 images (8 m) were being taken. The overall accuracy is as high as 96.8% after excluding the deviations near the boundary of each class caused by the different resolutions. With TLDT, the effectiveness of forest restoration at 30 sites are assessed using all available multispectral Formosat-2 images acquired between 2005 and 2016. The distinction between natural regeneration and regrowth enhanced by restoration efforts were also made by using the existing dataset and TLDT developed in this research. This work supports the use of multitemporal remote sensing imagery as a reliable source of data for assessing the effectiveness of forest restoration on a regular basis. This work also serves as the basis for studying the global trend of forest restoration in the future.

In recent years, recurring natural disasters such as wildfire and landslides have severely damaged forest ecosystems, resulting in a loss of habitat for numerous species, significant soil erosion and changes in land cover 1 . This condition is aggravated by the current climate variability and man-made impacts including timber harvest and fire succession that are detrimental to the forest environment on an even more massive scale 2 . The current scale of deforestation all around the world calls for an urgent need to restore biodiversity and the ecological structure and functioning, aiming toward long-term sustainability for the forest environment 3 . Although many countries' governments and various environmental management sectors have launched initiatives on restoring forest ecosystems, they lack a systematic and synoptic view for monitoring the effects of forest restoration 3 . Thus, it had become increasingly important to assess the effectiveness of forest restoration and monitor the posttreatment abiotic and biotic characteristics of the landscape 1,3 . Without an effective monitoring system, there is neither sufficient information provided on the impact on the restoration efforts, nor is there any basis for further improvements 3 .
Taiwan is located in the center of the East-Asian island arc formed by the slow collision of the Asian continental plate and the Philippine plate. The high mountains, broken terrain and frequent earthquakes, together with the heavy rainfall during the rainy and typhoon seasons, results in a very high erosion rates in the world 4 . Consequently, more than 90% of the country's population lives in areas that are at a relatively high risk of typhoons, earthquakes and landslides 5 . For example, Typhoon Morakot brought an extreme precipitation of 2,777 mm in less than a week in August 2009 6 , and triggered enormous landslides that caused massive destructions to the landscapes. To stabilize those areas affected by landslides found in the mountainous region, the Forestry Bureau of Taiwan (FBT) launched initiatives and utilized various techniques on restoring forest. The traditional method for evaluating forest restoration is in situ site surveying twice per year, once in the summer www.nature.com/scientificreports www.nature.com/scientificreports/ Material study sites. A total of 511 sites of forest restoration in Taiwan have been accomplished so far by FBT. Their geographical locations are labeled as blue polygons in Fig. 1(a). Considering the cost and availability, the major source of satellite imagery for this work is from Formosat-2, which was operated by National Space Organization of Taiwan from 2004 to 2016. Therefore, the timing of restoration in the study site has to be later than 2004 so the full process of restoration can be covered by the span of Formosat-2 mission. In addition, the area of restoration should be large enough to ensure that the gradual changes can be captured by the 2-m resolution of Formosat-2 imagery. Taking these conditions into consideration, a total of 30 study sites are selected and grouped into five zones, as marked as star symbols and blue boxes in Fig. 1(b). The detailed descriptions of each site including the restoration year and engineering method are listed in Table 1. Note that the results of restoration effectiveness assessed by this work are also listed in Table 1 for comparison. Together with the overall accuracy of 15 test sites, these results will be explained in detail later.
Remote sensing imagery. Formosat-2 is the first satellite with a high-spatial-resolution (2 m) sensor placed in a daily revisit orbit, as well as the second satellite that is owned and operated by the National Space Organization (NSPO), Taiwan 18,19 . The remote sensing instrument (RSI) onboard Formosat-2 acquires 2-m resolution panchromatic images and 8-m resolution multispectral images in four multispectral bands (blue, green, red and near-infrared) over 24 km swath width in the nadir direction. The spectral bands definition and the spectral radiances at the entrance aperture, including the mean radiance and the saturation radiance can be referred to Table 2 of Liu et al. 19 . RSI has a field of regard of ±45 deg for along-track and cross-track viewing, and our study areas are covered by Strip Nos. 3, 4 and 5 of Orbit 1. During twelve years of operations and services from 2005 to 2016, Formosat-2 has acquired a total of 1,453 strips covered our study areas. All quick-look images (Strip Nos. 3, 4 and 5 of Orbit 1) are annotated with the dates of acquisition and given in the supplement A. After excluding those images with mostly clouds and haze, there are 42, 37, 31, 28 and 26 available images for zone 1 to 5 respectively, as listed in Table 2. Note that four seasons are indeed covered in the entire time series of images for every zone. But it is not the case for every zone to have images for every season in every year.
All available Formosat-2 imagery of each zone are pre-processed by F-2 AIPS 18 . F-2 AIPS is able to digest the Gerald format of the raw data, apply the basic radiometric and geometric correction, output the level-1A product, conduct the rigorous band-to-band co-registration 20 , automatic orthorectification 21 , multitemporal image geometrical registration 22 , multitemporal image radiometric normalization 23 , spectral summation intensity modulation pan-sharpening 20 , and the absolute radiometric calibration 19 . One true color image of zone 2 taken by Formosat-2 on 9 June 2015 is shown in Fig. 2 as an example, of which the regions of site No. 355, 378 and 385 are annotated as blue, red and green boxes respectively. The time series of co-registered and radiance-normalized imagery are further divided into smaller rectangles so that each study site is fully enclosed by a small rectangle. Note that two landslides near site No. 378 are selected and labeled as the natural restoration site (NRS) 1 (white box) and 2 (yellow box) in Fig. 2. These landslides were also triggered by Typhoon Morakot but they were not selected as restoration sites by FBT. Therefore, the changes found at NRS 1 and 2 are purely the effect of natural regeneration, which can help us to gain insight into the difference between natural regeneration and restoration.
The general source of ground truth is usually collected on the ground by conducting a considerable amount of in situ survey. Interpreting aerial colored images to obtain the ground truth information offers an alternative to this general approach 24 . For 15 of 30 study sites, the aerial orthophotos with 25 cm resolution and four spectral bands (blue, green, red, and near-infrared) were provided by FBT. Note that the position, attitude values, as www.nature.com/scientificreports www.nature.com/scientificreports/ well as digital image of the photograph at the moment of being taken were decoded by using Global Positioning System (GPS) and Inertial Measurement Unit (IMU). Together with the Digital Terrain Model (DTM), the center projected aerial photographs were rectified pixel by pixel into orthophotos. These orthophotos are resampled to 2 m resolution and co-registered with their corresponding Formosat-2 imagery, to serve as an alternative source of ground truth. These 15 sites are marked as purple star symbols in Fig. 1(b) and used as test sites to validate our approach.

Method
For detecting the changes of land cover, and in our case, assessing the effectiveness of forest restoration, the radiometric normalization of multitemporal satellite optical images of the same terrain is often necessary 25 . Absolute radiometric correction requires an atmospheric correction algorithm and the associated atmospheric properties at the time of image acquisition, but this information is either unavailable or difficult to obtain 26 , particularly for Formosat-2 imagery that has only four spectral bands with broad bandwidth. Therefore, we attempt to employ the CART approach that requires neither image normalization nor atmospheric correction to determine thresholds for disturbance or regrowth. However, we also found that the process of radiance normalization provides an appropriate way of examining the quality of every image.  25 . Using the image taken on 7 July 2012 as an example, the DN values of four spectral bands collected at all searched PIFs before (cross marks) and after (circle marks) the process of radiance normalization are shown in Fig. 3(a). This scatter plot demonstrates that a consistent and robust correlation can be established from these PIFs to meet the requirement of radiance normalization for change detection 28 . By contrast, the image taken  Fig. 1(b). * A: staking and wattling; B: seeds spreading; C: drainage of longitudinal and transverse; D: tree planting; E: Cover grass net. ** I: restored well apparently; II: restored in a slow pace; III: restored inefectively; IV: no assessment. on 15 January 2013 gives an example that fails to give a one-to-one relationship because the image was taken with a narrow dynamic range, as shown in Fig. 3(b). Therefore, the process of radiance normalization provides an appropriate way to examine the quality of every image. All scatter plots of the other 34 images are given in the supplement B, of which a total of 19 images with poor quality are excluded from the time series analysis, as listed in Table 3.

Radiance normalization.
three-level decision tree approach. Automatic classification of landslides from multispectral imagery has been progressing rapidly and a few practical approaches have been proposed in the past few years. There are some similarities and differences among these approaches. Both Yang et al. 29 and Liu 30 classified the bare/ non-vegetated land by using the normalized difference vegetation index (NDVI). www.nature.com/scientificreports www.nature.com/scientificreports/ for its superiority to differentiate vegetation. Ma et al. 31 took the same idea but used the inverse normalized dif- Since the scattering of visible and near-infrared light are different, the partially-shaded topography regions are often misinterpreted as landslides for NDVI or NDVI * are distorted in these regions. Therefore, the normalized difference soil index (NDSI) 29

NDSI
Red Green Red Green , www.nature.com/scientificreports www.nature.com/scientificreports/ are employed to assist the interpretation. NDSI, NGRDI and NDSI * are essentially the spectral ratios of visible bands that have similar levels of scattering either inside or outside the partially-shaded topography regions. They are not as sensitive as NDVI or NDVI * , in terms of discriminating vegetation. However, they are ideal in identifying the partially-shaded topography regions from those misinterpreted bare/non-vegetated lands.
Regarding to classifying shadows, Yang et al. 29 employed the first principal component (PC1), while Liu 30 used the intensity of panchromatic band that is proportional to PC1 32 . Ma et al. 31 took a further step to calculate the shadow index (SI): nor nor where S and I are the saturation and intensity components respectively, after converting the Red-Green-Blue color model to the Hue-Saturation-Intensity color model 33 . PC1 nor is the normalized value of PC1 according to www.nature.com/scientificreports www.nature.com/scientificreports/ To illustrate how Eqs (6) and (7) are designed to enhance the contrast between shadow and non-shadow regions, the case of study site No. 385 (annotated as the green box in Fig. 2) is given in Fig. 4 as an example. The standard false color composite (Fig. 4a) clearly shows two large shadows on the upper-left and lower-right, as well as a considerable amount of shadows scattered throughout the entire area. By utilizing three properties of shadows: lower hue (closer to black, Fig. 4b), larger saturation (more scattering light in blue band, Fig. 4c), and lower intensity (blocked of light, Fig. 4d) 33 , the map of SI is calculated using Eqs (6) and (7) and presented in Fig. 4(e). Comparing to the map of PC1 (Fig. 3f), SI indeed enhances the difference between shadow and non-shadow regions, which facilitates the determination of threshold to detect shadow.
Considering the characteristic of each spectral index: NDVI is sensitive to vegetation, SI is sensitive to shadow, and NGRDI is sensitive to bare land, we adopt the concept of decision tree: building a classification model by breaking down a dataset into smaller and smaller subsets, and map FL, SL and BL sequentially from one optical image. The rest of the region can be regarded as LVL, under the assumption that there is no other land cover/land use in the study area. The flow chart of TLDT approach is shown in Fig. 5.
The case of site No. 378 (annotated as the red box in Fig. 2) is used as an example to illustrate the procedure of processing step-by-step, and the results are shown in Figs. 6 and 7, respectively. The histogram of NDVI exhibits a typical pattern of Gaussian distribution with single peak (Fig. 6a) since most of the study area is covered by forest (Fig. 7a). The threshold T NDVI can be determined by the inflection point without ambiguity to separate the first subset: FL (Fig. 7b) from the non-vegetated area. The second subset: SL is another feature that is frequently found in an optical imagery of mountainous areas, because the sun is not usually in the nadir direction while the image is acquired 14 . Note that SI is designed to enhance the difference between the shaded and non-shaded areas. The threshold T SI can also be determined by the inflection point in the histogram of SI, as long as there is a certain ratio of shadow in the non-vegetated area. This is exactly the case of Fig. 6(b), where a considerable part of topographic shadows can be seen in the study area (Fig. 7c). Under the assumption that there is no other land cover/ land use in the study area, the rest of the region can be regarded as BL and LVL. Because the difference between the amount of BL and LVL is not as large as what we have seen for the cases of FL and SL, there is no such a typical pattern of Gaussian distribution with single peak in the histogram of NGRDI (Fig. 6c). Instead of selecting the inflection point, the Otsu's method can be employed to determine the threshold T NGRDI 31 and separate BL (Fig. 7d) from LVL (Fig. 7e). The Otsu's method assumes that the image contains two classes and calculates the optimum threshold separating the two classes so that their inter-class variance is at the maximum. All boundaries of FL (green lines), SL (white lines), BL (yellow lines) and LVL (cyan lines) are overlaid on the true color composite and shown in Fig. 7(f) that gives a reasonable result by visual examination. To assess the effectiveness of forest restoration, this TLDT approach is applied to the time series of Formosat-2 images acquired at the same site. Then, the results are compared with the ones obtained independently from the high spatial resolution (25 cm) aerial orthophotos that were acquired approximately the same time in which the Formosat-2 images (8 m) were taken.

Validation.
To validate the results of FL, SL, BL and LVL mapped from the TLDT approach, the aerial orthophoto of site No. 378 taken on 3 July 2012 is used as an alternative source of ground truth, shown in Fig. 8(a). The same region is cut from the Formosat-2 image taken on 7 July 2012, shown in Fig. 8(b) for comparison. The results of classification of FL, SL, BL and LVL are shown in Fig. 9(a,b) and the confusion matrix is given in Table 4. The same assumption is made that the region of LVL is regarded as the rest from those of FL, SL and BL. Although the overall accuracy is as high as 86.4%, the commission error and omission error for LVL are both low (42.8%  www.nature.com/scientificreports www.nature.com/scientificreports/ and 50.6%). After carefully examining Fig. 9(a,b), we confirm that the deviations are mainly found near the boundary of each class. Because the spatial resolutions of two images are intrinsically different (aerial orthophoto: 25 cm; multispectral bands of Formosat-2 image: 8 m), it is expectable that they won't give the same details near the boundary of each class. To exclude the deviations caused by difference resolution, the buffer zone with a size of one pixel is taken at the boundary of each class, as shown as the white regions in Fig. 9(c,d). The union of both buffer zones is then generated and applied to both images, in order to mask out those suspicious pixels near the boundary of each class, as shown as the white regions in Fig. 9(e,f). The confusion matrix is calculated again and the result is shown in Table 5. Not only the overall accuracy increases to 96.8%, but the commission error and omission error for LVL are both improved to 23.0% and 35.8%. This provides an alternative way to validate the TLDT approach.
Quantitative assessment of forest restoration. For a successful case of forest restoration, we expect to see a gradual decrease of BL and a gradual increase of FL. Before performing the time series analysis, the same procedure of radiance normalization as described in section 3.2 is applied to the entire time series of Formosat-2 imagery to ensure that the classification of FL, SL, BL and LVL are comparable on various dates. For site No. 378 (located in Zone 2), the root mean square errors (RMSD) of each spectral band before (cross marks) and after (circle marks) the radiance normalization are shown in Fig. 10, using the one taken on 9 June 2015 as the base image of radiance normalization. Results show that RMSD is reduced significantly by the process of radiance www.nature.com/scientificreports www.nature.com/scientificreports/ normalization. For example, the RMSD of blue band of the image taken on 10 August 2010 is reduced from the highest value of 13.6 to less than 1.0. The radiance-normalized images, are employed for classification of FL, SL, BL and LVL. Based on the classified FL, SL, BL and LVL, the area ratio A X is defined as X where X represents the class of FL, SL, BL and LVL. Considering the fact that A SL is different in different images acquired in different seasons, we can use the shadow area ratio

SL FL BL LVL
to correct the influence of shadow on class X and calculate its corrected area ratio by Tracing the variation of A FL ⁎ , ⁎ A BL and ⁎ A LVL with time provides a good indicator of restoration. Figure 11 shows the time series of A FL , A SL , A BL , A LVL , ⁎ A FL , A BL ⁎ and A LVL ⁎ at site No. 378. Though A SL indeed varies with time, it occupies about merely 5% of the study area and results in a small value of SAR (Eq. 9). The influence of shadow is not that significant, and thus A X (solid lines) are very close to A X ⁎ (dotted lines) using Eq. (10), where X represents the class of FL, BL and LVL. Both the time series of A FL ⁎ and A BL ⁎ indicate that the landslide was triggered by Typhoon Morakot in early August of 2009 when the largest A BL is found. Ever since the restoration was implemented in early 2010, ⁎ A BL (orange dotted line) has been decreased and ⁎ A FL (green dotted line) has been increased gradually. Meanwhile, ⁎ A LVL (blue dotted line) has maintained at an approximately constant level that was slightly higher in the early stage and slightly lower in the later stage of restoration. This sheds some light on the process of forest restoration: A FL ⁎ is only one-third to a half of ⁎ A LVL in the beginning of restoration; A FL ⁎ keeps increasing while ⁎ A LVL stays at the same level; and eventually, ⁎ A FL could reach double or even triple of A LVL ⁎ in the later stage of restoration. Figure 11 also reveals the seasonal variations of ⁎ A FL , ⁎ A BL and ⁎ A LVL by annotating four seasons with four difference background colors. Generally speaking, the seasonal variations of ⁎ A FL are more apparent before Typhoon Morakot than in the later stage of restoration. Because the seasonal variations of ⁎ A FL is one of the methods of representing the function of plant metabolism, ⁎ A FL can be regarded as another indicator of forest restoration. In any case, the long-term trend of successful forest restoration at site No. 378 is clear, shown by the comparison of the annual data of the same season and the regression lines in Fig. 11. www.nature.com/scientificreports www.nature.com/scientificreports/

Results and Discussion
Following the same procedure described in section 3.3, the results of FL, SL, BL and LVL mapped from the TLDT approach at the other 14 test sites are compared to the results derived from the aerial orthophotos. The overall accuracy for each test site is listed in Table 1. As explained and emphasized in section 3.3, these values do not serve as an accuracy report, because the spatial resolutions of source images are intrinsically different. Despite those deviations found near the boundary of each class, the overall consistency of classification for all 14 test sites demonstrates that the TLDT approach is valid. The same approach is employed to assess the effectiveness of forest restoration for all 30 study sites, and the time series of ⁎ A BL are grouped, shown in Fig. 12 for comparison. The results of effectiveness assessment are also listed in Table 1, which indicate that the forests are restored efficiently at ten sites (No. 61, 378, 385, 402, 411, 430, 435, 442, 443 and 452) (Fig. 12a), but restored in a slow pace at four sites No. 440, 459, 460 and 482 (Fig. 12b). The other four sites, No. 445, 450, 461 and 462 show that the landslide areas have been stabilized or even expanded and there is no significant change of forests. Therefore, the forest restoration is ineffective (Fig. 12c). The rest of the 12 sites (No. 51, 54, 55, 355, 391, 400, 415, 424, 431, 451, 454 and 463) are all located in the shaded sides of mountain areas, where a considerable fraction of Formosat-2 image are severely covered by shadows (Fig. 12d). Consequently, no assessment of the effectiveness of forest restoration can be drawn from the time series of Formosat-2 images.
The most popular and widely-used engineering method of forest restoration in Taiwan is seeds spreading. In this research, a total of 27 sites adopted this method, though a few of them also employed other methods such as staking and wattling, drainage of longitudinal and transverse, tree planting, and covering the grass net. After excluding those 10 sites with no assessment due to shadows, 9 sites (52.9%) are restored well apparently, 4 sites (23.5%) are restored in a slow pace, and 4 sites (23.5%) are ineffective. In other words, the successful rate of forest restoration in Taiwan is 76.5%. The reason of unsuccessful restoration is also clarified by denoting the dates of www.nature.com/scientificreports www.nature.com/scientificreports/ serious earthquake, typhoon, and storm events in Fig. 12(c). The secondary or even tertiary disasters destroyed the vulnerable newly-restored land. As a result, all earlier efforts of restoration dissipated and BL was not converted to FL. Such a failure has nothing to do with the engineering method of seeds spreading. Nevertheless, a timely response to the secondary disaster can mitigate the damage. This can be achieved by acquiring the aftermath remote sensing imagery, analyzing with the approach proposed in this research, and comparing with the time series of observations at the same site. Because almost all study sites adopted the similar method of seeds spreading, there is no intention to specify how or whether the restoration technique affected recovery in this work. It is better not to include many types of restoration technique at the stage of developing and evaluating a new approach. Future works have been planned to employ this new approach to investigate how different restoration techniques affect recovery.
The distinction between natural regeneration and regrowth enhanced by restoration efforts can be made by using the existing dataset and method developed in this research. As shown in Fig. 2, the landslides at NRS 1 and 2 were also triggered by Typhoon Morakot but they were not selected as restoration sites by FBT. Therefore, the changes found at NRS 1 and 2 are purely the effect of natural regeneration, which can help us to gain insight into the difference between natural regeneration and restoration. The time series of A FL ⁎ , A LVL ⁎ and ⁎ A BL at sites No. 378, NRS 1 and NRS 2 are plotted in Fig. 13 for comparison, and their corresponding linear regression equations are listed in Table 6  www.nature.com/scientificreports www.nature.com/scientificreports/ highlights the importance of restoration for it accelerates the natural regeneration to at least four times, in terms of the slopes of the linear regression equation of ⁎ A FL listed in Table 6 (2.0 × 10 −4 vs. 4.492 × 10 −5 ). Considering the scale of landslides in Taiwan's mountain areas, the high spatial resolution satellite imagery or aerial orthophotos are preferable sources of data in the management of forest watershed. For example, the preparation of a long-term national landslide inventory is based on 2 m resolution Formosat-2 imagery 30 , while the events responding to landslides and debris flow are based on tens of centimeter resolution aerial orthophotos acquired from a low-cost unmanned aerial vehicle 34 . Though these optical imagery with multispectral bands (e.g. red, green, blue and near-infrared) have high spatial resolution, they lack the critical spectral bands to retrieve the amount of water vapor and aerosol. Consequently, it is very unlikely to conduct a rigorous atmospheric correction interference before conducting the time series analysis. This work demonstrates an alternative method, confirming that the relative calibration based on the technique of radiance normalization enables us to map FL, SL, BL and LVL using the TLDT approach and meets the requirement of time series analysis. The process of radiance normalization also provides an appropriate way of examining the quality of every image.
Geographical locations of forest restoration sites in Taiwan are all located in mountainous areas, where SL is one of the main features that are inevitably found in an optical imagery 14 . In this research, 12 of 30 sites (40%) are severely covered by shadows and no assessment of the effectiveness of forest restoration can be drawn from the time series of Formosat-2 images. Some shadows are casted by topographic relief and may not always be completely dark. They might be illuminated by scattering light and can be recovered by considering the topographic relief after rigorous atmospheric correction. Thanks to the success of Sentinel-2 mission 35 , high-temporal (5 days), -spatial (10 m), -spectral (13 bands), and -radiometric (12 bits), data are now available systematically and freely to all registered users. European Space Agency also developed and released an official tool, Sen2Cor, which is able to perform a rigorous atmospheric correction and produce the Bottom-Of-Atmosphere reflectance images 36 . Together with the usage of 90 m SRTM Digital Elevation Database, Sen2Cor is able to conduct terrain correction for rugged terrains. In light of the successful launch/operation of Sentinel-2A/B and the following new identical missions continuing to take the data record to the 2030 time frame 37 , future research is being proposed to investigate the feasibility of employing Sentinel-2 multitemporal imagery to assess the restoration of forest based on the experience acquired and technique developed in this work.

Conclusion
Climate variability and man-made impacts have severely damaged the forest around the world in recent years, which calls for an urgent need of restoration aiming toward long-term sustainability for the forest environment. To assess the effectiveness of forest restoration requires a systematic and synoptic view for monitoring the forest. However, this cannot be achieved by the traditional method of in situ site surveying, particularly for those inaccessible sites with large-scale restoration. Remote sensing imagery by far is more advantageous in assessing forest restoration, mainly due to its ability to detect changes in large areas over long periods of time that are www.nature.com/scientificreports www.nature.com/scientificreports/ difficult to observe from the ground. Considering the cost, restoration scale and data availability in Taiwan, the multitemporal, multispectral and high-spatial-resolution Formosat-2 imagery acquired between 2005 and 2016 are employed to assess the effectiveness of forest restoration at 30 sites. However, Formosat-2 imagery cannot go through a rigorous atmospheric correction due to the lack of critical spectral bands to retrieve the amount of water vapour and aerosol. This paper proposes a new TLDT approach to map FL, SL, BL and LVL sequentially   Table 4. Confusion matrix of validating the results of FL, BL and LVL mapped from the TLDT approach. † The aerial orthophoto of site No. 378 was taken on 3 July 2012, as shown in Fig. 8(a). ‡ Formosat-2 image was taken on 7 July 2012, as shown in Fig. 8 Table 5. Confusion matrix of validating the results of FL, BL and LVL mapped from the TLDT approach, after excluding the deviations near the boundary of each class caused by the different resolutions between Formosat-2 imagery and Aerial orthophoto. † The aerial orthophoto of site No. 378 was taken on 3 July 2012, as shown in Fig. 8(a). ‡ Formosat-2 image was taken on 7 July 2012, as shown in Fig. 8(b).  www.nature.com/scientificreports www.nature.com/scientificreports/ from the remote sensing imagery by integrating three spectral indices: NDVI, SI and NGRDI TLDT requires neither image normalization nor atmospheric correction, and improves on the other methods by introducing more levels of decision tree classification with inputs from the same multispectral imagery. This paper also demonstrates that the relative calibration based on the technique of radiance normalization meets the requirement of time series analysis of Formosat-2 imagery. The process of radiance normalization provides an appropriate way to examine the quality of every image. With TLDT, the effectiveness of forest restoration at 30 sites are assessed, using all available multispectral Formosat-2 images acquired between 2005 and 2016. The assessments at 15 test sites are compared with the results obtained independently from the high-spatial-resolution (25 cm) aerial www.nature.com/scientificreports www.nature.com/scientificreports/ orthophotos that were acquired approximately the same time in which the Formosat-2 images (8 m) were taken. Results show that the overall consistency is between 69.14% and 87.51%, while the kappa values are between 0.46 and 0.76. These comparisons are not intended to conclude the assessment accuracy but to clarify why and where the discrepancy is, since the spatial resolutions are rather different between the aerial orthophoto (25 cm) and Formosat-2 image (8 m). Among the 30 study sites, 10 have restored successfully, 4 are recovering slowly, and 4 have hardly re-vegetated. For the rest of the 12 sites, the shaded areas are too large to derive a detailed trend of restoration. But the effectiveness can still be assured by examining the pre-restoration and the most up-to-date Formosat-2 images. The distinction between natural regeneration and regrowth enhanced by restoration efforts were also made by using the existing dataset and TLDT developed in this research. This work supports the use of multitemporal remote sensing imagery as a reliable source of data for assessing the effectiveness of forest restoration on a regular basis. The global trend of forest restoration can be obtained by applying the same approach to Landsat dataset, the longest archive of space-based moderate-resolution land remote sensing data. Therefore, this work also serves as the basis for studying the global trend of forest restoration in the future.  Table 6.   Fig. 13.