Integrating Remote Sensing and Spatiotemporal Analysis to Characterize Artiﬁcial Vegetation Restoration Suitability in Desert Areas: A Case Study of Mu Us Sandy Land

: One of the major barriers to hindering the sustainable development of the terrestrial environment is the desertiﬁcation process, and revegetation is one of the most signiﬁcant duties in anti-desertiﬁcation. Desertiﬁcation deteriorates land ecosystems through species decline, and remote sensing is becoming the most effective way to monitor desertiﬁcation. Mu Us Sandy Land is the ﬁfth largest desert and the representative area under manmade vegetation restorations in China. Therefore, it is essential to understand the spatiotemporal characteristics of artiﬁcial desert transformation for seeking the optimal revegetation location for future restoration planning. However, there are no previous studies focusing on exploring regular patterns between the spatial distribution of vegetation restoration and human-related geographical features. In this study, we use Landsat satellite data from 1986 to 2020 to achieve annual monitoring of vegetation change by a threshold segmentation method, and then use spatiotemporal analysis with Open Street Map (OSM) data to explore the spatiotemporal distribution pattern between vegetation occurrence and human-related features. We construct an artiﬁcial vegetation restoration suitability index (AVRSI) by considering human-related features and topographical factors, and we assess artiﬁcial suitability for vegetation restoration by mapping methods based on that index and the vegetation distribution pattern. The AVRSI can be commonly used for evaluating restoration suitability in Sandy areas and it is tested acceptable in Mu Us Sandy Land. Our results show during this period, the segmentation threshold and vegetation area of Mu Us Sandy Land increased at rates of 0.005/year and 264.11 km 2 /year, respectively. Typically, we found the artiﬁcial restoration vegetation suitability in Mu Us area spatially declines from southeast to northwest, but eventually increases in the most northwest region. This study reveals the revegetation process in Mu Us Sandy Land by ﬁguring out its spatiotemporal vegetation change with human-related features and maps the artiﬁcial revegetation suitability.


Introduction
Drylands occupy 41% of the land surface on earth and their existence influences the life of 38% of the world's people and imposes a great threat to global land sustainability [1].
Cover (CCI-LC) and Moderate Resolution Imaging Spectroradiometer (MODIS), have been commonly used in monitoring the dynamics of the land surface [18][19][20][21]. Although these data products offer quality information in data-limited areas, the medium accuracy still cannot adapt to a smaller regional extent. Moreover, its insufficient time span brings a barrier to long-term earth surface monitoring.
In the past, there have mainly been two types of research on land surface monitoring by remote sensing in the Mu Us Sandy Land: (1) to directly classify the land cover classes, and then use these multi-temporal classification results for comparative analysis [4]; and (2) to compare spectral indices such as NDVI and Desertification Vulnerability Index (DVI) to reflect the gradual surface change. Both methods can effectively reflect the land cover status of the Mu Us Sandy Land to analyze and evaluate the local vegetation restoration [9]. However, vegetation and desertification indexes can only be used to reflect changes in vegetation growth or desertification level but its "non-vegetation" definition makes it difficult to use for spatiotemporal analysis or aggregation. In addition, previous studies have shown that the human and natural environmental system of the Mu Us Sandy Land has been undergoing frequent changes in recent years, which complicates the situation of restoring vegetation [22], but most existing remote sensing work on its land cover change can only cover several periods in past 30 to 40 years. These studies can generally reflect its overall changes in a long-term specific span, but for the short-term, such as on an annual or monthly basis, it can only provide a limited effect. Therefore, there is still a need to monitor the land surface on a higher time resolution for further understanding of the dynamic vegetation area in Mu Us Sandy Land.
In recent years, researchers have explored the causes of vegetation restoration in the Mu Us Sandy Land at multiple levels. Liu has used panel data to analyze and evaluate the relative role of climate, socio-economic conditions, and policy effects at a province level, showing the crucial positive effect of policy implementations on vegetation and ecological restoration [23,24]. Xiu and his group quantitatively analyzed the effect of the ecological engineering projects on the vegetation restoration of Mu Us and emphasized the importance of human activities to the growth of regional vegetation under the background of climate change [25]. Sun used linear regression methods to compare the vegetation recovery effects of climate factors and afforestation and found afforestation activities are more effective [26]. Most views believe that human interference imposed a significant effect on vegetation restoration, but there is no spatially explicit evidence to demonstrate the sensitivity of vegetation restoration to human activities, because most studies have made the "restoration process" a confirmed assumption and related it with vegetation changes [6,24,27]. Therefore, it is particularly important to explore the relationship between spatiotemporal vegetation changes and explicit geographical features of human activities.
A series of previous studies comprehensively assessed the vegetation suitability in Mu Us Sandy Land and worldwide dry areas by modeling and remote sensing methods. For example, Wang and his team integrated hydrological and ecological models to evaluate the sustainable vegetation restoration Loess Plateau of China, based on natural factors such as rainfall and soil [28]. McVicar and his group simulated vegetation suitability in Loess Plateau by considering land uses, topographical characteristics, and water use [29]. Peng and his research team exploited hidden vegetation suitability at the vegetation scale in terrestrial ecosystem models to apply to revegetation programs under diverse climate change scenarios [30]. In Africa, remote sensing-based vegetation indices were proven to be effective indicators to assess the reactions of terrestrial ecoregion drought [31]. However, to the best of our knowledge, there are no studies focusing on mapping vegetation suitability based on different revegetation convenience levels by considering the spatiotemporal distribution of human footprint features. The revegetation convenience refers to the different difficulties of construction in revegetation activities, such as different slopes causing different difficulties in restoration measurements. To bridge this knowledge gap, this study combines methods of satellite remote sensing and spatial and temporal analysis to explore the impact of human-related features on vegetation change.
There are three objectives on which our research is focused: (1) to achieve annual vegetation and non-vegetation classification in Mu Us Sandy Land; (2) to explore the relationship between spatial vegetation distribution patterns and human-related footprint features; (3) to construct a method to assess the artificial vegetation restoration suitability in desert areas and apply it in Mu Us case. Therefore, we designed our study contents based on the above three purposes, firstly, we exploit an Otsu-based threshold segmentation method to automatically extract annual vegetation pixels based on an improved Enhanced Normalized Difference Vegetation Index (ENDVI), for providing a data base of image classification and following spatiotemporal analysis. Secondly, we selected the geographic information features related to human activities and used the spatiotemporal analysis methods to perform vegetation distribution in these features to explore the regular pattern of vegetation change in Mu Us Sandy Land. Thirdly, we construct a novel AVRSI and assess the vegetation suitability of manmade restoration in the Mu Us Sandy Land by analyzing the influence of the dominant characteristics of human activities on the distribution of vegetation. Our research provides a scientific reference for making more sustainable vegetation restoration policies or projects in the Mu Us Sandy Land and other deserts.

