Fine-Resolution Mapping of Soil Total Nitrogen across China Based on Weighted Model Averaging

Accurate estimates of the spatial distribution of total nitrogen (TN) in soil are fundamental for soil quality assessment, decision making in land management, and global nitrogen cycle modeling. In China, current maps are limited to individual regions or are of coarse resolution. In this study, we compiled a new 90-m resolution map of soil TN in China by the weighted summation of random forest and extreme gradient boosting. After harmonizing soil data from 4022 soil profiles into a fixed soil depth (0–20 cm) by equal area spline, 18 environmental covariates were employed to characterize the spatial pattern of soil TN in topsoil across China. The accuracy assessments from independent validation data showed that the weighted model averaging gave the best predictions with an acceptable R2 (0.41). The prediction map showed that high-value areas of soil TN were mainly distributed in the eastern Tibetan Plateau, central Qilian Mountains and the north of the Greater Khingan Range. Climate factors had a considerable influence on the variation of the soil TN, and land-use types played a pivotal part in each climate zone. This high-resolution and high-quality soil TN data set in China can be very useful for future inventories of soil nitrogen, assessments of soil nutrient status, and management of arable land.


Introduction
Soil nitrogen is a common macronutrient needed by plants for vegetation growth [1].Changes in the soil nitrogen, an important part of the nitrogen cycle, can affect the stability and sustainability of global ecosystems [2,3].Via the biological processes of nitrification and denitrification, excessive soil nitrogen can diffuse into the atmosphere as a greenhouse gas (N 2 O and NO) [4,5].Dissolved nitrogen may seep into water bodies where it may contribute to eutrophication and trigger other ecosystem changes and responses [6].Therefore, an explicit understanding of soil nitrogen content and its spatial variation are of great importance for soil quality assessment, decision making in land management, and global nitrogen cycle modeling, which have grave impacts on various global issues.
Maps of soil attributes are conventionally complied by filling polygons with the average attribute values of a particular soil type or land use [7,8].The procedure is laborious, time-consuming, expensive,

Soil Data
The soil data are from the Second National Soil Survey [30][31][32][33][34][35], which was conducted from 1979 to 1985.While the geographic locations of the data points were recorded, the specific spatial coordinates are not recorded, which means that errors are inevitable.Despite this, we have used this dataset in this study because of its wide spatial coverage and range of attributes.The information about soil TN in this database was determined using the semi-micro Kjeldahl digestion procedure [36].
The contents of soil TN were recorded by soil genetic horizons, so we derived in the uppermost 20 cm soil TN using the equal-area smoothing spline function, which proposed by Bishop et al. [37], debuted in soil depth functions by Malone et al. [38], and widely employed in several studies [17,19].
If we suppose that a soil profile has n horizons, and that x is the sampling depth, then for layer i (i = 1, 2, . . ., n), the upper and lower boundaries are x i−1 and x i , and f (x) is depth function between the depth x and the soil TN measurement.The measurement of TN in layer i is, therefore y i and is calculated as: where f i is the f (x) mean value of layer i and f i = x i x i−1 f (x)d(x)/(x i −x i−1 ), e i is the measurement error and e i = 0. We can determine the formula for f (x) by minimizing as follows: where λ is the spline smoothing parameter.The smaller the value of λ, the tighter the fit between the spline and the input values.The optimal λ value is obtained from the smallest mean squared error.
After spline interpolating with Spline Tool v2 (www.asris.csiro.au),4022 soil TN records at 0-20 cm were obtained from soil profiles data (Figure 1).In this study, the soil TN data were logarithmically transformed so that the samples showed less deviation from the normal distribution.Later, the prediction results were back-transformed to the original values with an antilogarithm.where  is the spline smoothing parameter.The smaller the value of  , the tighter the fit between the spline and the input values.The optimal  value is obtained from the smallest mean squared error.
After spline interpolating with Spline Tool v2 (www.asris.csiro.au),4022 soil TN records at 0-20 cm were obtained from soil profiles data (Figure 1).In this study, the soil TN data were logarithmically transformed so that the samples showed less deviation from the normal distribution.Later, the prediction results were back-transformed to the original values with an antilogarithm.

