Tree Line Identiﬁcation and Dynamics under Climate Change in Wuyishan National Park Based on Landsat Images

: The alpine tree line ecotone, reﬂecting interactions between climate and ecology, is very sensitive to climate change. To identify tree line responses to climate change, including intensity and local variations in tree line advancement, the use of Landsat images with long-term data series and ﬁne spatial resolution is an option. However, it is a challenge to extract tree line data from Landsat images due to classiﬁcation issues with outliers and temporal inconsistency. More importantly, direct classiﬁcation results in sharp boundaries between forest and non-forest pixels / segments instead of representing the tree line ecotone (three ecological regions—tree species line, tree line, and timber line—are closely related to the tree line ecotone and are all signiﬁcant for ecological processes). Therefore, it is important to develop a method that is able to accurately extract the tree line from Landsat images with a high temporal consistency and to identify the appropriate ecological boundary. In this study, a new methodology was developed based on the concept of a local indicator of spatial autocorrelation (LISA) to extract the tree line automatically from Landsat images. Tree line responses to climate change from 1987 to 2018 in Wuyishan National Park, China, were evaluated, and topographic e ﬀ ects on local variations in tree line advancement were explored. The ﬁndings supported the methodology based on the LISA concept as a valuable classiﬁer for assessing the local spatial clusters of alpine meadows from images acquired in nongrowing seasons. The results showed that the automatically extracted line from Landsat images was the timber line due to the restriction in spatial autocorrelation. The results also indicate that parts of the tree line in the study area shifted upward vertically by 50 m under a 1 ◦ C temperature increase during the period from 1987 to 2018, with local variations inﬂuenced by slope, elevation, and interactions with aspect. Our study contributes a novel result regarding the response of the alpine tree line to global warming in a subtropical region. Our method for automatic tree line extraction can provide fundamental information for ecosystem managers.