Study Area
Mu Us Sandy Land, lying on the Ordos Plateau in China, administratively belonged to the Inner Mongolia Autonomous Region, Shaanxi Province, and the Ningxia Hui Autonomous Region. Its district is mainly located between 107 • 21 E and 111 • 30 E, and 37 • 27 N and 39 • 22 N, with an elevation from 0 m to 1614 m a.s.l ( Figure 1) and an area of about 42 thousand km 2 [6]. Continental semi-arid climate is its representative climate characteristic. The annual temperature ranges from 6 • C to 8.5 • C, with monthly mean temperatures of 22 • C in July and −11 • C in January. Annual precipitation varies from 440 mm in the southeast to 250 mm in the northwest, which is up to 80%, particularly from June to August.
to explore the impact of human-related features on vegetation change.
There are three objectives on which our research is focused: (1) to achieve annua vegetation and non-vegetation classification in Mu Us Sandy Land; (2) to explore the re lationship between spatial vegetation distribution patterns and human-related footprin features; (3) to construct a method to assess the artificial vegetation restoration suitability in desert areas and apply it in Mu Us case. Therefore, we designed our study content based on the above three purposes, firstly, we exploit an Otsu-based threshold segmenta tion method to automatically extract annual vegetation pixels based on an improved En hanced Normalized Difference Vegetation Index (ENDVI), for providing a data base o image classification and following spatiotemporal analysis. Secondly, we selected the ge ographic information features related to human activities and used the spatiotempora analysis methods to perform vegetation distribution in these features to explore the regu lar pattern of vegetation change in Mu Us Sandy Land. Thirdly, we construct a nove AVRSI and assess the vegetation suitability of manmade restoration in the Mu Us Sand Land by analyzing the influence of the dominant characteristics of human activities on th distribution of vegetation. Our research provides a scientific reference for making mor sustainable vegetation restoration policies or projects in the Mu Us Sandy Land and othe deserts.

