- Email: [email protected]
Application of a Digital Soil Mapping Method in Producing Soil Orders on Mountain Areas of Hong Kong Based on Legacy Soil Data
Pedosphere 21(3): 339–350, 2011 ISSN 1002-0160/CN 32-1315/P c 2011 Soil Science Society of China Published by Elsevier B.V. and Science Press
Application of a Digital Soil Mapping Method in Producing Soil Orders on Mountain Areas of Hong Kong Based on Legacy Soil Data∗1 SUN Xiao-Lin1,2,3 , ZHAO Yu-Guo2,3 , ZHANG Gan-Lin2,3 , WU Sheng-Chun1,3 , MAN Yu-Bon1 and WONG Ming-Hung1,3,∗2 1
Croucher Institute for Environmental Sciences, and Department of Biology, Hong Kong Baptist University, Hong Kong Special Administrative Region (China) 2 State Key Laboratory of Soil and Sustainable Agriculture, Institute of Soil Science, Chinese Academy of Sciences, Nanjing 210008 (China) 3 Joint Open Laboratory of Soil and Environment, Institute of Soil Science, Chinese Academy of Sciences and Hong Kong Baptist University (China) (Received September 30, 2010; revised March 22, 2011)
ABSTRACT Based on legacy soil data from a soil survey conducted recently in the traditional manner in Hong Kong of China, a digital soil mapping method was applied to produce soil order information for mountain areas of Hong Kong. Two modeling methods (decision tree analysis and linear discriminant analysis) were used, and their applications were compared. Much more effort was put on selecting soil covariates for modeling. First, analysis of variance (ANOVA) was used to test the variance of terrain attributes between soil orders. Then, a stepwise procedure was used to select soil covariates for linear discriminant analysis, and a backward removing procedure was developed to select soil covariates for tree modeling. At the same time, ANOVA results, as well as our knowledge and experience on soil mapping, were also taken into account for selecting soil covariates for tree modeling. Two linear discriminant models and four tree models were established finally, and their prediction performances were validated using a multiple jackknifing approach. Results showed that the discriminant model built on ANOVA results performed best, followed by the discriminant model built by stepwise, the tree model built by the backward removing procedure, the tree model built according to knowledge and experience on soil mapping, and the tree model built automatically. The results highlighted the importance of selecting soil covariates in modeling for soil mapping, and suggested the usefulness of methods used in this study for selecting soil covariates. The best discriminant model was finally selected to map soil orders for this area, and validation results showed that thus produced soil order map had a high accuracy. Key Words:
decision tree analysis, linear discriminant analysis, soil covariate selection
Citation: Sun, X. L., Zhao, Y. G., Zhang, G. L., Wu, S. C., Man, Y. B. and Wong, M. H. 2011. Application of a digital soil mapping method in producing soil orders on mountain areas of Hong Kong based on legacy soil data. Pedosphere. 21(3): 339–350.
INTRODUCTION Soil information is urgently needed globally for applications in agriculture (Bishop and McBratney, 2001), ecological (Dobos et al., 2001; Scull et al., 2005) and hydrological process modeling (Zhu et al., 1997). This is particularly true for Hong Kong, with ∗1
rapid industrialization and intensive urbanization happening during the past decades. The local agriculture has shrunk to a negligible size. The fast growing population of 7 million people residing on the 1 104 km2 territory is always posing huge pressure on the local environment and ecology (Li et al., 2001; Warren-Rhodes and Koenig, 2001). Due to the inherent heterogeneity
Supported by the Public Policy Research of the Research Grants Council of Hong Kong, China (No. 2002-PPR-3), the Knowledge Innovation Program of the Chinese Academy of Sciences (No. KZCX2-YW-409), the National Natural Science Foundation of China (Nos. 40625001 and 40771092), and the Mini-AOE (Area of Excellence) Fund from the Hong Kong Baptist University, China (No. RC/AOE/08-09/01). ∗2 Corresponding author. E-mail: [email protected].
340
of local soil, it is difficult to predict slope stability which is threatened by highly erratic nature of local hydrological influences given the local rugged terrain and intense seasonal rainfall (Au, 1998). The available soil information of Hong Kong is far from sufficient to meet various needs. The five small pieces of 1:25 000 soil maps which were drawn a half century ago (Grant, 1962) only focused on the local agricultural areas, about 27% of the whole territory. Several other soil maps covering the whole territory were made from a soil survey conducted recently in the traditional manner using the soil surveyors’ mental model (Luo et al., 2007). However, digital soil mapping has been considered more promising for future soil mapping (McBratney et al., 2003; Grinand et al., 2008; Lagacherie, 2008), and it is now being widely applied to update traditional soil maps based on legacy soil data (Kempen et al., 2009; Stoorvogel et al., 2009). Digital soil mapping is based on soil genesis, and therefore soil information is also inferred from soil covariates on the five soil forming factors, i.e., climate, parent material, relief, organisms and time (McBratney et al., 2003). According to Luo et al. (2007), the local soils are Ferrosols, Cambosols, Ferralosols, Argosols, Gleyosols, Primosols and Anthrosols based on the Chinese Soil Taxonomy (CRGCST, 2001). These soil orders are, respectively, similar to Acrisols, Cambisols, Ferralsols, Luvisols, Gleysols, Leptisols and Anthrosols of the World Reference Base for Soil Resources (IUSS Working Group WRB, 2006). Primosols and Anthrosols usually occur on local flat areas where soil formation is mainly due to human activities, and Gleyosols exclusively occurs on swamps and wetlands of the northern border. Therefore, these three soil orders can be easily identified according to land use. The remaining soil orders are mainly distributed on local mountain areas and their distributions are quite complex being influenced by parent material, relief and vegetation. Thus, this study is aimed at mapping these soil orders. According to Grant (1962), spatial distribution of soils in Hong Kong is mainly controlled by parent material and relief whereas vegetation is not a good indicator for soil change. Luo et al. (2007) indicated that distributions of the targeted soil orders are quite different by parent material and relief, but are very similar by vegetation. Hence, local parent material and relief information can be used as soil covariates to map the soil orders on local mountain areas. In digital soil mapping, prediction is implemented
X. L. SUN et al.
by applying a model describing the relationship between soil and soil covariates. Modeling the relationship is the key step which deserves special attention in digital soil mapping (Kempen et al., 2009). Various modeling methods have been used in digital soil mapping (McBratney et al., 2003). Among them, decision tree analysis is frequently used for predicting soil classes (Henderson et al., 2005; Scull et al., 2005; Grinand et al., 2008). Compared with other methods, this method has advantages of being non-parametric, quantitative and categorical variables compatible, highly efficient, etc. (Grinand et al., 2008). Linear discriminant analysis is also often used for predicting soil classes, and its performance has been considered as easy, efficient and successful (Thomas et al., 1999; Dobos et al., 2000). Although the two methods can choose variables automatically, it is crucial to provide suitable variables because a model for soil mapping should not be only statistically sound but also pedologically plausible (Kempen et al., 2009). Therefore, in this study of digital soil mapping, selecting soil covariates for modeling was emphasized. The major objective of the present study was to refine a soil order map for mountain areas of Hong Kong based on legacy soil data from the recent soil survey of Luo et al. (2007). Two modeling methods, decision tree and linear discriminant analysis, were employed, and different combinations of soil covariates from local dominant soil forming factors, i.e., parent material and relief, were tried. Application of the modeling methods and combinations of soil covariates on mapping were compared. MATERIALS AND METHODS Study area Hong Kong (22◦ 8 –22◦ 36 N, 113◦ 50 –114◦ 23 E) is located in the southeastern tip of China with a total land area of 1 104 km2 (Fig. 1). The climate is subtropical and temperate. The annual air temperature is about 23.0 ◦ C, and the mean annual rainfall is 2 214 mm. About 75% of the land is mountainous. The mountains mainly trend from northeast to southwest. The highest mountain is Tai Mo Shan with an elevation of 957 m, followed by Lan Tau Island with an elevation of 934 m (Fig. 1). Slopes on the mountains are usually very rugged. On foot slopes of the mountains, vegetation is mostly comprised of old Banyan trees (Ficus
DIGITAL SOIL ORDER MAPPING OF HONG KONG
Fig. 1
341
Study area and sample locations in Hong Kong, China.
retusa). On mountainsides not subjected to erosion, vegetation is generally grasses. On higher slopes, vegetation changes to grasses with patches of melastomaceous and ericaceous undershrubs. On the top of high mountains, such as Tai Mo Shan and Lan Tau Island (Fig. 1), the vegetation is mainly turf, largely consisting of Ischaemum spp., containing ground orchids, balsams and mountain Compositae (Grant, 1962). The remaining 25% of the land comprises many flat patches. These areas are separately distributed along the coastlines, in alluvial lowlands, in valleys and in basins. Most of the flat areas are very intensively urbanized, and the remaining only a few farmlands, less than 5% of the whole land area, are adjacent to and separated by the urban area in the New Territories (Fig. 1). On the mid-northern boundary of the New Territories, there are a lot of wetlands and swamps (mainly mangroves) (Fig. 1). Fig. 2 shows the local parent materials based on the 1:200 000 Hong Kong geologic map (The Hong Kong Geological Survey, 1999). According to this map, parental materials are tuff (48%), granite (26%), alluvium/colluvium (11%), sedimentaries (6%), lava (3%), tuff and lava (1%), and fill (5%). Tuff dominates the mountain areas, especially high mountains, followed by granite on relatively lower mountains. Sedimentaries are mixture of sandstone, siltstone, mudstone, graphitic and conglomerate, occurring on lower places. Alluvium/colluvium occurs on flat areas (Fig. 2), some of which are covered by patches of urban areas. Lava occurs on the eastern part of the Lan Tau Island, as well as mixture of tuff and lava. Fill means the material added by human when reclaiming land from the
sea, and it mainly occurs in urban areas. Soil sampling, classification and mapping In the soil survey of Luo et al. (2007), soil sampling was conducted in a traditional manner. The dominant soil forming factors on mountain areas, parent material and relief (Fig. 2), were mainly used to guide soil sampling. After studying the information of parent material and relief, soil surveyors empirically divided Hong Kong soil into 11 kinds of soil association, and several representative soil profiles for each kind of soil association were examined. Totally, 51 soil profiles were examined, and 41 of them were on mountain areas (Fig. 1). The morphological characteristics, such as soil horizon depth, Munsell color, structure, etc., of all examined profiles were recorded. A number of physical and chemical properties of the profiles related to soil classification were analyzed in order to classify the profiles. Details about these analyses refer to Luo et al. (2007). Based on measured soil properties, the profiles were classified as Ferrosols, Ferralosols, Argosols and Cambosols according to the Chinese Soil Taxonomy (CRGCST, 2001), as indicated earlier. Ferrosols are characterized as pH < 5.0, base saturation less than 50%, having a lower activity-ferric horizon and an allitic property. Ferralosols are characterized as acid, alumino-silicate minerals accumulated, having a lower activity-ferric horizon and high clay content. Argosols are characterized as weak acid, illuvial accumulations of clay in B horizon, a low ratio of SiO2 /Al2 O3 and free iron less than 40% of total iron in soil body. Cambosols are characterized as lower accumulations of clay in B horizon and a lower ratio of SiO2 /Al2 O3 in sub-
342
Fig. 2
X. L. SUN et al.
Information related to dominated soil forming factors: parent material and relief.
soil. The 4 soil orders are similar to Ultisols, Oxisols, Alfisols and Inceptisols of the U.S. Soil Taxonomy. Soil mapping in this survey was based on the hypothesis that similar soil formation conditions will generate similar soils. Information of the two dominant soil forming factors (Fig. 2) was analyzed, and Hong Kong area was divided into a number of landformparent material groups. By assigning representative soil taxa to each kind of landform-parent material group, a soil group map was established. In the present study, this soil group map was used to produce a soil order map (Fig. 3), by assigning a soil order of a soil group to the corresponding areas of that soil group. Soil covariates Soil covariates used in this study were also from the two dominant soil forming factors: parent mate-
rial and relief (Fig. 2). Based on the relief information which was a 15 m digital elevation model (DEM) constructed based on the 1:50 000 Hong Kong topographic map using ArcGIS software, many terrain attributes were derived: elevation (Z, m), slope (S, ◦ ), aspect (A, ◦ ), profile curvature [Kp , m (100 m)−1 ], plan curvature [Kc , m (100 m)−1 ], natural logarithm of specific catchment area [ln(As ), m2 m−1 ], stream power index (SPI) and topographic wetness index (TWI). These terrain attributes are commonly used in digital soil mapping (McBratney et al., 2003), and details about them refer to Florinsky et al. (2002). Considering the non-linear variation of aspect as a circular variant, the sine of aspect (sinA) was used instead of aspect itself. Decision tree analysis and linear discriminant analysis Decision tree analysis is commonly used in soil ma-
DIGITAL SOIL ORDER MAPPING OF HONG KONG
Fig. 3
343
Soil order map based on the soil group map from the soil survey of Luo et al. (2007).
pping (Dobos et al., 2001; Scull et al., 2005; Grinand et al., 2008). For example, Scull et al. (2005) used this method to predict soil taxonomic unit from soil covariates. In this study, C5.0 (Quinlan, 2002) was used to implement decision tree analysis. A tree model must be pruned before using it to classify new samples. In C5.0, a tree model was firstly pruned by combining a leaf node or a branch that was predicted to have a relatively high misclassification rate (%, percentage of misclassified objects over all classified objects) with another one. This pruning process was applied to every branch and every leaf node. Following this pruning process, globally pruning continued to prune the tree model by examining performance of the tree model as a whole. In this study, decision tree analysis was implemented based on all of the soil covariates mentioned in the section of soil covariates. Linear discriminant analysis is also a commonly used method for digital soil mapping (McBratney et al., 2003). For example, Thomas et al. (1999) used linear discriminant analysis to predict soil classes based on terrain attributes derived from a DEM. Linear discriminant function is determined by a measure of generalized squared distance based on the pooled covariance matrix (Sinowski and Auerswald, 1999). Scores of these discriminant functions on new cases were used to classify the new cases. In this study, this analysis was implemented by the DISCRIM procedure in SAS system. Since linear discriminant analysis can not handle categorical variables, only the terrain attributes stated in the section of soil covariates were used in this analysis.
Selection of soil covariates for decision tree analysis and linear discriminant analysis Selection of soil covariates used for digital soil mapping is an “acute” problem (Lark et al., 2007). Though a lot of data mining methods can automatically choose the best variables to build a model, they do this just based on the given set of variables by modelers. For different given sets of soil covariates, different modeling results will be generated. Therefore, in this study, different sets of soil covariates were tried. Before decision tree analysis and linear discriminant analysis, variances of those calculated terrain attributes between soil orders were tested by analysis of variance (ANOVA) in SAS system. Results of this analysis could indicate the importance of a terrain attribute in differentiating soil orders, which would be suggestive for selecting variables in linear discriminant analysis and decision tree analysis. Linear discriminant analysis was conducted through the stepwise DISCRIM procedure. Selecting variables in this procedure was based on statistics of Wilks’s Lambda and average squared canonical correlation. Meanwhile, results from ANOVA were also considered in building discriminants. For the 9 soil covariates mentioned above, there would be 512 combinations of them, which obviously would generate different mapping results. Assuming that the more closely a soil covariate is related to local soil formation, the more it will influence the results of decision tree analysis and vice versa, a backward removing procedure was developed to screen the most
344
and least influential soil covariates after trying all of the 9 soil covariates in decision tree analysis to find: (1) the most influential soil covariates by removing that attribute with the least usage in decision tree analysis, and (2) the least influential soil covariates by removing that attribute with the greatest usage in decision tree analysis. Subsequently, decision tree analysis was conducted again based on the remaining attributes. These processes were repeated until the performance of decision tree analysis became very poor (indicated by misclassification rate). Based on the results of backward removing variables in decision tree analysis, as well as the results of ANOVA, the most influential soil forming factors, i.e., the most closely related factors to soil formation, were identified and were selected as the best combinations of soil covariates to build a tree model. For comparison, some other sets of soil covariates were selected to build tree models, as well as one selected according to the knowledge on local soil genesis and our experience on soil mapping. Validation and the final soil order map Considering the very limited sample size in the legacy soil data, this study used the approach proposed by Bishop and McBratney (2001), i.e., multiple jackknifing, to validate prediction results. This approach involved random separating the whole sampling dataset in a ratio into two subsets, one for prediction and one for validation. The random separation was conducted a number of times. Each time, a model was constructed based on the prediction subset and then validated based on the validation subset. Mean accuracy from this cross-validation over a number of times was used to identify the best model. In this study, the sampling dataset was separated in the ratio of about 90% to 10%, i.e., 35 samples for prediction and 6 samples for validation. Accuracy of an established model was assessed in terms of misclassification rate. This validation procedure was randomly repeated 25 times. The best model was identified in terms of the mean misclassification rate over the 25 times. The best model of those established during the above validation was selected to make the final soil order map. The criterion for selecting the best model was prediction accuracy in terms of misclassification rate and kappa statistic (Sim and Wright, 2005). Kappa statistic is an index which compares the agreement against that which might be expected by chance. It can be thought of as the chance-corrected proportional agreement (K), and possible values range from 1 (per-
X. L. SUN et al.
fect agreement) to 0 (no agreement above that expected by chance) to −1 (complete disagreement). It was computed as: K=
Po − Pc 1 − Pc
(1)
Po =
np × 100% no
(2)
C nip × nio
Pc =
i=1
no no
(3)
where Po is the observed agreement, Pc is the chance agreement, C is the number of categories of the tested samples, no is the total number of tested samples, np is the total number of correctly predicted tested samples, nip is the total number of tested samples which are predicted into the ith category, nio is the total number of tested samples which actually are the ith category. The accuracy of the final map was also assessed in terms of producer accuracy and user accuracy (Scull et al., 2005). Producer accuracy is a measure of how well the test data are classified. User accuracy is a measure of how likely a test sample classified into a given category actually belongs to that category. The flowchart of making the final soil order map is shown in Fig. 4.
Fig. 4 The flowchart for mapping soil order on mountain areas of Hong Kong, China.
DIGITAL SOIL ORDER MAPPING OF HONG KONG
RESULTS AND DISCUSSION Variances of terrain attributes between soil orders Table I showed that only elevation had significant differences (P < 0.05) between soil orders, followed by ln(As ) with difference between soil orders approaching the significant level of 0.05 (P = 0.06). The variance of slope between soil orders was also wider than the others. The significant difference in elevation between soil orders indicated that soil orders over Hong Kong had an obvious vertical zonation. Average elevations of the four soil orders were 294, 216, 137 and 61 m, respectively, for Ferrosols, Argosols, Ferralosols and Cambosols. The relatively significant difference of ln(As ) (P = 0.06) indicated the importance of water flow in soil processes, as eluviation induced by water flow plays a major role in the formation of Ferrosols, Ferralosols and Argosols (CRGCST, 2001). The influence of slope on soil processes through redistributing materials and energy over space was slightly strong, as evident by its relatively more remarkable differences between soil orders than other terrain attributes, though it was not statistically significant (P = 0.1) (Table I). Other terrain attributes seemed to be less related to the distriTABLE I Analysis of variance of terrain attributesa) between soil orders on mountain areas of Hong Kong, China Item
Z
S
sinA Kp
Kc
ln(As ) SPI TWI
F value 3.33 2.17 0.12 1.02 0.67 2.62 P > F b) 0.03* 0.10 0.94 0.39 0.57 0.06
1.00 1.71 0.40 0.18
*Significant at 0.05 level. a) Z = elevation, S = slope, A = aspect, Kp = profile curvature, Kc = plan curvature, As = specific catchment area, SPI = stream power index, TWI = topographic wetness index. b) The confidence of F value.
345
bution of the four soil orders, according to the results shown in Table I. The results of ANOVA showed that the relationship between these soil orders and a single terrain attribute was not very significant. This may be due to the fact that formation of the soil orders was not governed by a single terrain attribute, but by a number of factors, including not only terrain attributes, but also other soil forming factors and their interactions (CRGCST, 2001). However, the results of ANOVA would be useful for identifying key soil covariates in the following linear discriminant analysis and decision tree analysis. Soil covariates in linear discriminant analysis Table II showed that five terrain attributes were statistically significant in constructing linear discriminants, according to the statistics of Wilks’ lambda and average squared canonical correlation. Elevation contributed the most to the discriminations, followed by ln(As ), slope, TWI and profile curvature. This sequence was the same as that presented in Table I because stepwise DISCRIM and ANOVA shared the same way to assess contributions of the terrain attributes, i.e., the way of analyzing variances of the terrain attributes among soil orders. However, ANOVA only considered one single terrain attribute at a time, while stepwise DISCRIM considered several at a time. Thus, more terrain attributes appeared as significant variables in linear discriminant analysis compared with ANOVA. Two sets of soil covariates were selected to model discriminant functions. The first set (DI) was selected according to the results of ANOVA that showed the most significant influences of elevation, ln(As ) and slope on soil formation. Therefore, the first set (DI) comprised these three terrain attributes. The second
TABLE II Results of linear discriminant analysis for the five terrain attributesa) through the stepwise DISCRIM procedure Step
Terrain attribute
Wilks’ lambda
P < Wilks’ lambda
Average squared canonical correlation (ASCC)
P < ASCC
1 2 3 4 5
Z ln(As ) S TWI Kp
0.78 0.65 0.55 0.49 0.47
0.03* 0.01* 0.01* 0.01* 0.03*
0.07 0.13 0.17 0.21 0.22
0.03* 0.01* 0.009** 0.01* 0.02*
*, **Significant at P = 0.05 and P = 0.01 levels, respectively. a) Z = elevation, As = specific catchment area, S = slope, TWI = topographic wetness index, Kp = profile curvature.
346
X. L. SUN et al.
results suggested that parent material, elevation, profile curvature and ln(As ) were the most influential soil covariates on decision tree modeling of the local soillandscape relationships. M1 to M8 listed in the middle part of Table III shows the results of removing the most important soil covariates. Number of leaf nodes and misclassification rates varied dramatically when parent material, elevation, slope, profile curvature, ln(As ), plan curvature and TWI were removed one by one. Comparatively, these indexes varied more dramatically in the middle part of Table III than in the upper part, indicating a greater change in the performances and structures of modeled trees. It reflected that the more important a variable was, the more it influenced tree modeling. As a consequence, plan curvature and SPI were determined to be the least influencing variables in decision tree modeling of the local soil-landscape relationships. The above results showed that slope was removed very early during the process of removing both the most and least important variables. Thus, it was not easy to determine whether slope was important or not. However, results of ANOVA in Table I showed that
set (DII) was selected according to the results of linear discriminant analysis in Table II, comprising elevation, ln(As ), slope, TWI and profile curvature. Two discriminant functions were modeled based on these two sets of soil covariates. On the whole training dataset, their misclassification rates were 46% and 44%, respectively. Soil covariates in decision tree analysis Table III shows the results from implementing the procedure of backward removing variables in decision tree analysis. L1 to L8 listed in the upper part of Table III shows the results of removing the least important soil covariates. In the first five trees, SPI, slope, TWI and plan curvature were removed one by one. All of the misclassification rates of these trees were 0%, and changes in the number of leaf nodes were few. The usages of most soil covariates also changed slightly. In the last three trees, sine of aspect and ln(As ) were removed one by one. Both misclassification rates and usages of the remaining soil covariates increased slightly, with the number of leaf nodes changed a little more. These TABLE III
Decision tree modeling results of removing the least and most important soil covariatesa) Setb)
Parent material
Z
S
sinA
Kp
Kc
TWI
SPI
ln(As )
Leaf nodes
Misclassification rate %
L1 L2 L3 L4 L5 L6 L7 L8
100 100 100 100 100 100 100 100
63 63 63 63 63 63 83 83
M1 M2 M3 M4 M5 M6 M7 M8
100 -
LM KE
100 100
a)
10 10 -
17 17 17 27 27 27 -
63 63 63 71 71 71 73 73
63 100 -
10 27 100 -
17 46 32 78 61 5 22 7
63 68 44 100 -
76 63
20 10
-
63 68
Least important 12 12 0 12 12 12 12 24 12 Most important 12 12 0 20 63 7 39 0 15 95 73 0 83 63 0 100 100 76 100 100 100 76 -
39
59
-
39 39 39 39 39 39 49 -
17 17 18 18 18 18 19 16
0.0 0.0 0.0 0.0 0.0 2.4 2.4 9.8
39 68 71 61 100 -
17 14 14 16 16 10 9 7
0.0 7.3 12.2 7.3 2.4 24.4 31.7 36.6
49
19 19
0.0 0.0
-
Z = elevation, S = slope, A = aspect, Kp = profile curvature, Kc = plan curvature, TWI = topographic wetness index, SPI = stream power index, As = specific catchment area. b) LM: removing the least and most important soil covariates together; KE: knowledge and experience in soil mapping.
DIGITAL SOIL ORDER MAPPING OF HONG KONG
slope was relatively closely related to distributions of soil orders (P = 0.1), and in fact slope was very frequently used for soil mapping (McBratney et al., 2003). Therefore, slope was finally used in the subsequent prediction of soil orders. Contrary to slope, sinA was removed very late during the both processes, but results of ANOVA showed that this soil covariate was scarcely related to distributions of soil orders (P = 0.9). Therefore, sinA was not used in the subsequent prediction of soil orders. Based on the above analysis, the first set of soil covariates for decision tree modeling was designed to include all of the most influential soil covariates: parent material, elevation, slope, profile curvature and ln(As ), i.e., LM listed in Table III. With the objective of generating a good quality soil order map for this area, other three sets of soil covariates, i.e., L1 , L5 and KE in Table III, were also selected to construct decision tree models. L1 and L5 were chosen because they were two extremes in the process of backward removing the most important variables. L1 contained all soil covariates, while L5 contained the fewest but keeping the lowest misclassification rate. KE was established according to the common knowledge about local soil genesis and our experience in soil mapping. As mentioned earlier, parent material and relief played dominant roles in the local soil development. McSweeney et al. (1994) contended that on the watershed scale, soil formation and development were mainly influenced by five environmental variables: elevation, slope, profile curvature, plan curvature and TWI. However, ln(As ), depicting the contribution area of water flowing across a point, would be more appropriate in reflecting the influence of water on soil processes than TWI, which depicts stationary water content in soil. This is because water flow passing through soil bodies or on soil surfaces affects soil more directly and more strongly than stationary water in soil bodies (Florinsky et al., 2002; Thompson et al., 2006). Particularly formations of Ferrolsols, Ferralosols and Argosols are more related to water movement than water content (CRGCST, 2001). Furthermore, results of ANOVA shown in Table I also indicated that ln(As ) was more related to distributions of soil orders than TWI. As a consequence, KE was designed to comprise parent material, elevation, slope, profile curvature, plan curvature and ln(As ). Based on the four designed sets of soil covariates, four decision trees were modeled respectively. Performances of these four tree models are presented in Table III. Misclassification rates of them were all 0%, and
347
their numbers of leaf nodes were close, i.e., 17, 18, 19 and 19. Mapping soil orders based on selected soil covariates Results from validating the above selected two discriminant functions and four decision tree models using the multiple jackknifing validation approach are summarized in Table IV. Generally, the two discriminant functions performed better than the four tree models, although they did not explore the usefulness of parent material for predicting soil orders. Between the two discriminant functions, model DI performed much better than model DII, reflecting that a more straightforward but good enough discriminant function was much more useful. Among the four tree models, model LM which was selected through backward removing procedure performed best, followed by model KE which was selected according to knowledge and experience on soil mapping, model L5 and model L1 . Results from validating the four tree models reflected that: 1) it was important to select soil covariates for tree modeling, even though decision tree analysis can automatically choose best variables for modeling; 2) backward removing procedure improved tree modeling, as well as knowledge and experience on soil mapping; 3) a more straightforward but good enough tree model was much more useful, like that for discriminant functions. TABLE IV Summaries on misclassification rate of selected two discriminant functions (DI and DII) and four decision tree models (L1 , L5 , LM and KE) using the multiple jackknifing validation approach Model Mean Minimum Maximum Standard deviation % DI DII L1 L5 LM KE
49 59 71 69 67 68
33 33 33 33 33 33
83 100 100 100 100 100
14.55 16.75 15.18 14.29 15.64 14.77
According to mean misclassification rates listed in Table IV, model DI with the lowest mean misclassification rate was identified as the best model for producing the final soil order map. Among all the 25 models which were established using this model during jackknifing validation, the one having the lowest misclassification rate (33%) and the highest kappa statistic
348
X. L. SUN et al.
(0.63) was selected to generate the final soil order map (Fig. 5). User accuracy and producer accuracy of this map from jackknifing validation are presented in Table V. The table shows that the final map was very accurate to reflect the distribution of Ferrosols and Ferralosols. Distribution of Argosols was very poorly predicted and distribution of Cambosols was the worst predicted. Weighted by the area percentages of soil orders on the map, average producer accuracy and average user accuracy of this map were 98.7% and 82.6%, respectively, showing that the final soil order map had a high accuracy. The final digital soil order map (Fig. 5) was very different from the traditional one (Fig. 3). Apparently, the digital soil order map was much more detailed in reflecting the distribution of soil orders. Comparison between the two maps showed that Cambosols in the digital map were much less than in the traditional map and were separately distributed among other soil orders.
Ferrolosols were also reduced substantially on the digital map compared with the traditional map. The two soil orders seemed to be replaced by Argosols in Fig. 5, while some Argosols were replaced by Ferralosols. Totally, the digital map reproduced only 18.4% of the traditional one. This big gap between them is due first to the different models used to produce the two maps and second but more crucially, it is most likely due to the limited sample size of the legacy soil data on which this study was based. CONCLUSIONS Linear discriminant models performed better than decision tree models. ANOVA was very useful for selecting soil covariates for building a discriminant. The removing procedure developed in this study for selecting soil covariates for tree modeling was useful to improve the performance of tree modeling. This was also true
Fig. 5 The final soil order map of Hong Kong predicted by model DI, which included three soil covariates, i.e., elevation, natural logarithm of specific catchment area and slope, for modeling discriminant functions. TABLE V Error matrix of final soil order map Predicted soil order
Argosol Cambosol Ferralosol Ferrosol Producer accuracy (%)
Observed soil order
User accuracy
Argosols
Cambosols
Ferralosols
Ferrosols
1 0 0 0 100
2 0 0 0 0
0 0 1 0 100
0 0 0 2 100
% 33 0 100 100
DIGITAL SOIL ORDER MAPPING OF HONG KONG
for knowledge and experience on soil mapping. Clearly, results of this study highlighted the importance of selecting soil covariates in modeling for soil mapping, and suggested the usefulness of the methods used in this study for selecting soil covariates. A soil order map was finally made using the best model established in this study, with a high accuracy in terms of misclassification rate, kappa statistics, producer accuracy and user accuracy. Comparison between the digital soil order map and the traditional one showed that the former just reproduced 18.4% of the later. The big gap between them can be attributed to the very limited sample size of the legacy soil data on which this study was based. Thus, more soil samples in future soil surveys for this area will be necessary in order to improve the final soil order map, as well as to produce more detailed soil maps. ACKNOWLEDGEMENT The authors are grateful to anonymous reviewers for their valuable comments and suggestions, and Dr. A. O. W. Leung at the Croucher Institute for Environmental Sciences, Hong Kong Baptist University, Hong Kong Special Administrative Region, China for editing the language of this manuscript. REFERENCES Au, S. W. C. 1998. Rain-induced slope instability in Hong Kong. Eng. Geol. 51: 1–36. Bishop, T. F. A. and McBratney, A. B. 2001. A comparison of prediction methods for the creation of field-extent soil property maps. Geoderma. 103: 149–160. Cooperative Research Group on Chinese Soil Taxonomy (CRGCST). 2001. Chinese Soil Taxonomy. Science Press, Beijing and New York. Dobos, E., Micheli, E., Baumgardner, M. F., Biehl, L. and Helt, T. 2000. Use of combined digital elevation model and satellite radiometric data for regional soil mapping. Geoderma. 97: 367–391. Dobos, E., Montanarella, L., N`egre, T. and Micheli, E. 2001. A regional scale of soil mapping approach using intergrated AVHRR and DEM data. Int. J. Appl. Earth Observ. Geoinf. 3: 30–42. Florinsky, I. V., Eilers, R. G., Manning, G. R. and Fuller, L. G. 2002. Prediction of soil properties by digital terrain modelling. Environ. Modell. Softw. 17: 295–311. Grant, C. J. 1962. The Soils and Agriculture of Hong Kong. Government Press, Hong Kong. Grinand, C., Arrouays, D., Laroche, B. and Martin, M. P. 2008. Extrapolating regional soil landscapes from an existing soil map: Sampling intensity, validation pro-
349 cedures, and integration of spatial context. Geoderma. 143: 180–190. Henderson, B. L., Bui, E. N., Moran, C. J. and Simon, D. A. P. 2005. Australia-wide predictions of soil properties using decision trees. Geoderma. 124: 383–398. IUSS Working Group WRB. 2006. World reference base for soil resources 2006. World Soil Resources Reports No. 103. FAO, Rome. Kempen, B., Brus, D. J., Heuvelink, G. B. M. and Stoorvogel, J. J. 2009. Updating the 1:50 000 Dutch soil map using legacy soil data: A multinomial logistic regression approach. Geoderma. 151: 311–326. Lagacherie, P. 2008. Digital soil mapping: a state of the art. In Hartemink, A. E., McBratney, A. B. and MendoncaSantos, M. L. (eds.) Digital Soil Mapping with Limited Data. Springer, Dordrecht. pp. 3–14. Lark, R. M., Bishop, T. F. A. and Webster, R. 2007. Using expert knowledge with control of false discovery rate to select regressors for prediction of soil properties. Geoderma. 138: 65–78. Li, X., Poon, C. and Liu, P. S. 2001. Heavy metal contamination of urban soils and street dusts in Hong Kong. Appl. Geochem. 16: 1361–1368. Luo, Y. M., Li, Z. G., Wu, L. H., Wu, S. C., Zhang, G. L., Zhou, S. L., Zhao, Y. G., Zhao, Q. G., Wong, M. H. and Zhang, H. B. 2007. Hong Kong Soils and Environment (in Chinese). Science Press, Beijing. McBratney, A. B., Mendonca-Santos, M. L. and Minasny, B. 2003. On digital soil mapping. Geoderma. 117: 3–52. McSweeney, K., Gessler, P. E., Slater, B. K., Hammer, R. D., Bell, J. and Petersen, G. W. 1994. Towards a new framework for modeling the soil-landscape continuum. In Amundson, R, G., Harden, J. W. and Singer, M. (eds.) Factors of Soil Formation: A Fiftieth Anniversary Retrospective. Soil Science Society of America, Madison, WI. pp. 127–145. Quinlan, J. R. 2002. Data Mining Tools See5 and C5.0. Available online at http://www.rulequest.com/see5info.html (verified on January, 2011). Scull, P., Franklin, J. and Chadwick, O. A. 2005. The application of classification tree analysis to soil type prediction in a desert landscape. Ecol. Model. 181: 1–15. Sim, J. and Wright, C. C. 2005. The kappa statistic in reliability studies: use, interpretation, and sample size requirements. Phys. Ther. 85: 257–268. Sinowski, W. and Auerswald, K. 1999. Using relief parameters in a discriminant analysis to stratify geological areas with different spatial variability of soil properties. Geoderma. 89: 113–128. Stoorvogel, J. J., Kempen, B., Heuvelink, G. B. M. and de Bruin, S. 2009. Implementation and evaluation of existing knowledge for digital soil mapping in Senegal. Geoderma. 149: 161–170. The Hong Kong Geological Survey. 1999. Geological Map
350 of Hong Kong in 1:200 000. 2nd Edition. Survey and Mapping Office, Lands Department, the Government of Hong Kong Special Administrative Region, Hong Kong. Thomas, A. L., King, D., Dambrine, E., Couturie, A. and Roque, J. 1999. Predicting soil classes with parameters derived from relief and geologic materials in a sandstone region of the Vosges mountains (Northeastern France). Geoderma. 90: 291–305. Thompson, J. A., Pena-Yewtukhiw, E. M. and Grove, J. H.
X. L. SUN et al. 2006. Soil-landscape modeling across a physiographic region: topographic patterns and model transportability. Geoderma. 133: 57–70. Warren-Rhodes, K. and Koenig, A. 2001. Ecosystem appropriation by Hong Kong and its implications for sustainable development. Ecol. Econ. 39: 347–359. Zhu, A. X., Band, L., Vertessy, R. and Dutton, B. 1997. Derivation of soil properties using a soil land inference model (SoLIM). Soil Sci. Soc. Am. J. 61: 523–533.
Application of a Digital Soil Mapping Method in Producing Soil Orders on Mountain Areas of Hong Kong Based on Legacy Soil Data∗1 SUN Xiao-Lin1,2,3 , ZHAO Yu-Guo2,3 , ZHANG Gan-Lin2,3 , WU Sheng-Chun1,3 , MAN Yu-Bon1 and WONG Ming-Hung1,3,∗2 1
Croucher Institute for Environmental Sciences, and Department of Biology, Hong Kong Baptist University, Hong Kong Special Administrative Region (China) 2 State Key Laboratory of Soil and Sustainable Agriculture, Institute of Soil Science, Chinese Academy of Sciences, Nanjing 210008 (China) 3 Joint Open Laboratory of Soil and Environment, Institute of Soil Science, Chinese Academy of Sciences and Hong Kong Baptist University (China) (Received September 30, 2010; revised March 22, 2011)
ABSTRACT Based on legacy soil data from a soil survey conducted recently in the traditional manner in Hong Kong of China, a digital soil mapping method was applied to produce soil order information for mountain areas of Hong Kong. Two modeling methods (decision tree analysis and linear discriminant analysis) were used, and their applications were compared. Much more effort was put on selecting soil covariates for modeling. First, analysis of variance (ANOVA) was used to test the variance of terrain attributes between soil orders. Then, a stepwise procedure was used to select soil covariates for linear discriminant analysis, and a backward removing procedure was developed to select soil covariates for tree modeling. At the same time, ANOVA results, as well as our knowledge and experience on soil mapping, were also taken into account for selecting soil covariates for tree modeling. Two linear discriminant models and four tree models were established finally, and their prediction performances were validated using a multiple jackknifing approach. Results showed that the discriminant model built on ANOVA results performed best, followed by the discriminant model built by stepwise, the tree model built by the backward removing procedure, the tree model built according to knowledge and experience on soil mapping, and the tree model built automatically. The results highlighted the importance of selecting soil covariates in modeling for soil mapping, and suggested the usefulness of methods used in this study for selecting soil covariates. The best discriminant model was finally selected to map soil orders for this area, and validation results showed that thus produced soil order map had a high accuracy. Key Words:
decision tree analysis, linear discriminant analysis, soil covariate selection
Citation: Sun, X. L., Zhao, Y. G., Zhang, G. L., Wu, S. C., Man, Y. B. and Wong, M. H. 2011. Application of a digital soil mapping method in producing soil orders on mountain areas of Hong Kong based on legacy soil data. Pedosphere. 21(3): 339–350.
INTRODUCTION Soil information is urgently needed globally for applications in agriculture (Bishop and McBratney, 2001), ecological (Dobos et al., 2001; Scull et al., 2005) and hydrological process modeling (Zhu et al., 1997). This is particularly true for Hong Kong, with ∗1
rapid industrialization and intensive urbanization happening during the past decades. The local agriculture has shrunk to a negligible size. The fast growing population of 7 million people residing on the 1 104 km2 territory is always posing huge pressure on the local environment and ecology (Li et al., 2001; Warren-Rhodes and Koenig, 2001). Due to the inherent heterogeneity
Supported by the Public Policy Research of the Research Grants Council of Hong Kong, China (No. 2002-PPR-3), the Knowledge Innovation Program of the Chinese Academy of Sciences (No. KZCX2-YW-409), the National Natural Science Foundation of China (Nos. 40625001 and 40771092), and the Mini-AOE (Area of Excellence) Fund from the Hong Kong Baptist University, China (No. RC/AOE/08-09/01). ∗2 Corresponding author. E-mail: [email protected].
340
of local soil, it is difficult to predict slope stability which is threatened by highly erratic nature of local hydrological influences given the local rugged terrain and intense seasonal rainfall (Au, 1998). The available soil information of Hong Kong is far from sufficient to meet various needs. The five small pieces of 1:25 000 soil maps which were drawn a half century ago (Grant, 1962) only focused on the local agricultural areas, about 27% of the whole territory. Several other soil maps covering the whole territory were made from a soil survey conducted recently in the traditional manner using the soil surveyors’ mental model (Luo et al., 2007). However, digital soil mapping has been considered more promising for future soil mapping (McBratney et al., 2003; Grinand et al., 2008; Lagacherie, 2008), and it is now being widely applied to update traditional soil maps based on legacy soil data (Kempen et al., 2009; Stoorvogel et al., 2009). Digital soil mapping is based on soil genesis, and therefore soil information is also inferred from soil covariates on the five soil forming factors, i.e., climate, parent material, relief, organisms and time (McBratney et al., 2003). According to Luo et al. (2007), the local soils are Ferrosols, Cambosols, Ferralosols, Argosols, Gleyosols, Primosols and Anthrosols based on the Chinese Soil Taxonomy (CRGCST, 2001). These soil orders are, respectively, similar to Acrisols, Cambisols, Ferralsols, Luvisols, Gleysols, Leptisols and Anthrosols of the World Reference Base for Soil Resources (IUSS Working Group WRB, 2006). Primosols and Anthrosols usually occur on local flat areas where soil formation is mainly due to human activities, and Gleyosols exclusively occurs on swamps and wetlands of the northern border. Therefore, these three soil orders can be easily identified according to land use. The remaining soil orders are mainly distributed on local mountain areas and their distributions are quite complex being influenced by parent material, relief and vegetation. Thus, this study is aimed at mapping these soil orders. According to Grant (1962), spatial distribution of soils in Hong Kong is mainly controlled by parent material and relief whereas vegetation is not a good indicator for soil change. Luo et al. (2007) indicated that distributions of the targeted soil orders are quite different by parent material and relief, but are very similar by vegetation. Hence, local parent material and relief information can be used as soil covariates to map the soil orders on local mountain areas. In digital soil mapping, prediction is implemented
X. L. SUN et al.
by applying a model describing the relationship between soil and soil covariates. Modeling the relationship is the key step which deserves special attention in digital soil mapping (Kempen et al., 2009). Various modeling methods have been used in digital soil mapping (McBratney et al., 2003). Among them, decision tree analysis is frequently used for predicting soil classes (Henderson et al., 2005; Scull et al., 2005; Grinand et al., 2008). Compared with other methods, this method has advantages of being non-parametric, quantitative and categorical variables compatible, highly efficient, etc. (Grinand et al., 2008). Linear discriminant analysis is also often used for predicting soil classes, and its performance has been considered as easy, efficient and successful (Thomas et al., 1999; Dobos et al., 2000). Although the two methods can choose variables automatically, it is crucial to provide suitable variables because a model for soil mapping should not be only statistically sound but also pedologically plausible (Kempen et al., 2009). Therefore, in this study of digital soil mapping, selecting soil covariates for modeling was emphasized. The major objective of the present study was to refine a soil order map for mountain areas of Hong Kong based on legacy soil data from the recent soil survey of Luo et al. (2007). Two modeling methods, decision tree and linear discriminant analysis, were employed, and different combinations of soil covariates from local dominant soil forming factors, i.e., parent material and relief, were tried. Application of the modeling methods and combinations of soil covariates on mapping were compared. MATERIALS AND METHODS Study area Hong Kong (22◦ 8 –22◦ 36 N, 113◦ 50 –114◦ 23 E) is located in the southeastern tip of China with a total land area of 1 104 km2 (Fig. 1). The climate is subtropical and temperate. The annual air temperature is about 23.0 ◦ C, and the mean annual rainfall is 2 214 mm. About 75% of the land is mountainous. The mountains mainly trend from northeast to southwest. The highest mountain is Tai Mo Shan with an elevation of 957 m, followed by Lan Tau Island with an elevation of 934 m (Fig. 1). Slopes on the mountains are usually very rugged. On foot slopes of the mountains, vegetation is mostly comprised of old Banyan trees (Ficus
DIGITAL SOIL ORDER MAPPING OF HONG KONG
Fig. 1
341
Study area and sample locations in Hong Kong, China.
retusa). On mountainsides not subjected to erosion, vegetation is generally grasses. On higher slopes, vegetation changes to grasses with patches of melastomaceous and ericaceous undershrubs. On the top of high mountains, such as Tai Mo Shan and Lan Tau Island (Fig. 1), the vegetation is mainly turf, largely consisting of Ischaemum spp., containing ground orchids, balsams and mountain Compositae (Grant, 1962). The remaining 25% of the land comprises many flat patches. These areas are separately distributed along the coastlines, in alluvial lowlands, in valleys and in basins. Most of the flat areas are very intensively urbanized, and the remaining only a few farmlands, less than 5% of the whole land area, are adjacent to and separated by the urban area in the New Territories (Fig. 1). On the mid-northern boundary of the New Territories, there are a lot of wetlands and swamps (mainly mangroves) (Fig. 1). Fig. 2 shows the local parent materials based on the 1:200 000 Hong Kong geologic map (The Hong Kong Geological Survey, 1999). According to this map, parental materials are tuff (48%), granite (26%), alluvium/colluvium (11%), sedimentaries (6%), lava (3%), tuff and lava (1%), and fill (5%). Tuff dominates the mountain areas, especially high mountains, followed by granite on relatively lower mountains. Sedimentaries are mixture of sandstone, siltstone, mudstone, graphitic and conglomerate, occurring on lower places. Alluvium/colluvium occurs on flat areas (Fig. 2), some of which are covered by patches of urban areas. Lava occurs on the eastern part of the Lan Tau Island, as well as mixture of tuff and lava. Fill means the material added by human when reclaiming land from the
sea, and it mainly occurs in urban areas. Soil sampling, classification and mapping In the soil survey of Luo et al. (2007), soil sampling was conducted in a traditional manner. The dominant soil forming factors on mountain areas, parent material and relief (Fig. 2), were mainly used to guide soil sampling. After studying the information of parent material and relief, soil surveyors empirically divided Hong Kong soil into 11 kinds of soil association, and several representative soil profiles for each kind of soil association were examined. Totally, 51 soil profiles were examined, and 41 of them were on mountain areas (Fig. 1). The morphological characteristics, such as soil horizon depth, Munsell color, structure, etc., of all examined profiles were recorded. A number of physical and chemical properties of the profiles related to soil classification were analyzed in order to classify the profiles. Details about these analyses refer to Luo et al. (2007). Based on measured soil properties, the profiles were classified as Ferrosols, Ferralosols, Argosols and Cambosols according to the Chinese Soil Taxonomy (CRGCST, 2001), as indicated earlier. Ferrosols are characterized as pH < 5.0, base saturation less than 50%, having a lower activity-ferric horizon and an allitic property. Ferralosols are characterized as acid, alumino-silicate minerals accumulated, having a lower activity-ferric horizon and high clay content. Argosols are characterized as weak acid, illuvial accumulations of clay in B horizon, a low ratio of SiO2 /Al2 O3 and free iron less than 40% of total iron in soil body. Cambosols are characterized as lower accumulations of clay in B horizon and a lower ratio of SiO2 /Al2 O3 in sub-
342
Fig. 2
X. L. SUN et al.
Information related to dominated soil forming factors: parent material and relief.
soil. The 4 soil orders are similar to Ultisols, Oxisols, Alfisols and Inceptisols of the U.S. Soil Taxonomy. Soil mapping in this survey was based on the hypothesis that similar soil formation conditions will generate similar soils. Information of the two dominant soil forming factors (Fig. 2) was analyzed, and Hong Kong area was divided into a number of landformparent material groups. By assigning representative soil taxa to each kind of landform-parent material group, a soil group map was established. In the present study, this soil group map was used to produce a soil order map (Fig. 3), by assigning a soil order of a soil group to the corresponding areas of that soil group. Soil covariates Soil covariates used in this study were also from the two dominant soil forming factors: parent mate-
rial and relief (Fig. 2). Based on the relief information which was a 15 m digital elevation model (DEM) constructed based on the 1:50 000 Hong Kong topographic map using ArcGIS software, many terrain attributes were derived: elevation (Z, m), slope (S, ◦ ), aspect (A, ◦ ), profile curvature [Kp , m (100 m)−1 ], plan curvature [Kc , m (100 m)−1 ], natural logarithm of specific catchment area [ln(As ), m2 m−1 ], stream power index (SPI) and topographic wetness index (TWI). These terrain attributes are commonly used in digital soil mapping (McBratney et al., 2003), and details about them refer to Florinsky et al. (2002). Considering the non-linear variation of aspect as a circular variant, the sine of aspect (sinA) was used instead of aspect itself. Decision tree analysis and linear discriminant analysis Decision tree analysis is commonly used in soil ma-
DIGITAL SOIL ORDER MAPPING OF HONG KONG
Fig. 3
343
Soil order map based on the soil group map from the soil survey of Luo et al. (2007).
pping (Dobos et al., 2001; Scull et al., 2005; Grinand et al., 2008). For example, Scull et al. (2005) used this method to predict soil taxonomic unit from soil covariates. In this study, C5.0 (Quinlan, 2002) was used to implement decision tree analysis. A tree model must be pruned before using it to classify new samples. In C5.0, a tree model was firstly pruned by combining a leaf node or a branch that was predicted to have a relatively high misclassification rate (%, percentage of misclassified objects over all classified objects) with another one. This pruning process was applied to every branch and every leaf node. Following this pruning process, globally pruning continued to prune the tree model by examining performance of the tree model as a whole. In this study, decision tree analysis was implemented based on all of the soil covariates mentioned in the section of soil covariates. Linear discriminant analysis is also a commonly used method for digital soil mapping (McBratney et al., 2003). For example, Thomas et al. (1999) used linear discriminant analysis to predict soil classes based on terrain attributes derived from a DEM. Linear discriminant function is determined by a measure of generalized squared distance based on the pooled covariance matrix (Sinowski and Auerswald, 1999). Scores of these discriminant functions on new cases were used to classify the new cases. In this study, this analysis was implemented by the DISCRIM procedure in SAS system. Since linear discriminant analysis can not handle categorical variables, only the terrain attributes stated in the section of soil covariates were used in this analysis.
Selection of soil covariates for decision tree analysis and linear discriminant analysis Selection of soil covariates used for digital soil mapping is an “acute” problem (Lark et al., 2007). Though a lot of data mining methods can automatically choose the best variables to build a model, they do this just based on the given set of variables by modelers. For different given sets of soil covariates, different modeling results will be generated. Therefore, in this study, different sets of soil covariates were tried. Before decision tree analysis and linear discriminant analysis, variances of those calculated terrain attributes between soil orders were tested by analysis of variance (ANOVA) in SAS system. Results of this analysis could indicate the importance of a terrain attribute in differentiating soil orders, which would be suggestive for selecting variables in linear discriminant analysis and decision tree analysis. Linear discriminant analysis was conducted through the stepwise DISCRIM procedure. Selecting variables in this procedure was based on statistics of Wilks’s Lambda and average squared canonical correlation. Meanwhile, results from ANOVA were also considered in building discriminants. For the 9 soil covariates mentioned above, there would be 512 combinations of them, which obviously would generate different mapping results. Assuming that the more closely a soil covariate is related to local soil formation, the more it will influence the results of decision tree analysis and vice versa, a backward removing procedure was developed to screen the most
344
and least influential soil covariates after trying all of the 9 soil covariates in decision tree analysis to find: (1) the most influential soil covariates by removing that attribute with the least usage in decision tree analysis, and (2) the least influential soil covariates by removing that attribute with the greatest usage in decision tree analysis. Subsequently, decision tree analysis was conducted again based on the remaining attributes. These processes were repeated until the performance of decision tree analysis became very poor (indicated by misclassification rate). Based on the results of backward removing variables in decision tree analysis, as well as the results of ANOVA, the most influential soil forming factors, i.e., the most closely related factors to soil formation, were identified and were selected as the best combinations of soil covariates to build a tree model. For comparison, some other sets of soil covariates were selected to build tree models, as well as one selected according to the knowledge on local soil genesis and our experience on soil mapping. Validation and the final soil order map Considering the very limited sample size in the legacy soil data, this study used the approach proposed by Bishop and McBratney (2001), i.e., multiple jackknifing, to validate prediction results. This approach involved random separating the whole sampling dataset in a ratio into two subsets, one for prediction and one for validation. The random separation was conducted a number of times. Each time, a model was constructed based on the prediction subset and then validated based on the validation subset. Mean accuracy from this cross-validation over a number of times was used to identify the best model. In this study, the sampling dataset was separated in the ratio of about 90% to 10%, i.e., 35 samples for prediction and 6 samples for validation. Accuracy of an established model was assessed in terms of misclassification rate. This validation procedure was randomly repeated 25 times. The best model was identified in terms of the mean misclassification rate over the 25 times. The best model of those established during the above validation was selected to make the final soil order map. The criterion for selecting the best model was prediction accuracy in terms of misclassification rate and kappa statistic (Sim and Wright, 2005). Kappa statistic is an index which compares the agreement against that which might be expected by chance. It can be thought of as the chance-corrected proportional agreement (K), and possible values range from 1 (per-
X. L. SUN et al.
fect agreement) to 0 (no agreement above that expected by chance) to −1 (complete disagreement). It was computed as: K=
Po − Pc 1 − Pc
(1)
Po =
np × 100% no
(2)
C nip × nio
Pc =
i=1
no no
(3)
where Po is the observed agreement, Pc is the chance agreement, C is the number of categories of the tested samples, no is the total number of tested samples, np is the total number of correctly predicted tested samples, nip is the total number of tested samples which are predicted into the ith category, nio is the total number of tested samples which actually are the ith category. The accuracy of the final map was also assessed in terms of producer accuracy and user accuracy (Scull et al., 2005). Producer accuracy is a measure of how well the test data are classified. User accuracy is a measure of how likely a test sample classified into a given category actually belongs to that category. The flowchart of making the final soil order map is shown in Fig. 4.
Fig. 4 The flowchart for mapping soil order on mountain areas of Hong Kong, China.
DIGITAL SOIL ORDER MAPPING OF HONG KONG
RESULTS AND DISCUSSION Variances of terrain attributes between soil orders Table I showed that only elevation had significant differences (P < 0.05) between soil orders, followed by ln(As ) with difference between soil orders approaching the significant level of 0.05 (P = 0.06). The variance of slope between soil orders was also wider than the others. The significant difference in elevation between soil orders indicated that soil orders over Hong Kong had an obvious vertical zonation. Average elevations of the four soil orders were 294, 216, 137 and 61 m, respectively, for Ferrosols, Argosols, Ferralosols and Cambosols. The relatively significant difference of ln(As ) (P = 0.06) indicated the importance of water flow in soil processes, as eluviation induced by water flow plays a major role in the formation of Ferrosols, Ferralosols and Argosols (CRGCST, 2001). The influence of slope on soil processes through redistributing materials and energy over space was slightly strong, as evident by its relatively more remarkable differences between soil orders than other terrain attributes, though it was not statistically significant (P = 0.1) (Table I). Other terrain attributes seemed to be less related to the distriTABLE I Analysis of variance of terrain attributesa) between soil orders on mountain areas of Hong Kong, China Item
Z
S
sinA Kp
Kc
ln(As ) SPI TWI
F value 3.33 2.17 0.12 1.02 0.67 2.62 P > F b) 0.03* 0.10 0.94 0.39 0.57 0.06
1.00 1.71 0.40 0.18
*Significant at 0.05 level. a) Z = elevation, S = slope, A = aspect, Kp = profile curvature, Kc = plan curvature, As = specific catchment area, SPI = stream power index, TWI = topographic wetness index. b) The confidence of F value.
345
bution of the four soil orders, according to the results shown in Table I. The results of ANOVA showed that the relationship between these soil orders and a single terrain attribute was not very significant. This may be due to the fact that formation of the soil orders was not governed by a single terrain attribute, but by a number of factors, including not only terrain attributes, but also other soil forming factors and their interactions (CRGCST, 2001). However, the results of ANOVA would be useful for identifying key soil covariates in the following linear discriminant analysis and decision tree analysis. Soil covariates in linear discriminant analysis Table II showed that five terrain attributes were statistically significant in constructing linear discriminants, according to the statistics of Wilks’ lambda and average squared canonical correlation. Elevation contributed the most to the discriminations, followed by ln(As ), slope, TWI and profile curvature. This sequence was the same as that presented in Table I because stepwise DISCRIM and ANOVA shared the same way to assess contributions of the terrain attributes, i.e., the way of analyzing variances of the terrain attributes among soil orders. However, ANOVA only considered one single terrain attribute at a time, while stepwise DISCRIM considered several at a time. Thus, more terrain attributes appeared as significant variables in linear discriminant analysis compared with ANOVA. Two sets of soil covariates were selected to model discriminant functions. The first set (DI) was selected according to the results of ANOVA that showed the most significant influences of elevation, ln(As ) and slope on soil formation. Therefore, the first set (DI) comprised these three terrain attributes. The second
TABLE II Results of linear discriminant analysis for the five terrain attributesa) through the stepwise DISCRIM procedure Step
Terrain attribute
Wilks’ lambda
P < Wilks’ lambda
Average squared canonical correlation (ASCC)
P < ASCC
1 2 3 4 5
Z ln(As ) S TWI Kp
0.78 0.65 0.55 0.49 0.47
0.03* 0.01* 0.01* 0.01* 0.03*
0.07 0.13 0.17 0.21 0.22
0.03* 0.01* 0.009** 0.01* 0.02*
*, **Significant at P = 0.05 and P = 0.01 levels, respectively. a) Z = elevation, As = specific catchment area, S = slope, TWI = topographic wetness index, Kp = profile curvature.
346
X. L. SUN et al.
results suggested that parent material, elevation, profile curvature and ln(As ) were the most influential soil covariates on decision tree modeling of the local soillandscape relationships. M1 to M8 listed in the middle part of Table III shows the results of removing the most important soil covariates. Number of leaf nodes and misclassification rates varied dramatically when parent material, elevation, slope, profile curvature, ln(As ), plan curvature and TWI were removed one by one. Comparatively, these indexes varied more dramatically in the middle part of Table III than in the upper part, indicating a greater change in the performances and structures of modeled trees. It reflected that the more important a variable was, the more it influenced tree modeling. As a consequence, plan curvature and SPI were determined to be the least influencing variables in decision tree modeling of the local soil-landscape relationships. The above results showed that slope was removed very early during the process of removing both the most and least important variables. Thus, it was not easy to determine whether slope was important or not. However, results of ANOVA in Table I showed that
set (DII) was selected according to the results of linear discriminant analysis in Table II, comprising elevation, ln(As ), slope, TWI and profile curvature. Two discriminant functions were modeled based on these two sets of soil covariates. On the whole training dataset, their misclassification rates were 46% and 44%, respectively. Soil covariates in decision tree analysis Table III shows the results from implementing the procedure of backward removing variables in decision tree analysis. L1 to L8 listed in the upper part of Table III shows the results of removing the least important soil covariates. In the first five trees, SPI, slope, TWI and plan curvature were removed one by one. All of the misclassification rates of these trees were 0%, and changes in the number of leaf nodes were few. The usages of most soil covariates also changed slightly. In the last three trees, sine of aspect and ln(As ) were removed one by one. Both misclassification rates and usages of the remaining soil covariates increased slightly, with the number of leaf nodes changed a little more. These TABLE III
Decision tree modeling results of removing the least and most important soil covariatesa) Setb)
Parent material
Z
S
sinA
Kp
Kc
TWI
SPI
ln(As )
Leaf nodes
Misclassification rate %
L1 L2 L3 L4 L5 L6 L7 L8
100 100 100 100 100 100 100 100
63 63 63 63 63 63 83 83
M1 M2 M3 M4 M5 M6 M7 M8
100 -
LM KE
100 100
a)
10 10 -
17 17 17 27 27 27 -
63 63 63 71 71 71 73 73
63 100 -
10 27 100 -
17 46 32 78 61 5 22 7
63 68 44 100 -
76 63
20 10
-
63 68
Least important 12 12 0 12 12 12 12 24 12 Most important 12 12 0 20 63 7 39 0 15 95 73 0 83 63 0 100 100 76 100 100 100 76 -
39
59
-
39 39 39 39 39 39 49 -
17 17 18 18 18 18 19 16
0.0 0.0 0.0 0.0 0.0 2.4 2.4 9.8
39 68 71 61 100 -
17 14 14 16 16 10 9 7
0.0 7.3 12.2 7.3 2.4 24.4 31.7 36.6
49
19 19
0.0 0.0
-
Z = elevation, S = slope, A = aspect, Kp = profile curvature, Kc = plan curvature, TWI = topographic wetness index, SPI = stream power index, As = specific catchment area. b) LM: removing the least and most important soil covariates together; KE: knowledge and experience in soil mapping.
DIGITAL SOIL ORDER MAPPING OF HONG KONG
slope was relatively closely related to distributions of soil orders (P = 0.1), and in fact slope was very frequently used for soil mapping (McBratney et al., 2003). Therefore, slope was finally used in the subsequent prediction of soil orders. Contrary to slope, sinA was removed very late during the both processes, but results of ANOVA showed that this soil covariate was scarcely related to distributions of soil orders (P = 0.9). Therefore, sinA was not used in the subsequent prediction of soil orders. Based on the above analysis, the first set of soil covariates for decision tree modeling was designed to include all of the most influential soil covariates: parent material, elevation, slope, profile curvature and ln(As ), i.e., LM listed in Table III. With the objective of generating a good quality soil order map for this area, other three sets of soil covariates, i.e., L1 , L5 and KE in Table III, were also selected to construct decision tree models. L1 and L5 were chosen because they were two extremes in the process of backward removing the most important variables. L1 contained all soil covariates, while L5 contained the fewest but keeping the lowest misclassification rate. KE was established according to the common knowledge about local soil genesis and our experience in soil mapping. As mentioned earlier, parent material and relief played dominant roles in the local soil development. McSweeney et al. (1994) contended that on the watershed scale, soil formation and development were mainly influenced by five environmental variables: elevation, slope, profile curvature, plan curvature and TWI. However, ln(As ), depicting the contribution area of water flowing across a point, would be more appropriate in reflecting the influence of water on soil processes than TWI, which depicts stationary water content in soil. This is because water flow passing through soil bodies or on soil surfaces affects soil more directly and more strongly than stationary water in soil bodies (Florinsky et al., 2002; Thompson et al., 2006). Particularly formations of Ferrolsols, Ferralosols and Argosols are more related to water movement than water content (CRGCST, 2001). Furthermore, results of ANOVA shown in Table I also indicated that ln(As ) was more related to distributions of soil orders than TWI. As a consequence, KE was designed to comprise parent material, elevation, slope, profile curvature, plan curvature and ln(As ). Based on the four designed sets of soil covariates, four decision trees were modeled respectively. Performances of these four tree models are presented in Table III. Misclassification rates of them were all 0%, and
347
their numbers of leaf nodes were close, i.e., 17, 18, 19 and 19. Mapping soil orders based on selected soil covariates Results from validating the above selected two discriminant functions and four decision tree models using the multiple jackknifing validation approach are summarized in Table IV. Generally, the two discriminant functions performed better than the four tree models, although they did not explore the usefulness of parent material for predicting soil orders. Between the two discriminant functions, model DI performed much better than model DII, reflecting that a more straightforward but good enough discriminant function was much more useful. Among the four tree models, model LM which was selected through backward removing procedure performed best, followed by model KE which was selected according to knowledge and experience on soil mapping, model L5 and model L1 . Results from validating the four tree models reflected that: 1) it was important to select soil covariates for tree modeling, even though decision tree analysis can automatically choose best variables for modeling; 2) backward removing procedure improved tree modeling, as well as knowledge and experience on soil mapping; 3) a more straightforward but good enough tree model was much more useful, like that for discriminant functions. TABLE IV Summaries on misclassification rate of selected two discriminant functions (DI and DII) and four decision tree models (L1 , L5 , LM and KE) using the multiple jackknifing validation approach Model Mean Minimum Maximum Standard deviation % DI DII L1 L5 LM KE
49 59 71 69 67 68
33 33 33 33 33 33
83 100 100 100 100 100
14.55 16.75 15.18 14.29 15.64 14.77
According to mean misclassification rates listed in Table IV, model DI with the lowest mean misclassification rate was identified as the best model for producing the final soil order map. Among all the 25 models which were established using this model during jackknifing validation, the one having the lowest misclassification rate (33%) and the highest kappa statistic
348
X. L. SUN et al.
(0.63) was selected to generate the final soil order map (Fig. 5). User accuracy and producer accuracy of this map from jackknifing validation are presented in Table V. The table shows that the final map was very accurate to reflect the distribution of Ferrosols and Ferralosols. Distribution of Argosols was very poorly predicted and distribution of Cambosols was the worst predicted. Weighted by the area percentages of soil orders on the map, average producer accuracy and average user accuracy of this map were 98.7% and 82.6%, respectively, showing that the final soil order map had a high accuracy. The final digital soil order map (Fig. 5) was very different from the traditional one (Fig. 3). Apparently, the digital soil order map was much more detailed in reflecting the distribution of soil orders. Comparison between the two maps showed that Cambosols in the digital map were much less than in the traditional map and were separately distributed among other soil orders.
Ferrolosols were also reduced substantially on the digital map compared with the traditional map. The two soil orders seemed to be replaced by Argosols in Fig. 5, while some Argosols were replaced by Ferralosols. Totally, the digital map reproduced only 18.4% of the traditional one. This big gap between them is due first to the different models used to produce the two maps and second but more crucially, it is most likely due to the limited sample size of the legacy soil data on which this study was based. CONCLUSIONS Linear discriminant models performed better than decision tree models. ANOVA was very useful for selecting soil covariates for building a discriminant. The removing procedure developed in this study for selecting soil covariates for tree modeling was useful to improve the performance of tree modeling. This was also true
Fig. 5 The final soil order map of Hong Kong predicted by model DI, which included three soil covariates, i.e., elevation, natural logarithm of specific catchment area and slope, for modeling discriminant functions. TABLE V Error matrix of final soil order map Predicted soil order
Argosol Cambosol Ferralosol Ferrosol Producer accuracy (%)
Observed soil order
User accuracy
Argosols
Cambosols
Ferralosols
Ferrosols
1 0 0 0 100
2 0 0 0 0
0 0 1 0 100
0 0 0 2 100
% 33 0 100 100
DIGITAL SOIL ORDER MAPPING OF HONG KONG
for knowledge and experience on soil mapping. Clearly, results of this study highlighted the importance of selecting soil covariates in modeling for soil mapping, and suggested the usefulness of the methods used in this study for selecting soil covariates. A soil order map was finally made using the best model established in this study, with a high accuracy in terms of misclassification rate, kappa statistics, producer accuracy and user accuracy. Comparison between the digital soil order map and the traditional one showed that the former just reproduced 18.4% of the later. The big gap between them can be attributed to the very limited sample size of the legacy soil data on which this study was based. Thus, more soil samples in future soil surveys for this area will be necessary in order to improve the final soil order map, as well as to produce more detailed soil maps. ACKNOWLEDGEMENT The authors are grateful to anonymous reviewers for their valuable comments and suggestions, and Dr. A. O. W. Leung at the Croucher Institute for Environmental Sciences, Hong Kong Baptist University, Hong Kong Special Administrative Region, China for editing the language of this manuscript. REFERENCES Au, S. W. C. 1998. Rain-induced slope instability in Hong Kong. Eng. Geol. 51: 1–36. Bishop, T. F. A. and McBratney, A. B. 2001. A comparison of prediction methods for the creation of field-extent soil property maps. Geoderma. 103: 149–160. Cooperative Research Group on Chinese Soil Taxonomy (CRGCST). 2001. Chinese Soil Taxonomy. Science Press, Beijing and New York. Dobos, E., Micheli, E., Baumgardner, M. F., Biehl, L. and Helt, T. 2000. Use of combined digital elevation model and satellite radiometric data for regional soil mapping. Geoderma. 97: 367–391. Dobos, E., Montanarella, L., N`egre, T. and Micheli, E. 2001. A regional scale of soil mapping approach using intergrated AVHRR and DEM data. Int. J. Appl. Earth Observ. Geoinf. 3: 30–42. Florinsky, I. V., Eilers, R. G., Manning, G. R. and Fuller, L. G. 2002. Prediction of soil properties by digital terrain modelling. Environ. Modell. Softw. 17: 295–311. Grant, C. J. 1962. The Soils and Agriculture of Hong Kong. Government Press, Hong Kong. Grinand, C., Arrouays, D., Laroche, B. and Martin, M. P. 2008. Extrapolating regional soil landscapes from an existing soil map: Sampling intensity, validation pro-
349 cedures, and integration of spatial context. Geoderma. 143: 180–190. Henderson, B. L., Bui, E. N., Moran, C. J. and Simon, D. A. P. 2005. Australia-wide predictions of soil properties using decision trees. Geoderma. 124: 383–398. IUSS Working Group WRB. 2006. World reference base for soil resources 2006. World Soil Resources Reports No. 103. FAO, Rome. Kempen, B., Brus, D. J., Heuvelink, G. B. M. and Stoorvogel, J. J. 2009. Updating the 1:50 000 Dutch soil map using legacy soil data: A multinomial logistic regression approach. Geoderma. 151: 311–326. Lagacherie, P. 2008. Digital soil mapping: a state of the art. In Hartemink, A. E., McBratney, A. B. and MendoncaSantos, M. L. (eds.) Digital Soil Mapping with Limited Data. Springer, Dordrecht. pp. 3–14. Lark, R. M., Bishop, T. F. A. and Webster, R. 2007. Using expert knowledge with control of false discovery rate to select regressors for prediction of soil properties. Geoderma. 138: 65–78. Li, X., Poon, C. and Liu, P. S. 2001. Heavy metal contamination of urban soils and street dusts in Hong Kong. Appl. Geochem. 16: 1361–1368. Luo, Y. M., Li, Z. G., Wu, L. H., Wu, S. C., Zhang, G. L., Zhou, S. L., Zhao, Y. G., Zhao, Q. G., Wong, M. H. and Zhang, H. B. 2007. Hong Kong Soils and Environment (in Chinese). Science Press, Beijing. McBratney, A. B., Mendonca-Santos, M. L. and Minasny, B. 2003. On digital soil mapping. Geoderma. 117: 3–52. McSweeney, K., Gessler, P. E., Slater, B. K., Hammer, R. D., Bell, J. and Petersen, G. W. 1994. Towards a new framework for modeling the soil-landscape continuum. In Amundson, R, G., Harden, J. W. and Singer, M. (eds.) Factors of Soil Formation: A Fiftieth Anniversary Retrospective. Soil Science Society of America, Madison, WI. pp. 127–145. Quinlan, J. R. 2002. Data Mining Tools See5 and C5.0. Available online at http://www.rulequest.com/see5info.html (verified on January, 2011). Scull, P., Franklin, J. and Chadwick, O. A. 2005. The application of classification tree analysis to soil type prediction in a desert landscape. Ecol. Model. 181: 1–15. Sim, J. and Wright, C. C. 2005. The kappa statistic in reliability studies: use, interpretation, and sample size requirements. Phys. Ther. 85: 257–268. Sinowski, W. and Auerswald, K. 1999. Using relief parameters in a discriminant analysis to stratify geological areas with different spatial variability of soil properties. Geoderma. 89: 113–128. Stoorvogel, J. J., Kempen, B., Heuvelink, G. B. M. and de Bruin, S. 2009. Implementation and evaluation of existing knowledge for digital soil mapping in Senegal. Geoderma. 149: 161–170. The Hong Kong Geological Survey. 1999. Geological Map
350 of Hong Kong in 1:200 000. 2nd Edition. Survey and Mapping Office, Lands Department, the Government of Hong Kong Special Administrative Region, Hong Kong. Thomas, A. L., King, D., Dambrine, E., Couturie, A. and Roque, J. 1999. Predicting soil classes with parameters derived from relief and geologic materials in a sandstone region of the Vosges mountains (Northeastern France). Geoderma. 90: 291–305. Thompson, J. A., Pena-Yewtukhiw, E. M. and Grove, J. H.
X. L. SUN et al. 2006. Soil-landscape modeling across a physiographic region: topographic patterns and model transportability. Geoderma. 133: 57–70. Warren-Rhodes, K. and Koenig, A. 2001. Ecosystem appropriation by Hong Kong and its implications for sustainable development. Ecol. Econ. 39: 347–359. Zhu, A. X., Band, L., Vertessy, R. and Dutton, B. 1997. Derivation of soil properties using a soil land inference model (SoLIM). Soil Sci. Soc. Am. J. 61: 523–533.