Introduction
The alpine tree line ecotone provides important ecosystem services (hotspots of biodiversity [1], slope stability, high-quality water for downstream areas [2], nutrient input, and carbon sequestration for lower-elevation ecosystems [3]) and reflects interactions among climate, species ecology, physiography, and physiology [4,5]. It is also very sensitive to climate change [6][7][8] and is one of the most sensitive Among studies that have focused on the movement of tree lines in the past century, several have observed that tree lines move upward in mountainous regions [10,12,[16][17][18][19][20][21][22][23][24], while others have not observed any clear upward [25][26][27] or downward [2,28,29] movement of tree lines. Four reasons have been cited to explain these inconsistent research findings. First, the research time period is typically short. Time scale (short period, median period, and long-term) is very important for tree line advancement research [30]. As Lloyd [31] indicated, the mean time lag between initial tree recruitment and forest development is 200 years, but study duration have comprised a fraction of this. Second, the intensity of upward tree line shift varies between climate zones, largely due to the uneven global climate change effects [32], as well as significantly different dominant vegetation species and communities correlated with climate zones [33]. Third, historical anthropogenic disturbances (e.g., deforestation) have caused upward tree line shifts [16,24], where they may be carried out due to climate conditions [17]. This creates an enormous challenge in studying tree line upslope advance responses to climate change [2]. Finally, temperature has not been a primary factor dominating tree line upward advancement under long-term global warming until now [28]. Even within a single area, strong local variations in tree line shifts have been observed in many studies [10,12,20,34,35]. Therefore, it is important to improve our understanding of tree line shifts with climate change and to explore other environmental factors that can cause local variation in tree line advancement. The long-term historical mapping of tree lines with detailed local spatial information, as opposed to studies in several elevational transitions [12,23,36,37], is fundamental for modern studies of tree line shifting.  [15].
Among studies that have focused on the movement of tree lines in the past century, several have observed that tree lines move upward in mountainous regions [10,12,[16][17][18][19][20][21][22][23][24], while others have not observed any clear upward [25][26][27] or downward [2,28,29] movement of tree lines. Four reasons have been cited to explain these inconsistent research findings. First, the research time period is typically short. Time scale (short period, median period, and long-term) is very important for tree line advancement research [30]. As Lloyd [31] indicated, the mean time lag between initial tree recruitment and forest development is 200 years, but study duration have comprised a fraction of this. Second, the intensity of upward tree line shift varies between climate zones, largely due to the uneven global climate change effects [32], as well as significantly different dominant vegetation species and communities correlated with climate zones [33]. Third, historical anthropogenic disturbances (e.g., deforestation) have caused upward tree line shifts [16,24], where they may be carried out due to climate conditions [17]. This creates an enormous challenge in studying tree line upslope advance responses to climate change [2]. Finally, temperature has not been a primary factor dominating tree line upward advancement under long-term global warming until now [28]. Even within a single area, strong local variations in tree line shifts have been observed in many studies [10,12,20,34,35]. Therefore, it is important to improve our understanding of tree line shifts with climate change and to explore other environmental factors that can cause local variation in tree line advancement. The long-term historical mapping of tree lines with detailed local spatial information, as opposed to studies in several elevational transitions [12,23,36,37], is fundamental for modern studies of tree line shifting.
Remote-sensing approaches have the potential to capture spatial variations in tree line positions and enable long-term analyses of upward tree line advancement. Previous studies applying remote-sensing methods to tree line ecotones have followed one of three directions: evaluating biophysical parameters in tree line ecotones, measuring topographic factors that impact tree line dynamics, or measuring tree line position and shifting. Biophysical parameter research has used temporal canopy cover information [13,21,59], tree density [60], species abundance [61][62][63], fraction of absorbed photosynthetically active radiation (FAPAR), and vegetation water content [64] measured from remote-sensing data (or directly from temporal changes of vegetation indices [65,66]) as indicators for tree line dynamics. Second, digital elevation models (DEMs) have been applied to extract parameters to analyze topographic effects on tree line position and tree density along the tree line [13,67]. Finally, the classification of forest and non-forest areas from moderate-resolution imagery [13,16,17,24,68] and the detection of individual trees and tree islands from high-resolution imagery [21,[69][70][71][72], aerial photos, or Lidar data [73,74] have been reported in the literature with regard to measuring the spatial variation of tree line position and the temporal change of tree line advancement. Nevertheless, few studies have attempted to use a searching window to find the highest-elevation forest pixels in order to identify potential tree line locations, as determined by climate factors instead of anthropogenic disturbance [17].
However, it remains a challenge to extract temporal changes in tree line position using remotely sensed images due to mismatches between ecological conditions and remote-sensing outputs. One issue comes from spatial resolution and the length of time over which tree line dynamics are studied. Long-term remote-sensing imagery (e.g., Landsat with a 30 m spatial resolution) has a limited ability to analyze the upward advancement of a tree line due to climate change if the shift is less than a single pixel's breadth (e.g., the tree line of the Khibiny Mountains in Russia advanced about 27 m vertically over 55 years [21]). However, when the slope is less than 45 • , the horizontal distance is larger than vertical distance, which gives Landsat imagery the potential to monitor tree line advancement within its 30 m spatial resolution.
Another mismatch is that the direct output of remote-sensing classification from forest and non-forest pixels/segments is usually sharp boundaries between alpine meadows and forests instead of the location of a tree line ecotone, even when the spatial resolution of the image is high enough to detect tree islands. Therefore, it is unclear which ecotone type (i.e., tree species line, tree line, or timber line; Figure 1) is extracted from remote-sensing images. Moreover, results from conventional pixel-based classification approaches (e.g., supervised or decision-tree classification) normally include erroneous land cover, which is a challenge for long-term temporal change analysis [75]. Maintaining a consistent remote-sensing dataset of spectral reflectance is a further challenge that must be overcome to ensure long-term the consistency of remote-sensing applications in ecosystems, especially with a background of rapid global change [76]. Thus, remote sensing has its own particular challenges in measuring tree line advancement with climate change due to the heterogeneity of ecotones caused by individual trees and tree islands ( Figure 1). Therefore, it often requires expert knowledge to conduct post-classification Remote Sens. 2020, 12, 2890 4 of 25 analysis and change detection to identify land-cover errors in order to properly attribute land-cover change [75].
There is also inconsistency in detecting temporal change from remote-sensing images. Studies have used pixel clustering on homogeneous segments with thresholding algorithms prior to classification to reduce erroneous results [77]. Local indicator of spatial association (LISA), another method that captures spatial clusters, performs well as a classifier for ecological research (e.g., burn severity [78] and vegetation fragmentation in urban areas [79]). However, LISA, as a geostatistical technique, is still rarely used as a classifier of remotely-sensed images even though it has been shown to be more accurate than conventional classification approaches [77]. These studies have also shown that introducing spatial weights (a major parameter of LISA) would improve the temporal consistency of land-cover classification from remote-sensing images [80]. For tree line mapping, LISA may have great potential to distinguish alpine meadows from closed forests and tree line ecotones because high spatial autocorrelation is common within alpine meadows, in contrast to heterogeneous tree line ecotones and closed forests.
This study aimed to develop a new detection method based on the LISA concept in order to detect tree line position and dynamics from Landsat imagery. The specific objectives were (1) to use LISA to develop a new method of automatically identifying tree line location from Landsat imagery and explaining detected ecotone types, (2) to evaluate the dynamics of tree line positioning over 30 years, (3) to analyze the response of upward tree line shifts to climate change, and (4) to explore topographic effects on local variations of tree line advancement.

Study Area
The study area was in the boundary of Wuyishan National Park in the northern Fujian and Jiangxi Provinces, China (27 • 49 -27 • 53 N, 117 • 45 -117 • 49 E; Figure 2). This area, characterized as a subtropical evergreen broad-leaf forest, has been identified as important for its contribution to world biodiversity [81]. Wuyi Mountain (Wuyishan) is included in the World Natural and Cultural Heritage list (23rd Session of the World Heritage Committee; December 1999) [82]. Wuyishan was first declared a national park in 2017, and its area includes Wuyishan National Nature Reserve (NNR; 565.27 km 2 ), the Nine-Bend Stream Ecological Protection Area (NEPA; 353.32 km 2 ), and the Wuyishan National Scenic Area (NSA; 64 km 2 ) [83].
The region has a typical humid subtropical monsoon climate that is characterized by warm summer temperatures, cool winter temperatures, and abundant annual precipitation [84]. Across the elevation gradient in the study area, the mean annual rainfall is approximately 2050, 2530, and 2820 mm at 1200, 1600, and 2000 m a.s.l., respectively, with approximately 70% of the precipitation falling during the growing season (from March to August) [85]. Nearly all the precipitation is rainfall instead of snowfall, however, extreme weather events (e.g., El Nino effects) can cause snowfall in winter time (e.g., strong snowfall in 2008). The annual mean temperature is 13, 11, and 9 • C at 1200, 1600, and 2000 m a.s.l., respectively [86]. The month with the lowest temperature is January or February, and the average temperature of the coldest month is 7.6 • C. There is no prevailing wind direction at the top of the mountain due to valley wind from different directions, but westerly and northwesterly winds are more frequent in the winter. The soil type is classified as yellow-red soil based on the Chinese soil classification system, which is equivalent to Ultisol in the US Department of Agriculture (USDA) soil taxonomy [87]. Along the altitude gradient from high to low, there are five distinct vegetation communities (Figure 3): the alpine meadow (AM; Figure 2b), the subalpine dwarf forest (DF; Figure 2b), the subalpine coniferous forest (CF; Figure 2b), the mixed forest, and the evergreen broad-leaf forest [88,89]. Along the tree line, the alpine meadow is dominated by Miscanthus sinensis, Deyeuxia arundinacea, and Molinia japonica; the subalpine dwarf forest is dominated by Buxus sinica (Rehd. et Wils.) and Cheng ex M. Cheng subsp. sinica var. parvifolia M. Cheng, and the coniferous forest is dominated by Pinus taiwanensis Hayata, Tsuga chinensis (Franch.), Pritz. var. tchekiangensis (Flous) Cheng et L.K. Fu, and Pinus massoniana Lamb. Overall, the climate conditions characterized by low temperature, low pressure, and strong wind formed the alpine meadow in the study area. In some areas, the tree line is dominated by a transactional ecotone between the dwarf forest and the alpine meadow, while in other areas, the tree line is dominated by a transactional ecotone between the coniferous forest and the alpine meadow (Figures 2b and 3).
topographic effects on local variations of tree line advancement.

Study Area
The study area was in the boundary of Wuyishan National Park in the northern Fujian and Jiangxi Provinces, China (27°49′-27°53′N, 117°45′-117°49′E; Figure 2). This area, characterized as a subtropical evergreen broad-leaf forest, has been identified as important for its contribution to world biodiversity [81]. Wuyi Mountain (Wuyishan) is included in the World Natural and Cultural Heritage list (23rd Session of the World Heritage Committee; December 1999) [82]. Wuyishan was first declared a national park in 2017, and its area includes Wuyishan National Nature Reserve (NNR; 565.27 km 2 ), the Nine-Bend Stream Ecological Protection Area (NEPA; 353.32 km 2 ), and the Wuyishan National Scenic Area (NSA; 64 km 2 ) [83]. The region has a typical humid subtropical monsoon climate that is characterized by warm summer temperatures, cool winter temperatures, and abundant annual precipitation [84]. Across the elevation gradient in the study area, the mean annual rainfall is approximately 2050, 2530, and 2820 Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 26 mm at 1200, 1600, and 2000 m a.s.l., respectively, with approximately 70% of the precipitation falling during the growing season (from March to August) [85]. Nearly all the precipitation is rainfall instead of snowfall, however, extreme weather events (e.g., El Nino effects) can cause snowfall in winter time (e.g., strong snowfall in 2008). The annual mean temperature is 13, 11, and 9 °C at 1200, 1600, and 2000 m a.s.l., respectively [86]. The month with the lowest temperature is January or February, and the average temperature of the coldest month is 7.6 °C. There is no prevailing wind direction at the top of the mountain due to valley wind from different directions, but westerly and northwesterly winds are more frequent in the winter. The soil type is classified as yellow-red soil based on the Chinese soil classification system, which is equivalent to Ultisol in the US Department of Agriculture (USDA) soil taxonomy [87]. Along the altitude gradient from high to low, there are five distinct vegetation communities ( Figure 3): the alpine meadow (AM; Figure 2b), the subalpine dwarf forest (DF; Figure 2b), the subalpine coniferous forest (CF; Figure 2b), the mixed forest, and the evergreen broad-leaf forest [88,89]. Overall, the climate conditions characterized by low temperature, low pressure, and strong wind formed the alpine meadow in the study area. In some areas, the tree line is dominated by a transactional ecotone between the dwarf forest and the alpine meadow, while in other areas, the tree line is dominated by a transactional ecotone between the coniferous forest and the alpine meadow (Figures 2b and 3).

Remotely Sensed Imagery
A total of 54 Landsat Thematic Mapper (TM) images taken from 1987 to 2011, 33 Enhanced Thematic Mapper Plus (ETM+) images taken from 2000 to 2012, and 14 Operational Land Imager (OLI) images taken from 2013 to 2018 were used in this study to extract the tree line in Wuyishan National Park (Table 1). All images had been previously geometrically and atmospherically corrected (Landsat Level-2 from Landsat Collection 1: surface reflectance data, https://www.usgs.gov/land-resources/nli/ landsat/landsat-collection-1) when downloaded from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/). Based on the Landsat Collection 1 information, the algorithm used for generating surface-reflectance data of Landsat TM and ETM+ was the Landsat Ecosystem Distribution Adaptive Processing System (LEDAPS), and that of Landsat OLI was the Land Surface Reflectance Code (LaSRC). The spatial resolution of all images was 30 × 30 m, and the projection was UTM (Universal Transverse Mercator) Zone 50N, WGS1984 (World Geodetic System 1984). The digital numbers of the Landsat Level 2 imagery were the reflectance values, with a scale factor of 10,000 (e.g., if the digital number is 4230, the reflectance is 0.423).
In total, 31 historical RGB (Red, Green and Blue: true color composition) images were downloaded from Google Earth (the software used for downloading was BIGEMAP, Chengdu, China). Among these, 26 images with a 16.9 × 16.9 m spatial resolution were acquired from 1987-2016, within which the images from 1991, 1994, 1997, and 2013 were left out due to poor quality. All images were acquired for the same date ( DEM data were from the ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) Global Digital Elevation Model V002 (GDEM2). One DEM image of the study area was downloaded from the NASA website (https://reverb.echo.nasa.gov/reverb/). The spatial resolution for the DEM image was 30 × 30 m, and it was projected as UTM Zone 50N, WGS1984. The horizontal errors for GDEM2 data were −0.13 arc-seconds (east/west shift) and −0.19 arc-seconds (north/south shift). The vertical error was −0.74 m on average, based on Japanese studies, ranging from −5.58 to +15.45 m-with 8.68 m over forest regions and −1.09 m over farmland. The horizontal and vertical accuracy were drawn from an ASTER GDEM2 validation report (https://ssl.jspacesystems.or.jp/ersdac/ GDEM/ver2Validation/Summary_GDEM2_validation_report_final.pdf).

Methods
To develop a new method for extracting tree lines and identifying the type of extracted lines (i.e., tree line, timber line, or tree species line), potential available bands for ecotone detection were tested for automatic tree line identification from Landsat series images ( Figure 4). The optimal season for tree line extraction (i.e., acquisition time for Landsat images) was also evaluated. To achieve our objectives, climate impacts on upward tree line shifts, as well as local variations in the intensity of upward tree line movement, were measured. This research also explored the topographic effects on local variations in the intensity of upward tree line shifts.

Automatic Tree Line Extraction from Remote-Sensing Imagery
Preprocessing for all Landsat images, including TM, ETM+, and OLI images, involved clipping to the specific study area and conversion to point format from rasters. LISA (in this study, Anselin Local Moran's I in the ArcGIS Spatial Statistics toolbox was used as the LISA tool) was then applied to the point data for each Landsat image. LISA was used to capture hot spots (i.e., spatial clusters with high values clustering together) and cold spots (i.e., spatial clusters with low values clustering together) in images based on the concept of spatial autocorrelation between pixels and their

Automatic Tree Line Extraction from Remote-Sensing Imagery
Preprocessing for all Landsat images, including TM, ETM+, and OLI images, involved clipping to the specific study area and conversion to point format from rasters. LISA (in this study, Anselin Local Moran's I in the ArcGIS Spatial Statistics toolbox was used as the LISA tool) was then applied to the point data for each Landsat image. LISA was used to capture hot spots (i.e., spatial clusters with high values clustering together) and cold spots (i.e., spatial clusters with low values clustering together) in images based on the concept of spatial autocorrelation between pixels and their surrounding pixels. Taking the normalized difference vegetation index (NDVI) as an example, the alpine meadow had lower NDVI values than other vegetation communities (coniferous and dwarf forests). Therefore, the spatial arrangement for the alpine meadow in images was spatial clusters with low NDVI values clustered together, especially in the grass senescence season when the subalpine forest was evergreen. Thus, low-low clusters (spatial clusters with low values) were selected for the NDVI and near-infrared band (NIR) (the reflectance of NIR was lower for the alpine meadow than for the subalpine forest) of point data with a p-value equal to 0.002 (statistically significant). High-high clusters (spatial clusters with high values) were selected for the blue band (B), green band (G), red band (R), first shortwave infrared band (SWIR1), and second shortwave infrared band (SWIR2) of point data (p-value = 0.002) because the alpine meadow had a higher reflectance in those bands than the subalpine forest (coniferous and dwarf forests in the study area). Statistically, a p-value < 0.05 indicated significance for low-low or high-high clusters. We selected the lowest p-value (0.002 was the lowest p-value result from Anselin Local Moran's I) to capture the smoother alpine meadow in order to decrease the heterogeneity effects of isolated tree island regions in the tree line ecotone. Postprocessing steps included analyzing the point density (30 × 30 m spatial resolution, consistent with Landsat images) and selecting point densities larger than 0.0005. This processing step converted the selected spatial cluster points (low-low or high-high clusters) into the raster format. The 0.0005 threshold was set to eliminate some small patches of the alpine meadow in the study area; however, it could be adjusted. If it was set lower, more small patches of the alpine meadow would be extracted. Steps also included converting rasters to polygons, eliminating small polygons (area < 2000 m 2 ; this threshold setting, about the area of two pixels, was used to reduce the "salt-and-pepper" effect in images), smoothing polygons, and converting polygons to the tree line (line feature). The shadow effect was an issue when extracting tree lines from NDVI images, so NIR reflectance was used to control it (controlling NIR reflectance values higher than 0.035 to eliminate shadow effects; after exploring the spectral characteristics between mountain shadows and the alpine meadow in the study area, this threshold sufficiently reduced the noise effect of shadows). All steps for extracting the tree line from bands were coded as a Python script, which is available for download (see Supplementary Materials).
To test the impact of the spatial resolution of imagery on tree line extraction, RGB images acquired on 23 December 2014 with spatial resolutions of 1.06 × 1.06, 2.11 × 2.11, 4.22 × 4.22, and 8.42 × 8.42 m were tested and compared.

Selecting Optimum Band and Season for Tree Line Extraction
The strategy for selecting the optimum band and season for tree line extraction was to compare all extraction results with the reference RGB image acquired in the same year. The boundaries between different vegetation communities do not change much within a year, so the variations of tree line extraction from different bands acquired for the same year were caused by extraction accuracy. It is necessary to choose the optimum extraction band and season for consistent temporal analysis. Vegetation phenology causes a high variation of spectral signals for the same vegetation type (e.g., this was noteworthy for the alpine meadow-the spectral signals of evergreen vegetation communities do not change much within a year), which may result in different extraction results even with images acquired in the same year. In addition, an accuracy assessment with the reference RGB image (spatial resolution of 0.53 × 0.53 m, with individual trees visible) acquired on 23 December 2014 was conducted for tree line extraction from the Landsat OLI images acquired on 2 January, 17 October, and 20 December 2014 to select the optimum extraction season.

Identification of Type of Extracted Tree Line
The tree line ecotone in the study area is transitional (the mosaic region of the alpine meadow and tree islands or individual trees; see Figure 1) between the coniferous forest and the alpine meadow or between the dwarf forest and the alpine meadow. In other words, the tree line ecotone potentially has an upper boundary (tree species line; Figure 1), lower boundary (timber line; Figure 1), and tree line (boundary between tree islands and individual trees; Figure 1). However, it is interpreted as a boundary line between the alpine meadow and forest when using remote-sensing imagery. Therefore, it is important to identify which line (timber line, tree line, or tree species line) is being extracted based on the LISA concept used in this study. The extracted line was identified visually by comparison to RGB images in which individual trees and tree islands were clearly visible at a high spatial resolution (0.53 × 0.53 m).

Time-Series Analysis and Segmented Linear Regression for Climate Data
To accurately capture the annual trend of temperature change, the seasonal effects of monthly mean temperature were reduced by a time-series analysis (TSA) using the Extensible Time Series (xts) package in the R software (https://www.r-project.org/). From this, the annual temperature trend data were analyzed using the segmented linear relationship with the Segmented package in the R software. The segmented linear relationship, which is also called the broken-line relationship, is widely used in ecological research [90]. In this study, the segmented linear relationship was used to analyze thresholds in temperature change trends and to evaluate the alpine meadow area trend.

Topographic Effects on Local Variations of Upward Tree Line Shift
The topographic factors of aspect, slope, and contour were estimated from DEM data using surface analysis in the ArcGIS software. To analyze the topographic effects on local variations of tree line advancement along an altitudinal gradient with climate change, a t-test was used to compare whether there was a significant difference of slope and elevation between stable and rapidly upward shifting tree lines. A two-way ANOVA was conducted to determine whether there were significant differences of slope and elevation between different tree line types (stable and advanced tree line) and aspects (flat, north, northeast, east, southeast, south, southwest, west, and northwest). The samples used for testing topographic effects were from the pixel values (slope, elevation, or aspect) where the tree line was located for the stable tree line, as well as from all pixels in the location between the current advanced tree line and previous tree lines during the study period.

Automatic Extraction of Tree Line from Landsat Imagery
With our methodology for automatic tree line extraction, the NDVI performed best of all bands based on the visual comparison with the RGB reference images (16.9 × 16.9 m spatial resolution) for each year ( Figure 5). According to the spectral signatures of different vegetation communities around the tree line ecotone (alpine meadow, subalpine coniferous forest, and subalpine dwarf forest in the study area), the use of the NIR band to detect the alpine meadow was the least accurate. Therefore, significant low-low clusters (low values significantly clustered) from spatial autocorrelation analysis were selected for the alpine meadows. However, the results from the NIR bands were greatly influenced by mountain shadows ( Figure 5) because they had even lower reflectance in the NIR bands than the alpine meadow. Based on these constraints, the reflectance values of B, G, R, SWIR1, and SWIR2 in the alpine meadow were higher than those of other vegetation communities because dead materials cause higher reflectance in the B, G, and R bands than green vegetation and they hold less water content, resulting in higher reflectance values in the SWIR1 and SWIR2 bands. Compared to tree line extraction from the NDVI, extraction results from B, G, R, SWIR1, and SWIR2 showed patchiness, especially the B band, due to strong atmospheric effects. Regardless, the NDVI performed optimally in extracting the tree line and was only slightly affected by mountain shadows. Therefore, the shadow effects of extraction from the NDVI band were controlled by reflectance from the NIR band because its reflectance of shadow was much lower than in alpine meadows. Reflectance values larger than 0.035 (the default value set in Python script) in the NIR band were selected to reduce the shadow effects as an additional postprocessing step.
According to tree line extraction results from available cloud-free Landsat images, the optimal extraction time was from late October to early May ( Figure 6; see tree line extraction from the NDVI of all images from Landsat series in Supplementary Material A). From the images acquired during June to September, only the disturbed region (the area surrounding the tourist trail to the mountain top) was extracted ( Figure 6). This optimal extraction season is consistent with the senescent phenology period of the alpine meadow. The alpine meadow in the study area starts to green up in early May and senesce in late September. Therefore, the extracted meadow area from 13 May was significantly smaller than that from 3 May (Figure 6), and it gradually increased from 25 September to 21 October ( Figure 6). Based on the accuracy assessment results with the reference RGB image acquired on 23 December 2014 (0.53 × 0.53 m spatial resolution), the tree line extracted from the Landsat OLI images acquired on 2 January 2014 (overall accuracy = 98%; kappa coefficient = 95.77%) and 20 December 2014 (overall accuracy = 98%; kappa coefficient = 95.80%) was better than that acquired on 17 October 2014 (overall accuracy = 88.67%; kappa coefficient = 75.30%-see distribution of points and error matrix for accuracy assessment in Supplementary Material B). Therefore, images acquired in the completely senescent season of the alpine meadow were the best choice for tree line extraction for each year. Based on this criterion and the comparison with reference RGB images, one tree line extraction was selected for each year to study the temporal dynamics of tree line advancement (for more information on the best extractions of each year, see Supplementary Material C and D). According to tree line extraction results from available cloud-free Landsat images, the optimal extraction time was from late October to early May ( Figure 6; see tree line extraction from the NDVI of all images from Landsat series in Supplementary Material A). From the images acquired during June to September, only the disturbed region (the area surrounding the tourist trail to the mountain top) was extracted ( Figure 6). This optimal extraction season is consistent with the senescent phenology period of the alpine meadow. The alpine meadow in the study area starts to green up in early May and senesce in late September. Therefore, the extracted meadow area from 13 May was significantly smaller than that from 3 May (Figure 6), and it gradually increased from 25 September to 21 October (Figure 6). Based on the accuracy assessment results with the reference RGB image acquired on 23 December 2014 (0.53 × 0.53 m spatial resolution), the tree line extracted from the Landsat OLI images acquired on 2 January 2014 (overall accuracy = 98%; kappa coefficient = 95.77%) and 20 December 2014 (overall accuracy = 98%; kappa coefficient = 95.80%) was better than that acquired on 17 October 2014 (overall accuracy = 88.67%; kappa coefficient = 75.30%-see distribution of points and error matrix for accuracy assessment in Supplementary Material B). Therefore, images acquired in the completely senescent season of the alpine meadow were the best choice for tree line extraction for each year. Based on this criterion and the comparison with reference RGB images, one tree line extraction was selected for each year to study the temporal dynamics of tree line advancement (for more information on the best extractions of each year, see Supplementary Material C and D).  (Figure 7). However, the tree line was defined here as the elevational boundary between individual tree regions and isolated tree island regions of the alpine tree line ecotone. Therefore, tree line extraction depends on the spatial resolution of remotely sensed images. Based on the results of this study, images with a spatial resolution higher than 2.11 × 2.11 m have high potential for use to extract alpine tree lines. Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 26 From the extraction results of the RGB images with four spatial resolutions based on the methods introduced in this study, images with a 1.06 × 1.06 m spatial resolution contained the largest number of patches, including both tree islands within the alpine meadow and small meadow patches (  (Figure 7). However, the tree line was defined here as the elevational boundary between individual tree regions and isolated tree island regions of the alpine tree line ecotone. Therefore, tree line extraction depends on the spatial resolution of remotely sensed images. Based on the results of this study, images with a spatial resolution higher than 2.11 × 2.11 m have high potential for use to extract alpine tree lines.

Type of Line in Tree Line Ecotone Automatically Extracted from Landsat Imagery
After the cross-identification of the tree lines extracted using RGB images with a high spatial resolution (0.53 × 0.53 m) as reference at various spatial locations (Figure 8), the tree line identified using the methods based on the LISA concept from Landsat imagery was determined to actually be the timber line, which is the lower boundary of the tree line ecotone (Figure 1) and the upper

Type of Line in Tree Line Ecotone Automatically Extracted from Landsat Imagery
After the cross-identification of the tree lines extracted using RGB images with a high spatial resolution (0.53 × 0.53 m) as reference at various spatial locations (Figure 8), the tree line identified using the methods based on the LISA concept from Landsat imagery was determined to actually be the timber line, which is the lower boundary of the tree line ecotone (Figure 1) and the upper boundary of the closed forest (the subalpine dwarf and coniferous forests in this study area). In the largest mountaintop meadow region (Figure 8), four tree islands within the tree line ecotone with areas of 17,843.35, 4292.55, 3368.89, and 6158.43 m 2 were extracted. Many small tree islands in the tree line ecotone were not identified from Landsat imagery using the automatic extraction method developed in this study, and some small patches of the alpine meadow were also not successfully extracted ( Figure 8). However, the tree lines extracted from the images with high spatial resolutions (1.06 × 1.06 and 2.11 × 2.11 m; Figure 7) were the boundaries separating tree islands and the alpine meadow or tree islands within the tree line ecotone-not the timber lines extracted from Landsat images.

Climate Impact on Tree Line Upward Shift
The meteorological station records indicated periods of temperature stability in the study area  (Figure 9). In 2008, heavy snowfall destroyed a portion of the coniferous forest, resulting in an increased alpine meadow area.
Overall, the area of the alpine meadow decreased significantly from 1987 to 2012 (Figure 9c), which indicated that the tree line shifted upward along the elevation in the study area. However, the

Climate Impact on Tree Line Upward Shift
The meteorological station records indicated periods of temperature stability in the study area (1957-1984; Figure 9a (Figure 9). In 2008, heavy snowfall destroyed a portion of the coniferous forest, resulting in an increased alpine meadow area.
Overall, the area of the alpine meadow decreased significantly from 1987 to 2012 (Figure 9c), which indicated that the tree line shifted upward along the elevation in the study area. However, the intensity of the tree line advancement had high local variation ( Figure 10). In some locations (Figure 10a-c), the tree line was relatively stable from 1987 to 2018, but it advanced rapidly during 1987-2018 (Figure 10d-i). The tree line shifted approximately 210 m horizontally and 10 m vertically at the location noted in Figure 10d; approximately 220 m horizontally and 100 m vertically at the location noted in Figure 10e; approximately 130 m horizontally and 50 m vertically at the location noted in Figure 10f; approximately 160 m horizontally and 50 m vertically at the location noted in Figure 10g; approximately 100 m horizontally and 30 m vertically northeast, 90 m horizontally and 50 m vertically southeast, and 260 m horizontally and 65 m vertically west at the location noted in Figure 10h; and approximately 50 m horizontally and 3 m vertically at the location noted in Figure 10i. The average diameter at breast height (DBH), crown diameter, tree height, and canopy density in the coniferous forest were found to be 20 cm, 0.8 m, 15 m, and 70%, respectively. In the subalpine dwarf forest, trees are normally clustered together, and the average DBH, crown diameter, tree height, and canopy density for single trees were found to be 8 cm, 1 m, 5 m, and 85%, respectively. In the coniferous tree line ecotone of the study area, tree height and crown size were seen to decrease dramatically with increasing elevation, and shrubs are very common in the tree line ecotones.
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 26 Figure 9. (a) Temporal trend of monthly mean temperature in the study area from 1957 to 2017 (seasonal variation was removed from temperature data using time-series analysis in the R software), (b) temporal trend of monthly mean temperature in the study area from 1984 to 2017 (seasonal variation was removed from temperature data using time-series analysis in the R software), and (c) temporal trend of the area extracted for the alpine meadow using the LISA concept. Figure 9. (a) Temporal trend of monthly mean temperature in the study area from 1957 to 2017 (seasonal variation was removed from temperature data using time-series analysis in the R software), (b) temporal trend of monthly mean temperature in the study area from 1984 to 2017 (seasonal variation was removed from temperature data using time-series analysis in the R software), and (c) temporal trend of the area extracted for the alpine meadow using the LISA concept.

Topographic Effects on Spatial Variations of Tree Line Response to Climate Change
Slope, elevation, and aspect were significantly different between the advanced and stable tree lines (t-test, p < 0.001; see descriptive statistics and normality test in Supplementary Material E). The tree lines in locations with steeper slopes were more stable than that in areas with shallower slopes (Figure 11a). The tree line was more stable at higher elevations than at lower elevations (Figure 11c). The tree line advancement on slopes and elevations interacted with aspect (Figure 11b,d) (the results of two-way ANOVA showed significant differences for both slope and elevation; for more information on these results and Tukey HSD (Honestly Significant difference) test results, see Supplementary Material E). In all directions, tree lines with higher slopes were more likely to be stable; however, the impact of slope on the northern, northeastern, northwestern, and southeastern aspects was statistically significant (Figure 11b and Tables E9, E10, E12, and E14 in Supplementary Material E). Tree lines at higher elevations were generally more stable than those at lower elevations, while elevation was not a factor affecting tree line advancement on the eastern and southeastern aspects (Figure 11d and Tables E11 and 12 in Supplementary Material E). Figure 9. (a) Temporal trend of monthly mean temperature in the study area from 1957 to 2017 (seasonal variation was removed from temperature data using time-series analysis in the R software), (b) temporal trend of monthly mean temperature in the study area from 1984 to 2017 (seasonal variation was removed from temperature data using time-series analysis in the R software), and (c) temporal trend of the area extracted for the alpine meadow using the LISA concept.

Topographic Effects on Spatial Variations of Tree Line Response to Climate Change
Slope, elevation, and aspect were significantly different between the advanced and stable tree lines (t-test, p < 0.001; see descriptive statistics and normality test in Supplementary Material E). The tree lines in locations with steeper slopes were more stable than that in areas with shallower slopes (Figure 11a). The tree line was more stable at higher elevations than at lower elevations ( Figure 11c). The tree line advancement on slopes and elevations interacted with aspect ( Figure 11b,d) (the results of two-way ANOVA showed significant differences for both slope and elevation; for more information on these results and Tukey HSD (Honestly Significant difference) test results, see Supplementary Material E). In all directions, tree lines with higher slopes were more likely to be stable; however, the impact of slope on the northern, northeastern, northwestern, and southeastern aspects was statistically significant (Figure 11b and Tables E9, E10, E12, and E14 in Supplementary material E). Tree lines at higher elevations were generally more stable than those at lower elevations, while elevation was not a factor affecting tree line advancement on the eastern and southeastern aspects (Figure 11d and Tables E11 and 12 in Supplementary Material E).

Advantages and Disadvantages of Automatic Tree Line Extraction from Landsat Imagery
Though the remotely sensed datasets used in this study for the temporal changes of tree line advancement was from the same satellite series, the images (TM, ETM+, and OLI) were not exactly consistent in spectral characteristics during the study's time period. This was due to the slight

Advantages and Disadvantages of Automatic Tree Line Extraction from Landsat Imagery
Though the remotely sensed datasets used in this study for the temporal changes of tree line advancement was from the same satellite series, the images (TM, ETM+, and OLI) were not exactly consistent in spectral characteristics during the study's time period. This was due to the slight difference in sensor design for spectral band, sun, and viewing geometry. Our method was a qualitative analysis, which may have mitigated some of the impact of temporal spectral signal inconsistency [91] among the Landsat data series. Choosing the NDVI for the LISA analysis of each image to extract the tree line was an additional way to maintain temporal consistency in automatic tree line extraction. Previous research has demonstrated that NDVI variation is primarily controlled by land-cover type and vegetation phenology and secondarily affected by atmospheric conditions, thus accounting for inconsistencies between sensors [92]. The results of this study indicated that the NDVI performed best in extracting the timber line from among all bands in Landsat imagery.
One reason for this appears to be that the NDVI reduced the effect of poor Landsat image quality (e.g., the salt-and-pepper effect in Landsat TM and ETM+ images) to smooth the image for the LISA analysis. Another reason is that the NDVI minimized the difference between the individual tree region of the tree line ecotone and the alpine meadow, smoothing the two regions and helping to identify the timber line ( Figure 5). The potential of vegetation indices (other than the NDVI) for smoothing images and minimizing the difference between the alpine meadow and individual tree regions is worth exploring in future research. The extraction from B had the greatest patchiness, because B in Landsat images is the most influenced by the atmosphere (Figure 5). Tree lines can be extracted from G, R, SWIR1, and SWIR2 if the NDVI band is not available. However, lines extracted from these bands inaccurately identify the timber line, as they are affected by heterogeneity in the tree line ecotone (i.e., tree island size; Figures 5 and 7d).
The best extraction season in which to delineate the tree lines of alpine meadows is the vegetation senescence period. This is due to the NDVI being quite similar to that of other vegetation communities (Figure 3) during the growing season ( Figure 6). Overall, the methodology based on the LISA concept has the potential to distinguish tree islands, closed forests, and alpine meadows with high accuracy and temporal consistency from remote-sensing imagery with an appropriate spatial resolution. The weakness of this methodology is that it only extracts relatively homogeneous vegetation types with higher (high-high clusters) or lower (low-low clusters) reflectance or vegetation index values compared to other vegetation types. In this study, the alpine meadow was easy to extract because it was more homogeneous with a lower NDVI than coniferous, dwarf, or evergreen broadleaf forests, especially in the winter. However, the extraction results also depended on the spatial resolution. With an appropriate spatial resolution, when the alpine meadow had similar characteristics in the image bands to individual tree regions of the tree line ecotone (Figure 7), the alpine meadow and individual tree regions were distinguished from isolated tree islands of the alpine tree line ecotone and the closed forest, which led to the actual tree line location according to ecological concepts (Figure 1). With the spatial resolution of Landsat images, the alpine meadow had similar characteristics in the image bands to the tree line ecotone, including individual tree regions and tree islands. Thus, the timber line instead of the tree line was extracted from Landsat images based on the developed methodology.
The methods used here also had error sources that reduced the accuracy of tree line extraction from Landsat images. Spatial resolution, image quality, distance for spatial autocorrelation, smoothing polygons, and the elimination of small polygons were all error sources. Clouds and cloud shadows in images had a large impact on tree line extraction using this methodology (Supplementary Figure S1). The LISA methodology used here might be an alternative for extracting clouds and cloud shadows. Extraction results from B, G, R, SWIR1, and SWIR from cloudy Landsat images were part cloud with extensive patchiness (Supplementary Figure S1). This also indicated that the NDVI smoothed the image to some extent, which was helpful for tree line extraction ( Figure 5).
Postprocessing steps, including smoothing polygons and eliminating small polygons, introduced some errors in tree line extraction. Smoothing polygons might have caused some variation (about 30 m horizontal shift) in the final extracted tree line, and eliminating small polygons (a default set for the area of small polygons in Python script was 2000 m 2 ) might have lost some tree islands in the tree line ecotone. The advantage of the two steps (smoothing polygons and eliminating small polygons) was that they eliminated image quality effects (e.g., salt-and-pepper effect) in the classification results.

Identifying Which Line for the Tree Line Ecotone Was Automatically Extracted from Landsat Imagery
The ecological transition between the subalpine forest and the alpine meadow is an ecotone associated with tree species line, tree line, and timber line. However, remote-sensing techniques (land cover classification and change detection) interpret an ecotone as a sharp boundary line. Therefore, it is important to identify which line was automatically extracted from Landsat imagery via the methodology used in this study. The line extracted automatically based on the NDVI was the timber line, which lies above the elevational limit of the closed forest ( Figure 8). The NDVI contributed to timber line extraction because it smoothed the individual tree regions of the tree line ecotone and the alpine meadow. This differed from extraction using individual bands of Landsat images, which were very patchy due to the heterogeneity inherent in the tree line ecotone ( Figure 5). Therefore, the NDVI reduced heterogeneity effects caused by individual trees and small tree islands in the tree line ecotone.
Spatial resolution was the primary reason that Landsat images extracted the timber line. Tree islands were not obviously different from the alpine meadow in these images because Landsat's 30 × 30 m spatial resolution was not detailed enough to capture small tree islands in the ecotone ( Figure 8). However, tree line islands were extracted from remote-sensing images with appropriately high spatial resolutions (Figure 7). The boundary of tree islands in the tree line ecotone is not an ecotone in the ecological sense, where the tree line is the elevational boundary between individual trees and tree islands in the tree line ecotone (Figure 1e). The actual ecological tree line (Figure 1e) could be distinguished from the elevational boundary between regions, with identifiable tree islands and regions of the alpine meadow. Regions with individual trees could not be distinguished from the alpine meadow based on the methodology developed in this study; therefore, the tree species line (Figure 1e) was not identified in our study area.

Climate Impact on Upward Tree Line Shift
The annual mean temperature increased by 1 • C from 1987 (18.23 • C) to 2017 (19.23 • C) in the study area. With this degree of climate change, the tree line shifted upward to a maximum of 100 m, with an average of 50 m ( Figure 10). This finding was consistent with previous research that found that a 1 • C increase in mean annual temperature was associated with an upward tree line shift of about 70 m [18]. The majority of tree lines in the study area advanced upward 1.  [10]. In the subarctic climate region, tree line advancement has been found to be 0.5 m per year in Finland [12] and 0.8 m per year in a humid continental climate region [18]. In our study area, few anthropogenic factors affect the tree line ecotone. Anthropogenic disturbances in the high-altitude region are limited to hiking tourists along a trail, field surveys by researchers, and an abandoned military building (Figures 8 and 10), all of which only affect a small region along the trail in the alpine meadow. The productivity of the alpine meadow is lower in these anthropogenically disturbed regions, indicated by a low NDVI in the peak growing season ( Figure 6, Landsat 5, 4 June 1994, NDVI). The long, narrow alpine meadow region in the southwest of the study area (Figure 10f) became increasingly fragmented from 1987 to 2018, which indicated a strong vertical tree line advancement.
The results of this study indicate a relatively rapid upslope tree line advancement under increasing temperatures in this subtropical region. This differs from findings in studies on large upward shifts occurring in heavily anthropogenically disturbed areas of tropical highlands [24] without reaching the maximum potential tree line position [17]. In addition, the area of the alpine meadow continued to increase from 2003 to 2012, while the temperature decreased slightly during the same period, which might have indicated lag effects (nine year lag) of increasing temperature on tree line advancement (Figure 9). However, the inconsistent temporal trend between temperature and area of the alpine meadow during 2012-2018 may have been caused by heavy snowfall in 2008 or it might imply that temperature is not a driving factor of tree line shifts (Figure 9). At the regional scale, environmental factors that influence alpine tree line position and dynamics are more complex, because such factors related to tree line position vary from site to site and from species to species [38] or are due to interactions between climate factors and local environmental factors [16,41]. In the literature, many studies have observed local variations of tree line upward shifting [10,12,34,35], even in the same study area [20]. The local variation phenomena for tree line position and upward shifting also indicates that temperature is no longer the primary environmental parameter that promotes tree line upward advancement. Besides temperature, the dominant drivers of upward tree line shift at the regional scale that have been proven in the literature are aspect [13,60], soil moisture [27,41], soil fertilization [93], wind exposure and wind speed [22,52,94,95], N (Nitrogen) deficiency [96], radiation stress [97], biotic interactions [98], species composition [12], surrounding vegetation shielding tree seedlings [99], drought [1,100], volcanic eruptions [101], fire [58,[102][103][104], grazing [105][106][107][108], land-use change [13,16,19,73], and land management [53,109].

Local Variations in Tree Line Response to Climate Change
Two classes of tree lines exist in the study area: the dwarf forest and the coniferous forest. Based on the results of tree line advancement from this research, it seems that the dwarf forest tree line is quite stable (Figure 10a,b) and the coniferous forest tree line is not (Figure 10d,e). The habitat of the subalpine dwarf forest (Figure 2b) is in areas with less sunlight, where coniferous species are less likely to survive [110]. More importantly, long-term exposure to winter wind (in the north or northwest direction) makes trees with needle leaves unable to survive (the subalpine dwarf forest northern and northwestern aspects in Figure 10a,b, respectively) [110]. Another characteristic of the subalpine dwarf forest habitat is the abrupt change in topography between it and the alpine meadow (i.e., the dwarf forest in the western part of the alpine meadow ( Figure 2) is located in the fault zone), which may prevent the dwarf forest tree line from shifting upward along the elevation gradient. Therefore, different regeneration stages of dwarf tree species were not visible along the elevation gradient, unlike the coniferous tree line ecotone. In addition, in subalpine dwarf forests, various species of trees are clustered together instead of scattered individual trees or isolated tree islands in preferable sites in the mosaic region between the alpine meadow and the subalpine coniferous forest. For the coniferous forest, the stable tree line was found to be in locations with steeper slopes (Figure 10c), while the highly advanced tree line was found to be in areas with shallower slopes (Figure 10d-i).
In addition to slope, elevation is a primary topographic factor that results in local variations of climate impact on tree line advancement. The coniferous forest tree line tended to be more stable at high elevations (about 2063 m a.s.l.; Figure 11c) than that at lower elevations (about 1963 m, a.s.l.; Figure 11d). Most of the current tree line is located along an elevation gradient between 2000 and 2050 m ( Figure 10). Our results indicated that slope is one of the main factors affecting local variations in tree line advancement, even in the face of ecological responses to climate change. Indeed, the slope where tree lines were stable was steeper than that with the advancing tree line on all aspects (Figure 11b). Nevertheless, the impact of slope on tree line advancement did interact with aspect to a degree (Figure 11b). It was weaker in the southeast, southwest, and west aspects than in the other aspects (Figure 11b). In addition, elevation was no longer considered one of the primary topographic factors influencing local variations of tree line advancement in the southeast and east aspects (Figure 11d).

Conclusions
The main conclusions of this study are as follows. (1) The methodology based on the LISA concept has good potential for identifying ecological tree lines; we determined that the line extracted from Landsat images was actually the timber line. (2) The NDVI is better for tree line extraction than single bands from Landsat images, and the optimal season for tree line extraction in alpine meadows is the nongrowing season. (3) Most tree lines moved upward by 50 m under 1 • C warming during 1987-2018 with local variations. (4) Tree lines are more likely to move upward with increasing temperature on shallower slopes than on steeper slopes, and the impact of slope on local variations of tree line advancement decreases on the southeast, southwest, and west aspects. (5) Tree lines are more stable at higher elevations than at lower elevations, except for those on the east and southeast aspects.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/18/2890/s1. Figure S1. Cloud and cloud shadow effects on automatic tree line extraction from Landsat images; Supplement A: Tree line extraction by NDVI from all Landsat images; Supplement B: Accuracy assessment with reference to RGB image with a spatial resolution of 0.53 × 0.53 m; Supplement C: Best extraction of tree line for each year from Landsat images; Supplement D: Comparison of tree line extraction from Landsat and RGB images with a 16.9 m spatial resolution; Supplement E: Sample statistics for testing topographic effects, Python script, ArcToolbox.
Author Contributions: D.X. contributed by coming up with the initial ideas, developing the methodology, writing the Python script, analyzing the data, and writing the manuscript. Q.G. contributed by downloading and pre-processing images, improving the methodology and Python script, and revising the manuscript. C.J. and Z.X. contributed by providing information of the study area, explaining the research results, and writing the discussion. X.X. guided the organization of initial ideas and research directions. All authors have read and agreed to the published version of the manuscript.