Study Area
Mu Us Sandy Land, lying on the Ordos Plateau in China, administratively belonged to the Inner Mongolia Autonomous Region, Shaanxi Province, and the Ningxia Hui Au tonomous Region. Its district is mainly located between 107°21′E and 111°30′E, and 37°27′N and 39°22′N, with an elevation from 0 m to 1614 m a.s.l ( Figure 1) and an area o about 42 thousand km 2 [6]. Continental semi-arid climate is its representative climate char acteristic. The annual temperature ranges from 6 °C to 8.5 °C, with monthly mean tem peratures of 22 °C in July and −11 °C in January. Annual precipitation varies from 440 mm in the southeast to 250 mm in the northwest, which is up to 80%, particularly from Jun to August.  The main type of land surface with vegetation is sandy grassland and the main plantation for cultivation are corn, buckwheat, potato, etc. [4]. There are mainly four landscape types in Mu Us Sandy Land: sandstone hills, meadow steppe, and active and fixed sand dunes [32]. Several rivers lie across the southeast of the region and gather into the Yellow River. Lakes are mainly distributed in the internal land of the area, but most of them are sodic lakes with chloride, and only a few are freshwater lakes [9].

Data Preparations
The process of vegetation classification is conducted through the cloud-based platform Google Earth Engine (GEE) and Pixel Information Expert Engine (PIE) (https://engine. piesat.cn/ (accessed on 9 September 2021)). Our study period lasts for 35 years; therefore, the historical images for monitoring vegetation areas are derived from Landsat Thematic Mapper (TM), Landsat Enhanced Thematic Mapper Plus (ETM+), and Landsat Operational Land Imager (OLI) with 30 m spatial resolution. All Landsat Tier1 level data are created using the best available processing level for each scene. The Landsat images with rare clouds at the Tier1 level were collected based on time and space conditions. According to the processing results, there are no data available covering the study area before 1986, so the images selected in the 35 years of data from 1986 to 2020 were counted. Considering climate factors and period consistency, images synthesized every summer are used as a study database. Totally 2785 scenes are used in our image synthesis process. The satellite data sources, period of selected scenes, and the number of scenes for each year are described in Table S1 (Supplementary File).
Human information features include roads, railways, human settlements, service stations, and waterways, which are collected from OSM and updated to the 2020 version. We also selected slope as a topographical feature for mapping the artificial suitability of vegetation restoration, which is derived from the Shuttle Radar Topography Mission (SRTM) at a scale of 30 m.

Methods
This research consists of two main sections: annual vegetation detection by remote sensing and spatiotemporal analysis with human-related factors. The final aim of this study is to map artificial vegetation restoration suitability based on the two above methods. Figure 2 indicates the data input, classification algorithm, related spectral index, and specific steps of each section. tation for cultivation are corn, buckwheat, potato, etc. [4]. There are mainly four landscape types in Mu Us Sandy Land: sandstone hills, meadow steppe, and active and fixed sand dunes [32]. Several rivers lie across the southeast of the region and gather into the Yellow River. Lakes are mainly distributed in the internal land of the area, but most of them are sodic lakes with chloride, and only a few are freshwater lakes [9].

Data Preparations
The process of vegetation classification is conducted through the cloud-based platform Google Earth Engine (GEE) and Pixel Information Expert Engine (PIE) (https://engine.piesat.cn/ (accessed on 9 September 2021)). Our study period lasts for 35 years; therefore, the historical images for monitoring vegetation areas are derived from Landsat Thematic Mapper (TM), Landsat Enhanced Thematic Mapper Plus (ETM+), and Landsat Operational Land Imager (OLI) with 30 m spatial resolution. All Landsat Tier1 level data are created using the best available processing level for each scene. The Landsat images with rare clouds at the Tier1 level were collected based on time and space conditions. According to the processing results, there are no data available covering the study area before 1986, so the images selected in the 35 years of data from 1986 to 2020 were counted. Considering climate factors and period consistency, images synthesized every summer are used as a study database. Totally 2785 scenes are used in our image synthesis process. The satellite data sources, period of selected scenes, and the number of scenes for each year are described in Table S1 (Supplementary File).
Human information features include roads, railways, human settlements, service stations, and waterways, which are collected from OSM and updated to the 2020 version. We also selected slope as a topographical feature for mapping the artificial suitability of vegetation restoration, which is derived from the Shuttle Radar Topography Mission (SRTM) at a scale of 30 m.

Methods
This research consists of two main sections: annual vegetation detection by remote sensing and spatiotemporal analysis with human-related factors. The final aim of this study is to map artificial vegetation restoration suitability based on the two above methods. Figure 2 indicates the data input, classification algorithm, related spectral index, and specific steps of each section.

Annual Vegetation Detection by Remote Sensing
In this section, we first construct ENDVI by original NDVI, NDSI, and MNDWI and then use the Otsu method to classify vegetation and non-vegetation areas in Mu Us Sandy Land. Last, we made statistics on the analysis results. This analysis will provide a panel reference for past vegetation change in Mu Us Sandy Land and this result will be used for the following analysis of the relationship between vegetation change and human features in the next section.
After image pre-processing, such as radiation correction, image mosaic, and band synthesis, we calculate three spectral indices NDVI, Modified Normalized Difference Water Index (MNDWI), and Normalized Difference Soil Index (NDSI) that, respectively, characterize vegetation, water, and bare land. Then, we unified these three indexes as values from −1 to 1 and mask the area of identified water bodies. To expand the background value difference between vegetation and bare land, we subtract NDSI from the value of NDVI, constructing a new ENDVI threshold under the background of weakened bare soil. This Otsu-based segmentation method with spectral index is commonly used in earth surface detection, but a proper index for differentiating vegetation and non-vegetation is still under exploration in desert areas. Therefore, after the above processing of coupled indices, the difference of ENDVI between vegetation and non-vegetation will become larger, which makes it more adaptive in arid areas. Lastly, we normalize the ENDVI value from 0 to 1 before processing it into vegetation classification. The related spectral index equations and derivation process of the ENDVI formula are shown in the following: Among them, Red, NIR, and SWIR refer to the red band, near-infrared red band, and shortwave infrared band of satellite images, respectively.
The Otsu method is a nonparametric and unsupervised method of automatic threshold selection for picture segmentation, which is used to maximize the difference between vegetation and non-vegetation. Assuming that image pixels can be divided into two parts based on distinct gray levels: background and objects. When the variance of the gray value between the objects and the background reaches the maximum, the threshold can be the optimal segmentation threshold. Currently, the difference between the target and the background is the largest, and the segmentation is the most effective [33]. It can be briefly introduced that: Assuming that the gray level of an image is L, which presents as a level range [0,1,2...L − 1]. The total number of pixels in the image is N, n i represents the number of pixels with gray level i, namely: P i represents the probability of the expected point with gray level i, namely: If the pixels are divided into two classes, object α 1 and background α 2 , by the threshold t, α 1 is composed of pixels whose gray value is between [0, t], and α 2 is composed of pixels whose gray value is between [t + 1, L − 1]. Then the probabilities of these two classes w 1 and w 2 are: The means of these two classes can be computed as: Among them: If σ 2 b is used to represent the between-class variance between the two classes, then it is: Let the t values be taken in sequence from the range of [0, L − 1], when σ 2 b becomes largest, the corresponding t value is the optimal threshold obtained by the Otsu algorithm [34,35].
The Otsu algorithm is to process the pixels of the entire area, and then calculate the optimal threshold, which would be ineffective if an imbalanced proportion of the object and background pixels occurs. In this study, to realize a more equal ratio of the object and background pixels, canny edge detection was performed, and buffer zones were identified based on the contour line before employing Otsu segmentation. NDVI value is a general indicator to demonstrate the greenness of vegetation visually and numerically on an observed land surface. In general, NDVI values below 0.1 signify a low level of vegetation cover, namely barren areas such as rock, sand, and snow. Moderate values from 0.2 to 0.3 correspond to shrub and grassland, and relatively high levels from 0.6 to 0.8 always indicate more bushy forests [36]. In this study, vegetation and non-vegetation area are divided by the Otsu ENDVI threshold. Therefore, to detect the distribution of vegetation as much as possible, we assume vegetation areas possess two attributes: the ENDVI value is greater than the Otsu ENDVI threshold and the original NDVI is over 0.1.

Spatiotemporal Analysis for Vegetation Suitability Mapping
By reviewing works of literature, we found that the exact actions of ecological restorations mainly include constructing windscreens to secure sand and planting shrubs and grasses [6], and the effective implementation of these measures is strongly related to the existence of transportation and settlement services and water resources; thus, we selected five human-related features to reflect the ecological restoration footprint features, including roads, railways, human settlements, service stations, and waterways. Euclidean distance refers to the shortest straight-line distance between pixels in the geographical study, which is commonly used for referring distance convenience between two locations [37,38]. Due to the limitations of the availability and format of analysis data, we adopted the Euclidean distance to conduct long time-series spatial statistics to explore the rules of spatial accessibility between each pixel and human-related features in the study area.
To explore the relationship between vegetation distribution and the above geographical features of human activities, we use the Euclidean distance method to analyze the average distance from each vegetation pixel to these features. First, we use the int function in ArcGIS Pro to convert human-related features in vector form into the same raster form. Then, we calculated the Euclidean distance from each grid to these human-related features. The selected human-related features for time-series analysis are only 2020-based data; therefore, if the analyzed period for the past is too long, these features might show a big difference from 2020-based data. Therefore, we chose the period of ten years only from 2011 to 2020 for analysis to avoid excessive feature changes from affecting the analysis.

Constructing an AVRSI
Typically, land suitability is defined in terms of specific activities or land uses [39][40][41]. Due to worldwide changing economic development and social construction, the theories, methods, and practices of land suitability assessment have made rapid progress, and the application ranges of land suitability evaluation are also constantly improving [41,42]. Based on previous studies, the main fields of land suitability evaluation application can be basically divided into four aspects: agricultural land, forestry and animal husbandry land, urban construction land, tourism development land, and reclamation and consolidation land [43][44][45][46][47].
We compute the AVRSI and generate a suitability map based on the following three principles by referring distance analysis of spatiotemporal vegetation change: (1) accessibility and rationality of human-related feature data-we select the geographical data which can be acquired and closely related to restoration activities; (2) correlation between distances from human-related features and vegetation suitability can be scaled-we assume that vegetation area is expanding in the study area and its growth will linearly decline with the distance from human-related features; and (3) the effect weights of each human-related features on vegetation restoration can be quantified.
In this study, we focus on combining the distribution of geographical features of ecological restoration activities and the law of land development and utilization to construct an AVRSI and map the suitability of artificial vegetation restoration, to achieve laborconsuming artificial production of vegetation in desert areas. The defined expression of AVRSI is described below: where n is the total number of selected human-related features; α i is the effect weights of each human-related feature; x i indicates the distance convenience for artificial restoration; T is the coefficient ratio of topographical factors on the targeted grid.

Mapping Artificial Vegetation Restoration Suitability
Based on the above AVRSI computing principles for mapping the suitability of anthropogenic vegetation restoration, we use the functions of GIS to spatially express the suitability of the Mu Us Sandy Land: (1) we compute the Euclidean distance of each grid for the five human-related features and use the Int function to define each grid value as the Euclidean distance value (five human-related features correspond to five sets of Euclidean distance values), and we use the natural breaks (Jenks) method to classify the distance convenience, which is arranged as very high, high, moderate, low, and very low (quantitatively defined as 9:7:5:3:1 in descending order); (2) we determine the weights of each characteristic factor through the buffer analysis of vegetation distribution, and then we define the ratio of their weights as the ratio of their rounded regression equation slopes, namely, railway: roads: waterways: human settlements: service station is 1:2:3:2:1. In addition, we define the slope grades into 5 levels based on the Manual of Detailed Geomorphological Mapping [48]: and 36 • -83 • , and then we define the suitability coefficient ratio of these five grades as 1:0.9:0.7:0.5:0.1 based on our understanding of the influence of slope on human construction activities; and (3) we first calculate the weighted distance from each grid to these five human-related features through the raster calculator function, and then multiply the data of the modified weighted distance by the slope coefficient ratio of the corresponding grid, as shown in Equation (13).

ENDVI
The ENDVI values distribution and the corresponding number of pixels are shown in Figure S1 (Supplementary File) and Figure 3. From the two figures, we can find the overall ENDVI value in the Mu Us Sandy Land was relatively distributed between 0.4 and 0.6, and the number of pixels with ENDVI below 0.4 became less and the number of pixels with ENDVI over 0.6 became more. In addition, the ENDVI value points of which the count of pixels peaks every year are also increasing.
define the suitability coefficient ratio of these five grades as 1:0.9:0.7:0.5:0.1 based on our understanding of the influence of slope on human construction activities; and (3) we first calculate the weighted distance from each grid to these five human-related features through the raster calculator function, and then multiply the data of the modified weighted distance by the slope coefficient ratio of the corresponding grid, as shown in Equation (13).

ENDVI
The ENDVI values distribution and the corresponding number of pixels are shown in Figure S1 (Supplementary File) and Figure 3. From the two figures, we can find the overall ENDVI value in the Mu Us Sandy Land was relatively distributed between 0.4 and 0.6, and the number of pixels with ENDVI below 0.4 became less and the number of pixels with ENDVI over 0.6 became more. In addition, the ENDVI value points of which the count of pixels peaks every year are also increasing.       (Figure 4b). Regardless of the temporal change, the spatial pattern of vegetation distribution shows a general increase trend from southeast to northwest.

Spatiotemporal Analysis between Time-Series Vegetation Dynamic and Human-Related Factors
We computed the Euclidean distance from each grid to these human-related features, as shown in Figure 5, which we could use as spatial accessibility between vegetation and human-related features for restoration suitability. Mean Euclidean distances from the vegetation grid to human-related features are calculated annually as Figure 6.    (Figure 4b). Regardless of the temporal change, the spatial pattern of vegetation distribution shows a general increase trend from southeast to northwest.

Spatiotemporal Analysis between Time-Series Vegetation Dynamic and Human-Related Factors
We computed the Euclidean distance from each grid to these human-related features, as shown in Figure 5, which we could use as spatial accessibility between vegetation and human-related features for restoration suitability. Mean Euclidean distances from the vegetation grid to human-related features are calculated annually as Figure 6. From the results of the line chart in Figure 6, we found that the average Euclidean distance from the vegetation grid to these human-related features increased over the targeted ten years, showing that the overall boundaries of vegetation growth were expanding toward areas where it was not before, and the start points of these movements are human-related features. Therefore, the spatiotemporal characteristics of human ecological restoration measures are to restore from the areas closer to the human-related features to the area far from these features, which indicates that the difficulty of vegetation restoration increases with the distance between the targeted restoration location and the humanrelated features.  From the results of the line chart in Figure 6, we found that the average Euclidean distance from the vegetation grid to these human-related features increased over the targeted ten years, showing that the overall boundaries of vegetation growth were expanding toward areas where it was not before, and the start points of these movements are human-related features. Therefore, the spatiotemporal characteristics of human ecological restoration measures are to restore from the areas closer to the human-related features to the area far from these features, which indicates that the difficulty of vegetation restoration increases with the distance between the targeted restoration location and the human-related features. geted ten years, showing that the overall boundaries of vegetation growth were expanding toward areas where it was not before, and the start points of these movements are human-related features. Therefore, the spatiotemporal characteristics of human ecological restoration measures are to restore from the areas closer to the human-related features to the area far from these features, which indicates that the difficulty of vegetation restoration increases with the distance between the targeted restoration location and the humanrelated features. In addition, we analyze the relationship between vegetation distribution and humanrelated features in 2020 by establishing a buffer zone (10 km in total) with an interval of 1 km. The result is displayed in Figure 7. We found that as the distance from the geographic elements of human activities increased the proportion of vegetation grids in the total buffer grid decreased. In addition, we analyze the relationship between vegetation distribution and humanrelated features in 2020 by establishing a buffer zone (10 km in total) with an interval of 1 km. The result is displayed in Figure 7. We found that as the distance from the geographic elements of human activities increased the proportion of vegetation grids in the total buffer grid decreased. We assume that the distance from any locations in the Mu Us Sandy Land to these five geographical elements of human activities is linearly negatively correlated with the suitability of vegetation restoration. This result shows that human activities played a positive role in the process of vegetation restoration in Mu Us from the perspective of distance analysis. Based on the results shown in Figure 7, in 2020, we use the ratio of the slope of each trend of these features to represent the weight ratio of each element in mapping the suitability of vegetation restoration.
As one of the dominant topographical features, the slope reflects the steepness of the land surface, which has a significant impact on human-driven restoration activities [4,49,50]. As the slope in Mu Us sandy land shows nearly unchanged over the last 40 years, we count the slope values of each vegetation grid from 1986 to 2020 and calculate their averages for each year to explore the impact on vegetation distribution at Mu Us Sandy Land. The result is shown in Figure 8. We assume that the distance from any locations in the Mu Us Sandy Land to these five geographical elements of human activities is linearly negatively correlated with the suitability of vegetation restoration. This result shows that human activities played a positive role in the process of vegetation restoration in Mu Us from the perspective of distance analysis. Based on the results shown in Figure 7, in 2020, we use the ratio of the slope of each trend of these features to represent the weight ratio of each element in mapping the suitability of vegetation restoration.
As one of the dominant topographical features, the slope reflects the steepness of the land surface, which has a significant impact on human-driven restoration activities [4,49,50]. As the slope in Mu Us sandy land shows nearly unchanged over the last 40 years, we count the slope values of each vegetation grid from 1986 to 2020 and calculate their averages for each year to explore the impact on vegetation distribution at Mu Us Sandy Land. The result is shown in Figure 8.
As one of the dominant topographical features, the slope reflects the steepness of the land surface, which has a significant impact on human-driven restoration activities [4,49,50]. As the slope in Mu Us sandy land shows nearly unchanged over the last 40 years, we count the slope values of each vegetation grid from 1986 to 2020 and calculate their averages for each year to explore the impact on vegetation distribution at Mu Us Sandy Land. The result is shown in Figure 8. The spatial range of natural vegetation growth should be from the area with a low slope to the area with a higher slope. However, the slope of the vegetation grid increases over years from 1986 to 2020, following the changes in ENDVI threshold and vegetation area over the targeted years, which indicated that human-related restoration promoted The spatial range of natural vegetation growth should be from the area with a low slope to the area with a higher slope. However, the slope of the vegetation grid increases over years from 1986 to 2020, following the changes in ENDVI threshold and vegetation area over the targeted years, which indicated that human-related restoration promoted the vegetation expansion. From this perspective, we assume that the suitability of anthropogenic vegetation restoration decreases with the rising slope. This summary provides guidance for our subsequent vegetation restoration suitability.

Vegetation Restoration Suitability Mapping
The suitability results of anthropogenic vegetation restoration are shown in Figure 9. We can find that the overall spatial distribution of vegetation restoration suitability is decreasing from the southeast to the northwest, and after reaching the lowest point of suitability in the central-western area, it shows an upward trend in the northwest eventually. The probable causes of the restoration differences are attributed to the distribution of more dense waterways and flatter plains in the southeast, while the terrain in the central and western regions shows more large undulating slopes and is far from the water source. However, more dense road distribution and gentler slope lead to upward suitability in the most northwest region. the vegetation expansion. From this perspective, we assume that the suitability of anthropogenic vegetation restoration decreases with the rising slope. This summary provides guidance for our subsequent vegetation restoration suitability.

Vegetation Restoration Suitability Mapping
The suitability results of anthropogenic vegetation restoration are shown in Figure 9. We can find that the overall spatial distribution of vegetation restoration suitability is decreasing from the southeast to the northwest, and after reaching the lowest point of suitability in the central-western area, it shows an upward trend in the northwest eventually. The probable causes of the restoration differences are attributed to the distribution of more dense waterways and flatter plains in the southeast, while the terrain in the central and western regions shows more large undulating slopes and is far from the water source. However, more dense road distribution and gentler slope lead to upward suitability in the most northwest region. Integrating the annual changes in the vegetation distributions in the Mu Us Sandy Land with the regular patterns of artificial transformation found in result 3.3 that is, vegetation tends to be restored closer to waterways, roads, and railways, this suitability result is well-matched with the actual expectations, and it also shows that the vegetation suita- Integrating the annual changes in the vegetation distributions in the Mu Us Sandy Land with the regular patterns of artificial transformation found in result 3.3 that is, vegetation tends to be restored closer to waterways, roads, and railways, this suitability result is well-matched with the actual expectations, and it also shows that the vegetation suitability index designed in Section 2.3.3 is applicable to the progress of artificial vegetation restoration in the Mu Us Sandy Land over recent years.

Anti-Desertification Process in Mu Us Sandy Land
Before the 21st century, human activities, such as agricultural and animal husbandry activities, and mineral resources development, strongly disturbed the natural environment in the Mu Us region, leading to the instability of the "greenness degree" and vegetation area [22,51]. Since the beginning of this century, when human activities have become ecological restoration activities, desertification showed the characteristics of reversal in this desert. Under the governmental implementations of the multipolicy, such as the "Grain for Green Project", "The Three Northern Regions (northeastern, northwestern, and northern China) Shelter", and "mine ecological restoration", the vegetation area and the "greenness degree" have significantly increased in Mu Us Sandy Land [24,51,52]. This environmental benefit can be indicated from the result of Figure 4a and Figure S2 by our remote sensing work, the upward trend of the segmentation threshold and vegetation area indicates that the "greenness degree" and vegetation area in the Mu Us Sandy Land are continuously developed, which is basically consistent with the previous studies by many researchers [25,53,54].
The driving factors for the vegetation restoration in Mu Us Sandy Land mainly include natural driving factors and anthropogenic driving factors. The phenomenon of desertification existed throughout the quaternary period; thus, the impact of human activities on the landscape change before the modern era only accounted for a small effect of the entire desertification process [55]. Arid and windy extreme events deteriorate the stability of ecosystems, exacerbating soil barrenness, wind erosion, and vegetation degeneration [56], but the previous study clarified the climate factor is one of the minimally significant factors influencing the ecological restoration process in the Mu Us desert areas [24], which is the main reason our study exclude climate change factor as consideration for suitability mapping. On the other hand, along with the increasing frequency of human activities, the improvement of productivity levels, and the expansion of the scope and content of agricultural activities in arid areas, anthropogenic factors have also become crucial factors in the changes of desert vegetation [57]. These causes of vegetation degradation in arid areas include excessive grazing, deforestation, unreasonable farming, and overuse of water resources caused by the increasing population and inordinate utilization of natural resources. However, the impact of anthropogenic factors on desert vegetation can be bidirectional. We found that the changes in the "greenness degree" and vegetation area in Mu Us show a time-series fluctuating increase. The fluctuation shows that the unreasonable socioeconomic activities that still exist in the Mu Us area, mainly occurring the unreasonable agricultural and animal husbandry development, indiscriminate logging, deforestation, intensive mining activities, and possible sand dune movement, while the results show an overall increasing trend indicates that ecological restoration activities dominate the change of Mu Us vegetation [51].

Future Restoration
Combining the suitability mapping results in Figure 10 with the vegetation cover in 2020 (the most recent targeted year in this study), we can obtain the distribution of the difficulty level of implementation of anthropogenic vegetation restoration, as the results shown in Figure 10.

Future Restoration
Combining the suitability mapping results in Figure 10 with the vegetation cover in 2020 (the most recent targeted year in this study), we can obtain the distribution of the difficulty level of implementation of anthropogenic vegetation restoration, as the results shown in Figure 10. If the vegetation restoration measures in the Mu Us Sandy Land continue to be promoted in the future, the non-vegetation areas in the southeast and most northwest areas show better adaptability, indicating that the obstacles to vegetation restoration are relatively small and it would be more efficient to achieve significant results there. Therefore, these areas with higher suitability should be restored to vegetation first, rather than those areas with lower suitability.
Our study outcomes are aiming to save labor and economical costs during future ecological; however, due to the lack of data for specific costs for restoration processes, such as transportation of plants and irrigation, we cannot give a quantitative analysis for evaluation of its economic feasibility. Therefore, making a more comprehensive assessment of artificial restoration suitability by collecting panel data for specific restoration costs is significant in our future work.

Limitations and Future Work
This study uses the, x i , α i , and T (Equation (13)) as influencing factor coefficients to a novel vegetation suitability index for the anthropogenically transforming desert, which acts as a novel method to achieve mapping vegetation suitability in the Mu Us region based on results of remote sensing interpretation with an originally designed ENDVI and spatiotemporal analysis.
However, there exist the following limitations in this study: (1) the spatial resolution of the remote sensing data used in this study is 30 m, which indicates that there are still some details that need to be explained more accurately; (2) due to the availability of data, the selected anthropogenic geographical features have relatively limited influence on the distribution of vegetation restoration and also probably change in future, causing the impracticality of suitability mapping; (3) meanwhile, the overlapping effects of multiple factors cannot be quantitatively valued; (4) although the Mu Us Sandy Land is a representative study area due to its special agropastoral location and on-going vegetation restoration activities, it is inevitable that regional differences still exist. Typically, since a previous study has explained the climate factor accounts for limited weight in Mu Us vegetation change [24], this study is a special case for transforming desert which is less influenced by climate conditions; (5) since the restoration measures and financial costs of different vegetation species are different, the lack of subclassification of vegetation classes will lead to a certain difference between the modeled suitability and the actual suitability.
We are supposed to solve the above limitations from the following aspects: (1) to use satellite images with a higher spatial resolution for monitoring; (2) to investigate more data sources and update them for optimizing the completeness of selected human-related features and topographical factors, and if deems necessary, consider more additional other related factors; (3) to use proper weight analysis methods to rank vegetation change driven factors; (4) when applying the vegetation suitability mapping methods to tropical and subtropical regions, the climatic and geomorphological characteristics, such as surface temperature, precipitation, soil moisture content, soil texture, and wind force, are supposed to be additionally considered; and (5) to subdivide vegetation types by higher resolution remote sensing images or ground observation, and use this more detailed vegetation classification results to comprehensively make artificial restoration suitability for different vegetations.

Conclusions
This study uses a remote sensing method based on an automatic threshold segmentation algorithm to monitor annual vegetation change in Mu Us Sandy Land from 1986 to 2020 by an introduced novel enhanced normalized vegetation index. Based on analysis of the threshold segmentation and vegetation classification results, we have found, albeit with some fluctuations, that both the "greenness degree" and vegetation area show an obvious overall upward trend, especially after 2000. Then, we use distance analysis methods to investigate regular patterns between time-series vegetation distribution and human activity footprint features. These features include roads, railways, human settlements, service stations, waterways, and topographic factors. The results show that the vegetation restoration is expanding from the area with a lower slope near human-related features to the area with a higher slope far from these features, and this pattern effectively explains the vegetation restoration did cause by human activities from the perspective of spatiotemporal analysis. Moreover, we exploit the buffer zone function to analyze the relationship between vegetation distribution and human-related features, of which results indicate that the area portion of vegetation occurrence decreases with the distance from human-related features to the target location, and waterways, roads, and human settlements act as more influential geospatial features. Last, we construct an AVRSI for manmade desert transformation to map vegetation suitability based on distance analysis of human-related features and topographical factors. In the case of the Mu Us Sandy Land, we find artificial restoration vegetation suitability spatially decreases from southeast to northwest, but it shows a considerably increasing trend in the most northwest region. This study offers an effective means for quantitatively evaluating the resistance strength and expected effect of manmade desert transformation and provides data support and reference suggestions for relevant governments and organizations to formulate future policies and plans for artificial desert transformation.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/rs14194736/s1, Table S1. Data sources, period of selected scenes, and the number of scenes for each year from 1986 to 2020; Table S2. Annual segmentation thresholds for Otsu method; Figure S1. Relationship between ENDVI value and number of corresponding pixels; Figure S2. Annual detected vegetation in the Mu Us region from 1986 to 2020.

Data Availability Statement:
The data presented in this study are cited within the article.