Environmental Covariates
In this study, we collected 18 environmental variables from available remote-sensing and ancillary data that may have influenced the distribution of soil TN (Table 1).
Terrain information at a resolution of 90 m was obtained from the shuttle radar topographic mission (SRTM, https://www2.jpl.nasa.gov/srtm/).Other terrain attributes, including slope, aspect, curvature, terrain ruggedness index (TRI), topographic wetness index (TWI), and multi-resolution valley-bottom flatness (MrVBF) were calculated from DEM in SAGA (System for Automated Geoscientific Analyses) GIS [39].The main terrain features were displayed in Figure S1 as well as the climate zones.
We obtained values of net primary productivity (NPP) and normalized difference vegetation index (NDVI) from the Global Production Efficiency Model (GloPEM) and Global Inventory Modeling and Mapping Studies (GIMMS) datasets, respectively, held by the Global Land Cover Facility at the University of Maryland [40,41].For each grid point, the mean values of NPP and NDVI were calculated from 1981 to 1985 and used then as variables to aid the predictions.

Environmental Covariates
In this study, we collected 18 environmental variables from available remote-sensing and ancillary data that may have influenced the distribution of soil TN (Table 1).
Terrain information at a resolution of 90 m was obtained from the shuttle radar topographic mission (SRTM, https://www2.jpl.nasa.gov/srtm/).Other terrain attributes, including slope, aspect, curvature, terrain ruggedness index (TRI), topographic wetness index (TWI), and multi-resolution valley-bottom flatness (MrVBF) were calculated from DEM in SAGA (System for Automated Geoscientific Analyses) GIS [39].The main terrain features were displayed in Figure S1 as well as the climate zones.
We obtained values of net primary productivity (NPP) and normalized difference vegetation index (NDVI) from the Global Production Efficiency Model (GloPEM) and Global Inventory Modeling and Mapping Studies (GIMMS) datasets, respectively, held by the Global Land Cover Facility at the University of Maryland [40,41].For each grid point, the mean values of NPP and NDVI were calculated from 1981 to 1985 and used then as variables to aid the predictions.Land surface temperature, day time (LSTD), land surface temperature, night time (LSTN), and evapotranspiration (ET) were collected from two Moderate Resolution Imaging Spectroradiometer (MODIS) datasets (https://lpdaac.usgs.gov/),namely MOD11A2 and MOD16A3.To correct the non-negligible difference between the acquisition time of the profiles and satellite data, fusion data was calculated using collocated cokriging [42], with the ground meteorological station data for 1980 to 1985, downloaded from the National Meteorological Information Center (http://data.cma.cn/), as the main variable, and the MODIS data sets as auxiliary variables.Data for the mean annual temperature (MAT) and mean annual precipitation (MAP) were derived from the Resource and Environment Data Cloud Platform (REDCP, http://www.resdc.cn/).The mean annual solar radiation (MASR) data can be acquired from the National Earth System Science Data Sharing Infrastructure (http://www.geodata.cn).Data for land use and vegetation types in 1980 were also obtained from the REDCP.These two datasets contain six land-use types (25 subclasses) and 10 vegetation groups (54 subgroups).The Chinese 1:1,000,000 soil map was published after China's Second National Soil Inventory [43].We converted 61 soil groups based on the Genetic Soil Classification of China (GSCC) system into 13 soil types based on the World Reference Base for Soil Resources (WRB) system [44,45].Some groups with few samples and similar characteristics had been integrated into a great group to ensure that each category has enough samples to train a model.Finally, the TN contents at the great group level were summarized in Figure 2.
All of these raster layers were resampled to the 90 m grid.

Model Development
Two machine-learning approaches named extreme gradient boosting (XGBoost) and random forest (RF) were used to produce the predictions and their lower and upper 90%.These two predictions were then ensembled by WMA.

Machine-Learning Approaches
The Random Forest creates a series of randomly generated classification and regression trees based on a bagging algorithm [46].The basic idea of RF is that each tree grows independently by randomly sampling replacements from an original training set, with the unselected samples considered as out-of-bag data for testing the accuracy of the predictions.Then m covariates (m is less than the total number of the covariates) are randomly selected to split the tree at each node in an optimal way by minimizing the variance of the split groups.Each tree of the model will be at its maximum growth with no pruning, and the final prediction is the average of all fitted trees.
XGBoost is based on a gradient-boosting algorithm [47].This method attaches importance to incorrect predictions done by prior learners and converts a series of weak and inaccurate learners to a strong predictor by weighted summation.In XGBoost, the objective function that needs to be minimized is equal to the sum of two parts, the training loss and regularization loss.Regularization is a function defined by the complexity of a tree that can avoid overfitting; shrinkage and column subsampling also play the same role.

Weighted Model Averaging
Each machine learning algorithm has its strengths and weaknesses [25].Weighted model averaging is an approach combining the outcomes across different contributor models by the weighted summation.The principle of model averaging suggests the new prediction value is at least as good as any of the individual prediction value [25].The new prediction value of WMA y ˆ is calculated as:

Model Development
Two machine-learning approaches named extreme gradient boosting (XGBoost) and random forest (RF) were used to produce the predictions and their lower and upper 90%.These two predictions were then ensembled by WMA.

Machine-Learning Approaches
The Random Forest creates a series of randomly generated classification and regression trees based on a bagging algorithm [46].The basic idea of RF is that each tree grows independently by randomly sampling replacements from an original training set, with the unselected samples considered as out-of-bag data for testing the accuracy of the predictions.Then m covariates (m is less than the total number of the covariates) are randomly selected to split the tree at each node in an optimal way by minimizing the variance of the split groups.Each tree of the model will be at its maximum growth with no pruning, and the final prediction is the average of all fitted trees.
XGBoost is based on a gradient-boosting algorithm [47].This method attaches importance to incorrect predictions done by prior learners and converts a series of weak and inaccurate learners to a strong predictor by weighted summation.In XGBoost, the objective function that needs to be minimized is equal to the sum of two parts, the training loss and regularization loss.Regularization is a function defined by the complexity of a tree that can avoid overfitting; shrinkage and column subsampling also play the same role.

Weighted Model Averaging
Each machine learning algorithm has its strengths and weaknesses [25].Weighted model averaging is an approach combining the outcomes across different contributor models by the weighted summation.The principle of model averaging suggests the new prediction value is at least as good as any of the individual prediction value [25].The new prediction value of ŷWMA is calculated as: where M is the number of predictors, ŷmi is the prediction from the predictor m, and w m is the weight attributed to that model.In this study, M = 2, the value of ŷWMA is the weighted summation of the predictions from RF and XGBoost: ŷRF and ŷXGBoost .The least absolute shrinkage and selection operator (LASSO) regression was used to solve for the parameters: w RF and w XGBoost .The LASSO algorithm is better than the unpenalized least squares algorithms as it can enforce regularity rather than simply lead to putting all weight on the most complicated predictor, which had been proved by Hansen [18].

Model Calibration and Validation
After partitioning the soil samples into calibration data (90%, 3620) and validation data (10%, 402), we constructed prediction models using calibration data and their corresponding 18 covariates, and then the model performance was evaluated with 402 validation data.
For validation purposes, three indices including the determination coefficient (R 2 ), the root mean squared error (RMSE) and the mean error (ME) were selected to estimate the prediction accuracy by the aforementioned models.Their formulae are listed as follows.
where y is the mean of the measured data (equal-area spline interpolated values in our case), and y i and f i are measured and predicted values for sample i (i = 1, 2, . . ., n) respectively.

Uncertainty Assessment
To measure and quantify the uncertainty associated with the predictions, we generated 50 calibration datasets based on 50 times non-parametric bootstrapping.We constructed 50 RF models and 50 XGBoost models based on the 50 calibration datasets in R software with Caret (Classification and Regression Training) package [48,49].Then the WMA assigned weight parameters to the RF and XGBoost model of each iteration using the LASSO algorithm, with the TN measurement values as the dependent variable and the TN prediction values from RF and XGBoost as independent variables.The 50 WMA models were constructed with glmnet package [50].
After applying these models on each pixel (90 × 90 m), we derived 50 maps of soil TN for each single machine learning method.We calculated the weighted summation of two maps (i.e. a RF map and a XGBoost map) obtained from one iteration and then generated 50 WMA maps.We took the average of the 50 soil TN maps as one final map of China.For clarity, the process of spatial modeling was presented in Figure S2.
For each grid point, we calculated the 90% confidence intervals (CIs) of the prediction values from each algorithm, which indicated that the true value of the TN content has 90% possibility within the interval between upper and lower CI limits.The 90% CIs is calculated in Equation ( 7): where P is the average soil TN content of the n times predictions, and in this study, n = 50.s is the standard deviation of n times predictions.Finally, we can calculate the uncertainty of our estimates in Equation ( 8): where CI upper and CI lower are the lower and upper bounds of CIs.

One-Way Analysis of Variance (ANOVA) Test
One-way ANOVA (analysis of variance) test was used to test if the value of R 2 differs significantly between XGBoost, RF, and WMA models.

Importance of Covariates
The results of the average relative importance of the covariates through 50 iterations of XGBoost and RF are shown in Figure 3.The climate variables played a vital role among 18 covariates.Precipitation had relative values of about 28%, followed by land surface temperature (26%), radiation (22%), temperature (22%) and ET (21%).While elevation (24%) was the foremost topographic variable.The NDVI (26%) and NPP (22%) were of almost equal importance.The vegetation and land-use types were much less important than the other factors, with relative values of less than 12%.The importance of soil types is 23%.The prediction models assigned the highest importance to precipitation, land surface temperature and NDVI that these factors were pivotal to shape the spatial distribution pattern of soil TN.where P is the average soil TN content of the n times predictions, and in this study, n = 50.s is the standard deviation of n times predictions.Finally, we can calculate the uncertainty of our estimates in Equation ( 8

One-Way Analysis of Variance (ANOVA) Test
One-way ANOVA (analysis of variance) test was used to test if the value of R 2 differs significantly between XGBoost, RF, and WMA models.

Importance of Covariates
The results of the average relative importance of the covariates through 50 iterations of XGBoost and RF are shown in Figure 3.The climate variables played a vital role among 18 covariates.Precipitation had relative values of about 28%, followed by land surface temperature (26%), radiation (22%), temperature (22%) and ET (21%).While elevation (24%) was the foremost topographic variable.The NDVI (26%) and NPP (22%) were of almost equal importance.The vegetation and land-use types were much less important than the other factors, with relative values of less than 12%.The importance of soil types is 23%.The prediction models assigned the highest importance to precipitation, land surface temperature and NDVI that these factors were pivotal to shape the spatial distribution pattern of soil TN.

Evaluation of the Approaches
The weight parameters assigned to RF and XGBoost models through 50 iterations are summarized in Figure S3.Table 2 showed the average values of validation indices through 50 iterations of XGBoost, RF, and WMA algorithms are presented in Table 2.A plot between the measured and predicted log-transformed TN values from the three models is presented in Figure 4. Table S1 showed the output of the ANOVA analysis and whether there is a significant difference between the three models' average R 2 .

Evaluation of the Approaches
The weight parameters assigned to RF and XGBoost models through 50 iterations are summarized in Figure S3.Table 2 showed the average values of validation indices through 50 iterations of XGBoost, RF, and WMA algorithms are presented in Table 2.A plot between the measured and predicted log-transformed TN values from the three models is presented in Figure 4. Table S1 showed the output of the ANOVA analysis and whether there is a significant difference between the three models' average R 2 .   1 The root mean squared error. 2 The mean error.The ME values in the three models were very small, which indicated that the predictions were roughly unbiased.The R 2 of the Random Forest was always 0.04 higher, and the RMSE was 0.02 g•kg −1 lower than those for the XGBoost method.Therefore, the RF models were always preferred and assinged bigger weight parameters than XGBoost in constructing a WMA model (Figure S3). Figure 4a showed that both RF and XGBoost overestimated predictions in low values of targets and underestimated high values.Compared with each other, RF and XGBoost produced better predictions for the data points with low (Figure 4b) and high (Figure 4c) TN contents, respectively.After the two models were weight-averaged, consequently, the model performance was best with R 2 of 0.41 and RMSE value of 1.15 g•kg −1 ; this showed that the WMA could explain 41% of the spatial variation in the TN contents.The results of one-way ANOVA showed (p = 0.000) that there was a statistically significant difference in the mean R 2 values between WMA, RF and XGBoost.
The three uncertainty results for each validation data with the 50 XGBoost, 50 RF, and 50 WMA models are shown in Figure 5.The mean uncertainty value of WMA is 0.05, much smaller than that of XGBoost (0.14) and RF (0.08).For each validation data, the finest result always derived from the WMA method, except for the two points emphasized with orange color in Figure 5, which demonstrated that there was a high probability that the prediction uncertainty could be reduced by using the WMA approach.The ME values in the three models were very small, which indicated that the predictions were roughly unbiased.The R 2 of the Random Forest was always 0.04 higher, and the RMSE was 0.02 g•kg −1 lower than those for the XGBoost method.Therefore, the RF models were always preferred and assinged bigger weight parameters than XGBoost in constructing a WMA model (Figure S3). Figure 4a showed that both RF and XGBoost overestimated predictions in low values of targets and underestimated high values.Compared with each other, RF and XGBoost produced better predictions for the data points with low (Figure 4b) and high (Figure 4c) TN contents, respectively.After the two models were weight-averaged, consequently, the model performance was best with R 2 of 0.41 and RMSE value of 1.15 g•kg −1 ; this showed that the WMA could explain 41% of the spatial variation in the TN contents.The results of one-way ANOVA showed (p = 0.000) that there was a statistically significant difference in the mean R 2 values between WMA, RF and XGBoost.
The three uncertainty results for each validation data with the 50 XGBoost, 50 RF, and 50 WMA models are shown in Figure 5.The mean uncertainty value of WMA is 0.05, much smaller than that of XGBoost (0.14) and RF (0.08).For each validation data, the finest result always derived from the WMA method, except for the two points emphasized with orange color in Figure 5, which demonstrated that there was a high probability that the prediction uncertainty could be reduced by using the WMA approach.

Mapping of Soil TN and Its Uncertainty
A national scale map of soil TN was constructed by WMA at 90 m resolution (Figure 6), and its increased details also were displayed at a gradually finer scale.This map clearly shows the total geographical distribution of soil TN.The TN content in the topsoil (0-20 cm) ranged from 0.13 to 9.92 g•kg −1 and averaged 1.648 g•kg −1 across China.The TN contents varied considerably and were distributed unevenly.The TN contents were high (>4.8g•kg −1 ) in eastern areas of the Tibetan Plateau, the central Qilian Mountains, and the northern section of the Greater Khingan Range, and were low (<0.8 g•kg −1 ) between 35°N and 42.5°N, in the Tarim Basin, Qaidam Basin and western Inner Mongolia Plateau, particularly in the desert areas.The TN contents were also low on the Loess Plateau and the North China Plain.We mapped the spatial distribution of the uncertainty based on the 50 WMA soil TN map (Figure 7).The uncertainty values were low in the southeastern area where there were many soil profiles.Equally, the uncertainty values were high in the western area because of the low density of

Mapping of Soil TN and Its Uncertainty
A national scale map of soil TN was constructed by WMA at 90 m resolution (Figure 6), and its increased details also were displayed at a gradually finer scale.This map clearly shows the total geographical distribution of soil TN.The TN content in the topsoil (0-20 cm) ranged from 0.13 to 9.92 g•kg −1 and averaged 1.648 g•kg −1 across China.The TN contents varied considerably and were distributed unevenly.The TN contents were high (>4.8g•kg −1 ) in eastern areas of the Tibetan Plateau, the central Qilian Mountains, and the northern section of the Greater Khingan Range, and were low (<0.

Mapping of Soil TN and Its Uncertainty
A national scale map of soil TN was constructed by WMA at 90 m resolution (Figure 6), and its increased details also were displayed at a gradually finer scale.This map clearly shows the total geographical distribution of soil TN.The TN content in the topsoil (0-20 cm) ranged from 0.13 to 9.92 g•kg −1 and averaged 1.648 g•kg −1 across China.The TN contents varied considerably and were distributed unevenly.The TN contents were high (>4.8g•kg −1 ) in eastern areas of the Tibetan Plateau, the central Qilian Mountains, and the northern section of the Greater Khingan Range, and were low (<0.8 g•kg −1 ) between 35°N and 42.5°N, in the Tarim Basin, Qaidam Basin and western Inner Mongolia Plateau, particularly in the desert areas.The TN contents were also low on the Loess Plateau and the North China Plain.We mapped the spatial distribution of the uncertainty based on the 50 WMA soil TN map (Figure 7).The uncertainty values were low in the southeastern area where there were many soil profiles.Equally, the uncertainty values were high in the western area because of the low density of We mapped the spatial distribution of the uncertainty based on the 50 WMA soil TN map (Figure 7).The uncertainty values were low in the southeastern area where there were many soil profiles.Equally, the uncertainty values were high in the western area because of the low density of the soil samples.The large tracts of desert land and high-attitude depopulated zone make it difficult to sample by soil surveyors, which thus lead to the large uncertainty in these regions when using DSM technology.It is also noteworthy that the uncertainty was greatest in the Kunlun Mountains, Altun Mountains and the Three-River Source region, the source of the Yangtze, the Yellow and the Lantsang Rivers.These areas have a complex and highly fragmented landscape structure which can hardly be detected by the sparse soil data used in this study.
Remote Sens. 2019, 11, x FOR PEER REVIEW 10 of 18 the soil samples.The large tracts of desert land and high-attitude depopulated zone make it difficult to sample by soil surveyors, which thus lead to the large uncertainty in these regions when using DSM technology.It is also noteworthy that the uncertainty was greatest in the Kunlun Mountains, Altun Mountains and the Three-River Source region, the source of the Yangtze, the Yellow and the Lantsang Rivers.These areas have a complex and highly fragmented landscape structure which can hardly be detected by the sparse soil data used in this study.

Soil TN of Different Soil Types and Land-Use Types
Soil TN contents at a depth of 0-20 cm under different soil types and land-use patterns are shown in Table S2 and Figure 8.The soil TN varied considerably across the soil groups and ranked in the following order: Histosols > Phaeozems > Chernozems > Luvisols > Cryosols > Kastanozems > Acrisols > Anthrosols > Cambisols > Calcisols > Solonchaks > Arenosols > Fluvisols.The TN content also showed dramatic differences between land-use types.At the great group level, TN contents were highest in the forest, followed by grassland, cropland and construction land, while lowest in unused land.In terms of the subgroups, there was a direct relationship between high coverage and high soil TN content in the group of forest and grassland.Wetland held the highest TN content among the unused land.

Soil TN of Different Soil Types and Land-Use Types
Soil TN contents at a depth of 0-20 cm under different soil types and land-use patterns are shown in Table S2 and Figure 8.The soil TN varied considerably across the soil groups and ranked in the following order: Histosols > Phaeozems > Chernozems > Luvisols > Cryosols > Kastanozems > Acrisols > Anthrosols > Cambisols > Calcisols > Solonchaks > Arenosols > Fluvisols.The TN content also showed dramatic differences between land-use types.At the great group level, TN contents were highest in the forest, followed by grassland, cropland and construction land, while lowest in unused land.In terms of the subgroups, there was a direct relationship between high coverage and high soil TN content in the group of forest and grassland.Wetland held the highest TN content among the unused land.The average TN contents for the different land-use types in each climate zone are listed in Figure 9.The TN contents increased as precipitation increased and temperature decreased.In each climatic zone, the TN contents in arable land were noticeably lower than that of forest and grassland.The drier the climate was, the greater were the differences in soil TN contents between arable land and the other two vegetated land-use types.In arid regions, the soil TN contents in arable land are only about half of that in the forest.We converted TN contents of arable land into six grades (Table 3) according to the rules of the National Soil Survey Office [51].The TN contents were most in Grade 3 (32.3%)which between 1.0 and 1.5 g•kg −1 .Nitrogen-lacking soil with TN less than 1.0 g•kg −1 accounted for 41.3%.The average TN contents for the different land-use types in each climate zone are listed in Figure 9.The TN contents increased as precipitation increased and temperature decreased.In each climatic zone, the TN contents in arable land were noticeably lower than that of forest and grassland.The drier the climate was, the greater were the differences in soil TN contents between arable land and the other two vegetated land-use types.In arid regions, the soil TN contents in arable land are only about half of that in the forest.The average TN contents for the different land-use types in each climate zone are listed in Figure 9.The TN contents increased as precipitation increased and temperature decreased.In each climatic zone, the TN contents in arable land were noticeably lower than that of forest and grassland.The drier the climate was, the greater were the differences in soil TN contents between arable land and the other two vegetated land-use types.In arid regions, the soil TN contents in arable land are only about half of that in the forest.We converted TN contents of arable land into six grades (Table 3) according to the rules of the National Soil Survey Office [51].The TN contents were most in Grade 3 (32.3%)which between 1.0 and 1.5 g•kg −1 .Nitrogen-lacking soil with TN less than 1.0 g•kg −1 accounted for 41.3%.We converted TN contents of arable land into six grades (Table 3) according to the rules of the National Soil Survey Office [51].The TN contents were most in Grade 3 (32.3%)which between 1.0 and 1.5 g•kg −1 .Nitrogen-lacking soil with TN less than 1.0 g•kg −1 accounted for 41.3%.

Quality of the Prediction
In general, the prediction ability of RF was better than XGBoost in 50 iterations of bootstrapping.RF model can reduce over-learning and over-fitting [46], it was often the best one among multiple machine learning methods, which have been confirmed by many studies [22,52].However, RF model inherited the insensitivity to outliers from recursive partitioning and tree averaging, while XGBoost trained subsequent models using residuals and was more sensitive to outliers [46,47].Among soil data used in this study, some points with extremely high TN contents were considered outliers in statistics.We confirmed the authenticity of these data and kept the records.For these high values, XGBoost model had better prediction accuracy, as shown in Figure 4. Therefore, WMA ensembled two models, kept the advantages and discarded the inaccurate aspect of each algorithm [25].Consequently, WMA exhibited the best competence for capturing the spatial variation of soil TN and reducing prediction uncertainty as well.The results from many other studies have been similar, with researchers concluding that the results were better from an ensemble of different DSM models than from a single model [22,23].Hence, we suggest that it is often useful to combine the predictions of several models because the large variance can be balanced through averaging.
We compared our data with those from other published studies and found that our combined model performed better than other models applied at the national scale for China.Shangguan et al. [53] developed a 30 arc-second (1 km) resolution soil TN map with the polygon linkage method.Li et al. [54] produced a map of TN in topsoil by combining a multiple regression model and neural networks.The mean relative error for the predicted value was 61.06%, which was higher than that in our study (49.17%).Our models were suitable for producing a 90-m resolution national map of soil TN that was more reliable, provided detailed information about the spatial variation in TN, and followed the specifications of the GlobalSoilMap [55].

Spatial Distribution of Soil TN
The distribution map (Figure 5) revealed the spatial pattern of soil TN and gave us insights into how the environment influenced the TN in soil.As one of soil formation factors, topography plays an important role in TN modeling and can determine the hydrothermal conditions and distribution of soil-forming substances [10].Soil TN increased significantly as the elevation increased, possibly reflecting less human disturbance and better moisture-temperature conditions at high altitude [17,28,29].Precipitation and temperature acted as climate proxies and the most robust predictors of TN in our study as shown in Figures 3 and 9.They affected the spatial distribution of TN via changes in the soil moisture content and soil temperature [56].The activity and species of the bacteria that decompose and convert organic matter is limited in low temperature [19,57].The NDVI and NPP are also good predictors of the soil TN as shown in Figures 3 and 8, as the vegetation productivity and biomass are correlated with the amount of litter that is returned to the soil [20].Generally, the TN contents increased with increases in elevation, precipitation and vegetation cover density, decreases in temperature.Soil TN gradually accumulated and eventually reached high levels via slow decomposition, accumulation of biomass, and weak leaching.Thus, we could expect the geographic trend shown in Figure 6, with higher soil TN in Cryosols in the alpine meadow area of the eastern Tibet plateau, and in Phaeozems, Chernozems and Luvisols in the northern forest zone of the Greater Khingan Range, while the lower TN mainly in the desert area.These distribution patterns are consistent with those presented on the national TN maps produced by Li et al. [54] and Shangguan et al. [53].

Uncertainty in Soil TN Prediction
The amount of importance in assessing the uncertainty is almost as great as that in making a prediction map [58], uncertainty map helped researchers to identify the source of the uncertainty and propose solutions [59,60].To obtain a very reliable prediction, we can control the uncertainty of soil TN map from two aspects.
First, we can deal with the uncertainty introduced by a limited number of covariates by adding more variables that are either more relevant or more precise into the model's dataset.In this study, the elevation drops dramatically in the mountainous areas; in these places, the uncertainty is large, reflecting the high heterogeneity in the adjacent geographical factors, which is further translated into inaccurate descriptions of variables when mapped at coarse resolutions.When downscaling, these uncertainties will be further propagated to the soil TN predictions.Moreover, the drainage area of the Three-River Source region is 237,957 km 2 , or 65.9% of the total area [61].The soils here are susceptible to erosion by runoff [62], so considerable amounts of soil TN are probably redistributed [63], thereby weakening the relationship between soil TN and covariates.It would therefore be useful to add quantitative estimates of either soil erosion or the soil loss potential to improve the accuracy of the TN estimates [64].
Second, where the uncertainty originates from insufficient soil profiles, the density of the sampling network could be increased and more samples could be collected.With an adequate dataset, we can capture spatial variance explained by environmental covariates more precisely.When establishing a new database, areas with low certainty should be emphasized and should be a priority for further database investment [65,66].

Effect of Land Use on Soil TN Contents
Figures 8 and 9 demonstrated that lower TN contents in cultivated soils than in soils under natural vegetation cover, similar results were recorded by previous studies [20,67].Human activities play a dominant role in controlling the TN contents in a specific climate region [11].Tillage and other agricultural strategies will destroy the physical protective layer and promote soil respiration, thus accelerate the decomposition of organic matter [68,69].Zhao et al. [70] have certified that after 50 years of cultivation, the TN contents in 0-20 cm decreased by 67%-68%.Consequently, soil TN dropped dramatically when land use shifted from forest or grassland to arable land [67,[71][72][73].
The soil TN levels were low across large-scale cultivated areas, as shown in Table 3, which would have a serious constraint on agricultural yields.Many farmers chose to expand acreage to reach a high total value of crop production, thus China's landscape has undergone large changes since the 1980s, with great losses of woodland and grassland and corresponding gains in cropland [74].Blind reclamation resulted in soil nutrient loss, soil erosion, environmental deterioration and a vicious circle of poverty.A turnaround appeared from 1999 onwards with a rebound in forest and grassland, when the Returning Farmland to Forest Program (RFFP) was introduced.The RFFP compensated farmers for giving marginal cropland back to forest or grassland in order to increase forest cover, improve soil quality and alleviate poverty [75].Soil TN concentration increased after the establishment of the woodland or grassland on abandoned cropland [70,76].However, compared with rapid soil degradation caused by injudicious land use and management, the recovery of soil fertility would take a long period of time once the fragile environment was destroyed [70].Incontrovertibly, there is an important need to make policy decisions about land-use planning that take account of the local conditions rather than applying experience mechanically.These N-rich soils should be used for activities that are suitable for their soil formation characteristics and environmental histories such as forestry and animal husbandry, and these N-lacking soils in arable land can be improved by applying N fertilizer rationally and returning straw [77].

Limitations and Perspectives
In 1980, the annual consumption of nitrogen fertilizer in China had exceeded 9 million tons, and it increased dramatically to over 22 million tons in 2017 [78].The environmental cost of immoderate nitrogen fertilizer has become a dilemma that causes environmental degradation in China including soil acidification, air pollution, water pollution and ecosystem diversity decline [71,79,80].Soil TN content should be monitored so as to know how this has changed since the 1980s under human activities, especially for the nitrogen application.The accuracy of our TN map is limited by the size of the soil database which cannot completely capture the large variation in such a large country.Another limitation of map accuracy results from the relatively low quality of environmental covariates corresponding to the soil sampling year.However, the modeling and predicting procedure proposed in this study can be easily transferred to other time periods under GlobalSoilMap specifications and further applied to the monitoring of soil TN dynamics.
Several issues should be carefully addressed in subsequent works.First, the uncertainties from inaccuracy environmental covariates and inadequate soil profiles should be handled by integrating higher-quality data sources and supplementing samples as mentioned in Section 4.3.Second, some soil properties which correlate with soil TN (e.g., soil texture, pH, soil organic carbon) should be incorporated as environmental covariates to further improve map accuracy.Third, it is worth considering total nitrogen simulation in three dimensions.
Nevertheless, all of the data for mapping soil TN in this study are the most detailed dataset we can acquire for China.Our soil TN map can fill the gap in national TN inventorying at a high spatial resolution and thereby function as a baseline to monitor how the TN content and spatial distribution change with human-induced changes in land use, the development of industry and agriculture, population growth, land management, and climate change.Furthermore, governmental departments that make decisions can formulate corresponding management strategies to improve nitrogen sequestration and mitigate global warming, based on the soil-monitoring network.

Conclusions
Our study is the first contribution of mapping Chinese soil total nitrogen at 90 m resolution.We constructed the prediction model by combining a random forest model with an extreme gradient boosting model.The conclusions of this study are listed below: 1.
Using the weighted average of these two models, a reasonable result was obtained, with the lowest RMSE (1.15 g•kg −1 ) and the highest R 2 (0.41) compared with individual models, that explained 41% of the spatial discrepancy in the soil TN contents and reduced the prediction uncertainty as well.

2.
The TN map showed high spatial heterogeneity, with the spatial variation influenced by variables related to climate, relief and organisms.The spatial trends were similar to previous TN maps in coarser resolution, with high TN in the eastern Tibetan Plateau and north-eastern China, and low TN in the desert area.

3.
The uncertainty map can help policymakers and stakeholders to understand the reliability of the map produced in our study.It should be noted that the uncertainties can be reduced by using more covariates or supplementing the number of soil profiles in the future.
As the most quantitative and precise estimate of the soil TN contents at the national scale, our TN map can serve as a baseline to monitor soil TN dynamics, formulate land-use strategies and support biogeochemical simulations.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/12/1/85/s1: Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 18 where i f is the f(x) mean value of layer i and ( can determine the formula for f(x) by minimizing as follows:

Figure 2 .
Figure 2. Raincloud plot of soil TN content in different (a) land use types, (b) vegetation types and (c) soil types.The "water" group is omitted because it is empty.

Figure 2 .
Figure 2. Raincloud plot of soil TN content in different (a) land use types, (b) vegetation types and (c) soil types.The "water" group is omitted because it is empty.
and upper bounds of CIs.

Figure 3 .
Figure 3. Relative importance of covariates for predicting TN.

Figure 3 .
Figure 3. Relative importance of covariates for predicting TN.

Figure 4 .
Figure 4.The soil TN plot between the measured and averaged predicted values of validation data (log-transformed).For each data point, three prediction values from RF, XGBoost and WMA were linked with the grey line.(b,c) showed part of the detailed information of (a).

Figure 4 .
Figure 4.The soil TN plot between the measured and averaged predicted values of validation data (log-transformed).For each data point, three prediction values from RF, XGBoost and WMA were linked with the grey line.(b,c) showed part of the detailed information of (a).

18 Figure 5 .
Figure 5. Plot of the uncertainty values of 402 validation data points.

Figure 5 .
Figure 5. Plot of the uncertainty values of 402 validation data points.

18 Figure 5 .
Figure 5. Plot of the uncertainty values of 402 validation data points.

Figure 7 .
Figure 7. Uncertainty map of soil total nitrogen.

Figure 7 .
Figure 7. Uncertainty map of soil total nitrogen.

Figure 8 .
Figure 8. Mean soil TN content (g•kg −1 ) in different land-use types in China.

Figure 9 .
Figure 9.The average TN content under different land use types in the climate zones of China.The "others" group contains construction land and other unused land.

Figure 8 .
Figure 8. Mean soil TN content (g•kg −1 ) in different land-use types in China.

18 Figure 8 .
Figure 8. Mean soil TN content (g•kg −1 ) in different land-use types in China.

Figure 9 .
Figure 9.The average TN content under different land use types in the climate zones of China.The "others" group contains construction land and other unused land.

Figure 9 .
Figure 9.The average TN content under different land use types in the climate zones of China.The "others" group contains construction land and other unused land.

Figure S1 :
Topography and climate zones of China.
Figure S2: Flow chart showing the process of spatial modeling.

Figure S3 :
Weights assigned to RF and XGBoost models in 50 iterations.

Table 1 .
Environmental covariates used in modeling soil TN.

Table 2 .
Model diagnostics for independent validation.

Table 2 .
Model diagnostics for independent validation. R

Table S1 :
One-Way ANOVA Output.TableS2: Mean soil TN content (g•kg −1 ) in different soil types in China.