Mapping Fragmented Impervious Surface Areas Overlooked by Global Land-Cover Products in the Liping County, Guizhou Province, China

: Imperviousness is an important indicator for monitoring urbanization and environmental changes, and is evaluated widely in urban areas, but not in rural areas. An accurate impervious surface area (ISA) map in rural areas is essential to achieve environmental conservation and sustainable rural development. Global land-cover products such as MODIS MCD12Q1, ESA CCI-LC, and Global Urban Land are common resources for environmental practitioners to collect land-cover information including ISAs. However, global products tend to focus on large ISA agglomerations and may not identify fragmented ISA extents in less populated regions. Land-use planners and practitioners have to map ISAs if it is di ﬃ cult to obtain such spatially explicit information from local governments. A common and consistent approach for rural ISA mapping is yet to be established. A case study of the Liping County, a typical rural region in southwest China, was undertaken with the objectives of assessing the global land-cover products in the context of rural ISA mapping and proposing a simple and feasible approach for the mapping. This approach was developed using Landsat 8 imagery and by applying a random forests classiﬁer. An appropriate number of training samples were distributed to towns or villages across all townships in the study area for classiﬁcation. The results demonstrate that the global land-cover products identiﬁed major ISA agglomerations, speciﬁcally at the county seat; however, other fragmented ISAs over the study area were overlooked. In contrast, the map created using the developed approach inferred ISAs across all townships with an overall accuracy of 91%. A large amount of training samples together with geographic information of towns or villages is the key suggestion to identify and map ISAs in rural areas.


Introduction
The development of impervious surface areas (ISAs) has profound impacts on both urban and rural regions [1,2]. Impervious surfaces (ISs) usually refer to artificial surfaces such as pavements and roofs, through which a fluid cannot pass [3,4]. Imperviousness was originally defined as an artificial condition that affects water quality and urban runoff [5,6]. More recently, imperviousness was recognized as a key indicator to monitor urban extent [7,8] and understand human-environment interactions [4,9]. The expansion of ISAs has the potential to substantially impact local hydrology and water environment [10], accelerate loss of biodiversity [11], contribute to climate change [12], and threaten food security [13]. Spatially explicit information on ISAs is essential to support watershed management, urban/rural planning, and land use regulation [14,15]. Monitoring ISAs is essential for all practitioners, policymakers, and researchers who are concerned with anthropogenic impacts on land [4].

Study Area
The Liping County in the Guizhou Province is a typical mountainous rural county in southwest China in terms of economic status and geography [33]. Liping County consists of 26 townships administering 33 towns and 390 villages [34] (Figure 1). Most of the towns and villages are located in relatively lower flat areas in valleys and are surrounded by farmlands. It covers approximately 4433 km 2 with a population of 393,000 in 2017 [34]. With the process of rapid urbanization in China, many local young people have migrated to eastern coastal urbanized cities, resulting in a drastic depopulation of approximately 18% from 2006 to 2017 [34]. The migrants, however, face a social barrier in becoming citizens of their destined cities owing to the household registration system, which effectively makes them unwilling to give up their rural land quota and motivates them to build new bigger houses in their original hometown [17], causing land-cover conversions from agricultural to ISAs in rural areas. Furthermore, the local government in China might have huge ambitions to develop new constructions, which may lead to excessive development of ISAs [35]. In the absence of a spatially explicit ISA map for the Liping County, concerns on efficient rural land management arise. Therefore, practitioners and policymakers require an accurate ISA map for the Liping County to prepare a sustainable land-use plan that can be utilized to prevent excessive future land developments.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 14 The Liping County in the Guizhou Province is a typical mountainous rural county in southwest China in terms of economic status and geography [33]. Liping County consists of 26 townships administering 33 towns and 390 villages [34] (Figure 1). Most of the towns and villages are located in relatively lower flat areas in valleys and are surrounded by farmlands. It covers approximately 4433 km 2 with a population of 393,000 in 2017 [34]. With the process of rapid urbanization in China, many local young people have migrated to eastern coastal urbanized cities, resulting in a drastic depopulation of approximately 18% from 2006 to 2017 [34]. The migrants, however, face a social barrier in becoming citizens of their destined cities owing to the household registration system, which effectively makes them unwilling to give up their rural land quota and motivates them to build new bigger houses in their original hometown [17], causing land-cover conversions from agricultural to ISAs in rural areas. Furthermore, the local government in China might have huge ambitions to develop new constructions, which may lead to excessive development of ISAs [35]. In the absence of a spatially explicit ISA map for the Liping County, concerns on efficient rural land management arise. Therefore, practitioners and policymakers require an accurate ISA map for the Liping County to prepare a sustainable land-use plan that can be utilized to prevent excessive future land developments.

Global Land-Cover Products
To assess the performance of global land-cover products for ISA mapping in the Liping County, three global land-cover datasets were utilized: MODIS/Terra and Aqua Combined Land Cover Type Yearly Global 500 m SIN Grid Version 006 (MODIS MCD12Q1) [24], ESA Climate Change Initiative Land Cover Project (ESA CCI-LC) [25], and high-resolution multi-temporal mapping of Global Urban Land (GUL) [8], with spatial resolutions of 500 m, 300 m, and 30 m, respectively. Considering the mapped years that were available across these products, the year 2015 was selected to ensure consistency in comparison. Imperviousness-related classes are labeled as "urban and built-up" in

Global Land-Cover Products
To assess the performance of global land-cover products for ISA mapping in the Liping County, three global land-cover datasets were utilized: MODIS/Terra and Aqua Combined Land Cover Type Yearly Global 500 m SIN Grid Version 006 (MODIS MCD12Q1) [24], ESA Climate Change Initiative Land Cover Project (ESA CCI-LC) [25], and high-resolution multi-temporal mapping of Global Urban Land (GUL) [8], with spatial resolutions of 500 m, 300 m, and 30 m, respectively. Considering the mapped years that were available across these products, the year 2015 was selected to ensure consistency in comparison. Imperviousness-related classes are labeled as "urban and built-up" in Remote Sens. 2020, 12, 1527 4 of 14 MODIS MCD12Q1, "urban areas" in ESA CCI-LC, and "urban" in GUL. The legends of "urban" related classes in the three global products are all associated with impervious and/or artificial surfaces, providing a similarity of definition between them (Table 1). Therefore, these comparable classes were considered as ISAs, and the classification was simplified in these global land-cover products into two classes: "ISs" and "other" (i.e., land covers except for ISs).

Satellite Imagery
Landsat archive from United States Geological Survey (USGS)-a free, open, and long archiving satellite data-is widely used in land cover/use monitoring [37]. As the Landsat archive contains data over the last 30 years, it has a huge potential for use in rural ISA mapping for the past as well. Landsat 8 at 30 m spatial resolution for the year 2015 was selected for a consistent comparison with the other global land-cover products. To demonstrate the availability of the Landsat archive, the Landsat 8 Surface Reflectance Tier 1 was used to produce the ISA map. This dataset is the atmospherically corrected surface reflectance data from the Landsat 8 OLI/TIRS sensors [38]. Using the cloud mask provided by quality control band values found in the dataset, the annual median composite of imagery in 2015 with bands 1-7 (b1: ultra-blue, b2: blue, b3: green, b4: red, b5: near-infrared, b6: shortwave infrared 1, and b7: shortwave infrared 2) at the spatial resolution of 30 m was preprocessed on the Google Earth Engine Platform (GEE) [39].

Supervised Classification
The random forests (RF) classifier was used to classify land covers from Landsat 8 surface reflectance imagery into two classes: "ISs" and "other". "IS(s)" is defined as an artificial land covered with impenetrable materials such as asphalt, concrete, tile, brick, and stone, corresponding with the definitions used in the global land-cover products (Table 1) [8,24,25]. "Other" refers to a land covered with surfaces other than impervious surfaces such as vegetation, soils, and water.
RF, a non-parametric and robust supervised classification method, is widely applied in land-cover classification [40,41]. RF is developed from the decision tree classifier by randomly selecting a subset of variables through replacement [42]. A decision tree divides a dataset into smaller subdivisions based on the defined threshold at each branch [43]. However, a decision tree is not always robust as different training data may result in different trees. RF has been developed to address this issue by randomly selecting a subset of training data and growing many trees, like that of a "forest".

Training Data Collection
A training sample was collected for the RF classifier with the use of a false-color composite image with the combination of RGB = 7-6-4 bands to highlight the ISAs visually. For example, in Figure 2, the colors represent different land types: purple, white, green, and cyan indicate roof materials, pavements, forests, and cropland, respectively. Thus, the pixels were randomly selected and manually labeled with colors. The accuracy of classification results improves with the increase in the size of the training sample [44]. The training sample collection fulfills sufficiency criteria as a considerable amount of data was collected using a polygon-based sampling method with sufficient spatial coverage within the study area. The polygon-based training method performs better than the pixel-based method if the classes are spectrally heterogeneous [45]. Thus, a total of 624 sample polygons, including 217 for "ISs" and 407 for "other" (Figure 3), were created. It was assumed that the ISAs were distributed in all townships. Sample collection was designed such that it included all townships, that is, each township has at least one polygon. The existence of villages and towns was supplementally confirmed using OpenStreetMap. cover classification [40,41]. RF is developed from the decision tree classifier by randomly selecting a subset of variables through replacement [42]. A decision tree divides a dataset into smaller subdivisions based on the defined threshold at each branch [43]. However, a decision tree is not always robust as different training data may result in different trees. RF has been developed to address this issue by randomly selecting a subset of training data and growing many trees, like that of a "forest".

Training Data Collection
A training sample was collected for the RF classifier with the use of a false-color composite image with the combination of RGB = 7-6-4 bands to highlight the ISAs visually. For example, in Figure 2, the colors represent different land types: purple, white, green, and cyan indicate roof materials, pavements, forests, and cropland, respectively. Thus, the pixels were randomly selected and manually labeled with colors. The accuracy of classification results improves with the increase in the size of the training sample [44]. The training sample collection fulfills sufficiency criteria as a considerable amount of data was collected using a polygon-based sampling method with sufficient spatial coverage within the study area. The polygon-based training method performs better than the pixel-based method if the classes are spectrally heterogeneous [45]. Thus, a total of 624 sample polygons, including 217 for "ISs" and 407 for "other" (Figure 3), were created. It was assumed that the ISAs were distributed in all townships. Sample collection was designed such that it included all townships, that is, each township has at least one polygon. The existence of villages and towns was supplementally confirmed using OpenStreetMap.

Model Setup
The RF classifier was implemented using randomForest function in the randomForest package of R software [46]. The variables for the RF classifier are shown in Table 2. The predictor variables were seven bands from Landsat 8 Surface Reflectance Tier 1, which were processed on GEE, and the response variable was a binary variable defined by the land covers: "ISs" (coded as 1) and "other" (coded as 0). RF has a good performance in coping with complex interactions and highly correlated predictor variables by randomly selecting a subset of predictor variables [42]. As a result, the trees are less similar and less correlated, and the final prediction can be more accurate. It is noted that the multi.collinear function in rfUtilities package confirmed no multi-collinearity among the predictor variables [47,48].

Model Setup
The RF classifier was implemented using randomForest function in the randomForest package of R software [46]. The variables for the RF classifier are shown in Table 2. The predictor variables were seven bands from Landsat 8 Surface Reflectance Tier 1, which were processed on GEE, and the response variable was a binary variable defined by the land covers: "ISs" (coded as 1) and "other" (coded as 0). RF has a good performance in coping with complex interactions and highly correlated predictor variables by randomly selecting a subset of predictor variables [42]. As a result, the trees are less similar and less correlated, and the final prediction can be more accurate. It is noted that the multi.collinear function in rfUtilities package confirmed no multi-collinearity among the predictor variables [47,48]. Two essential parameters ntree and mtry need to be set for the RF classifier [41]; ntree indicates the number of trees to grow and mtry denotes the number of predictor variables randomly sampled at each split of tree node. Figure 4 illustrates the result of some preliminary tests based on out-of-bag (OOB) errors with different parameters. The OOB error is an effective internal estimate to monitor the error of RF, which helps tune parameters [42]. This test indicated that mtry = 2 with all the ntree Remote Sens. 2020, 12, 1527 7 of 14 numbers resulted in the lowest OOB error, and larger ntree tended to reduce the OOB error. From this exploratory test, we set the parameters of ntree as 500 and mtry as 2 for the RF classifier.
Two essential parameters ntree and mtry need to be set for the RF classifier [41]; ntree indicates the number of trees to grow and mtry denotes the number of predictor variables randomly sampled at each split of tree node. Figure 4 illustrates the result of some preliminary tests based on out-of-bag (OOB) errors with different parameters. The OOB error is an effective internal estimate to monitor the error of RF, which helps tune parameters [42]. This test indicated that mtry = 2 with all the ntree numbers resulted in the lowest OOB error, and larger ntree tended to reduce the OOB error. From this exploratory test, we set the parameters of ntree as 500 and mtry as 2 for the RF classifier.

Accuracy Assessment
The ISA classification map was assessed by a confusion matrix and the overall, precision, and recall accuracy measures [49], which were calculated as follows: where tp, fp, tn, and fn indicate the number of cases of true positive, false positive, true negative, and false negative, respectively. Here, the class of ISs is positive, while the class of "other" is negative in this study. Stratified random sampling was employed for collecting validation sample points over the study area [50]. The targeted sample size is calculated as follows: where ( ) is the standard error of the estimated overall accuracy that is required, is the mapped proportion of the area of the class , and is the standard deviation of stratum . A target standard error for the overall accuracy was set at 0.02, and it was conjectured that the recall accuracy will be 0.85 for "ISs" and 0.9 for "other". On this basis, the total validation sample size was determined as

Accuracy Assessment
The ISA classification map was assessed by a confusion matrix and the overall, precision, and recall accuracy measures [49], which were calculated as follows: Overall accuracy = tp + tn tp + tn + f p + f n where tp, fp, tn, and fn indicate the number of cases of true positive, false positive, true negative, and false negative, respectively. Here, the class of ISs is positive, while the class of "other" is negative in this study. Stratified random sampling was employed for collecting validation sample points over the study area [50]. The targeted sample size is calculated as follows: where S(ô) is the standard error of the estimated overall accuracy that is required, w i is the mapped proportion of the area of the class i, and s i is the standard deviation of stratum i. A target standard error for the overall accuracy was set at 0.02, and it was conjectured that the recall accuracy will be 0.85 for "ISs" and 0.9 for "other". On this basis, the total validation sample size was determined as 227 for this study. Proportional allocation of the validation sample size is a popular sample allocation strategy, but produces inaccurate recall accuracy for rare classes [50]. Hence, to investigate "ISs" more intensively, sample points were allocated equally, as 113 for "ISs" and 114 for "other".
The validation sample was visually recorded and collected from very fine spatial resolution imagery in Google Earth to assure the independence of the accuracy assessment from the training scheme. It should be noted that very fine spatial resolution images of some parts of the Liping County were not available in Google Earth for 2015. To deal with this issue, images from 2013 and 2017 were used complementarily to conjecture the 2015 land covers in these areas by confirming a stable land cover over the period. In this case, it is a reasonable approach to obtain sample data in regions that are not covered by the available very fine spatial resolution images, as IS tends to be developed quickly and is not often converted into other land-cover types [51].

Results
ISA maps for the Liping County produced using the global land-cover products (MODIS MCD12Q1, ESA CCI-LC, and GUL) and the RF classifier are shown in Figure 5. The global products depicted ISAs at the county seat, the largest ISA in the Liping County, but underestimated them in the other regions. ESA CCI-LC represented a few distributions of small ISAs in the eastern regions; however, the map was far from accurate. As illustrated in Figure 1, there are 390 villages and 33 towns located in the Liping County, and each township has at least four villages according to the statistical yearbook published by the Liping County [34]. As the three "urban" classes (Table 1) are defined as all major artificial ISs, the ISAs are underestimated by the global land-cover products in almost all townships. Moreover, the sizes of ISA agglomerations at the county seat are different among the three global products. For example, MODIS MCD12Q1 identified the smallest total ISA (375 ha) among the three global products, while ESA CCI-LC showed the largest ISA (517 ha). GUL inferred a total of 467 ha of ISAs. As the estimation errors were clearly significant, accuracy assessments were not conducted for the global product maps. In contrast, RF inferred a wide range of sizes of ISAs in the semi-urban area around the county seat and rural areas across all the townships.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 14 227 for this study. Proportional allocation of the validation sample size is a popular sample allocation strategy, but produces inaccurate recall accuracy for rare classes [50]. Hence, to investigate "ISs" more intensively, sample points were allocated equally, as 113 for "ISs" and 114 for "other". The validation sample was visually recorded and collected from very fine spatial resolution imagery in Google Earth to assure the independence of the accuracy assessment from the training scheme. It should be noted that very fine spatial resolution images of some parts of the Liping County were not available in Google Earth for 2015. To deal with this issue, images from 2013 and 2017 were used complementarily to conjecture the 2015 land covers in these areas by confirming a stable land cover over the period. In this case, it is a reasonable approach to obtain sample data in regions that are not covered by the available very fine spatial resolution images, as IS tends to be developed quickly and is not often converted into other land-cover types [51].

Results
ISA maps for the Liping County produced using the global land-cover products (MODIS MCD12Q1, ESA CCI-LC, and GUL) and the RF classifier are shown in Figure 5. The global products depicted ISAs at the county seat, the largest ISA in the Liping County, but underestimated them in the other regions. ESA CCI-LC represented a few distributions of small ISAs in the eastern regions; however, the map was far from accurate. As illustrated in Figure 1, there are 390 villages and 33 towns located in the Liping County, and each township has at least four villages according to the statistical yearbook published by the Liping County [34]. As the three "urban" classes (Table 1) are defined as all major artificial ISs, the ISAs are underestimated by the global land-cover products in almost all townships. Moreover, the sizes of ISA agglomerations at the county seat are different among the three global products. For example, MODIS MCD12Q1 identified the smallest total ISA (375 ha) among the three global products, while ESA CCI-LC showed the largest ISA (517 ha). GUL inferred a total of 467 ha of ISAs. As the estimation errors were clearly significant, accuracy assessments were not conducted for the global product maps. In contrast, RF inferred a wide range of sizes of ISAs in the semi-urban area around the county seat and rural areas across all the townships.  Table 3 summarizes the results from the four maps obtained using the RF classifier and global land-cover products. The RF method estimated that the ISAs represented 1.3% of the land cover in the Liping County, an order of magnitude higher than that of the global products, all of which identified ≤0.12%. Table 3. Total IS areas and their proportion of total area in the Liping County estimated from the global land-cover products and the random forests classifier for 2015.

MODIS MCD12Q1 ESA CCI-LC GUL
ISAs (ha) 5797 375 517 467 The proportion of ISAs in the total area (%) 1.3 0.08 0.12 0.11 Figure 6 presents the variable importance of the RF classifier. Mean decrease accuracy indicates how much the accuracy decreases if a variable is excluded, while mean decrease Gini suggests the decrease of node impurity when a variable is selected. Both of the variable importance measures indicated that band 5, near-infrared (0.851-0.879 µm), contributed the most to the final accuracy, which is reasonable as vegetation, water, and impervious areas are likely to be distinct in this wavelength range. Table 4 summarizes the OOB assessment of the RF classifier. The OOB-based overall accuracy achieved 0.99, indicating the goodness of fit of the model. The confusion matrix assessed by independent validation samples indicated that the RF map achieved an overall accuracy of 0.91, a precision of 0.84, and a recall of 0.97, indicating relatively high accuracy (Table 5). These assessments indicated a similar tendency of the accuracy, suggesting robustness of the model.  Table 3 summarizes the results from the four maps obtained using the RF classifier and global land-cover products. The RF method estimated that the ISAs represented 1.3% of the land cover in the Liping County, an order of magnitude higher than that of the global products, all of which identified ≤0.12%. Table 3. Total IS areas and their proportion of total area in the Liping County estimated from the global land-cover products and the random forests classifier for 2015.

MODIS MCD12Q1 ESA CCI-LC GUL
ISAs (ha) 5797 375 517 467 The proportion of ISAs in the total area (%) 1.3 0.08 0.12 0.11 Figure 6 presents the variable importance of the RF classifier. Mean decrease accuracy indicates how much the accuracy decreases if a variable is excluded, while mean decrease Gini suggests the decrease of node impurity when a variable is selected. Both of the variable importance measures indicated that band 5, near-infrared (0.851-0.879 μm), contributed the most to the final accuracy, which is reasonable as vegetation, water, and impervious areas are likely to be distinct in this wavelength range.  Table 4 summarizes the OOB assessment of the RF classifier. The OOB-based overall accuracy achieved 0.99, indicating the goodness of fit of the model. The confusion matrix assessed by independent validation samples indicated that the RF map achieved an overall accuracy of 0.91, a precision of 0.84, and a recall of 0.97, indicating relatively high accuracy (Table 5). These assessments indicated a similar tendency of the accuracy, suggesting robustness of the model. Table 4. Out-of-bag (OOB) confusion matrix for the random forests classifier.

Discussion
The results show considerable differences between global land-cover products and an RF classification map in terms of the ISA estimation in the Liping County in China. All global land-cover products underestimated the extent of the ISAs and only inferred a few ISA agglomerations including the largest ISA at the county administrative center. In contrast, the RF estimated the distribution of ISAs in the lower and flat areas of valleys in a wide range of regions including all the townships. The classification accuracy also indicated an acceptable result. Even though a relatively precise ISA map was produced, some classification errors still exist. The aim of this study was not to develop the best classification model in terms of accuracy, but to demonstrate a simple and feasible classification method that can be utilized by experts as well as non-experts who have an interest in producing ISA maps of rural areas. The RF classifier was found to provide acceptable performance in mapping ISAs in the Liping County. It is evident that such a simple approach performs better than the existing global land-cover products in the context of mapping ISAs in less populated areas.

Reasons for Global Products Overlooking Rural ISAs
Unlike the large extents of ISs in urban areas, the ISAs in rural areas are likely to be small and fragmented. Global land-cover products tend to have coarse resolution and mainly focus on the large urban ISA extents. These characteristics might be the foremost reasons for the lack of recognition of ISAs in rural areas by the global land-cover products. As summarized in Table 1, the MODIS MCD12Q1 product employed supervised classification [24], while ESA CCI-LC uses spectral unsupervised and supervised algorithms that input variables associated with spectral properties and then integrates the two classification results based on objective rules for land-cover classification [25]. Both products employ per-pixel classification with a coarse spatial resolution (over 300-500 m), and hence, some rural ISAs may be smaller than such a minimum mapping unit, resulting in a failure of mapping. Vintrou et al. [52] reported a similar problem when using coarse-resolution global land-cover products to map the small extents of croplands in Mari. Therefore, spatial resolution is important in determining the mapping performance. In this sense, GUL has a resolution of 30 m [8] and the accuracy pertaining to spatial resolution is less. However, the protocols used for GUL may not perform well in mapping rural ISAs. This product mainly focuses on large urban agglomerations with the use of The Defense Meteorological Program (DMSP) operational line-scan system (OLS) (DMSP-OLS) nighttime lights, whose spatial resolution is 1 km to mask urban coverage, leading to the disregard of small impervious agglomerations in rural areas. Moreover, as global land-cover products collect training or validation data at a global level, they are more likely to neglect local scale entities. In contrast to the global products, in this study, one particular area was focused on; thus, reference data that were distributed more widely throughout the area of interest could be collected rather than focusing only on large urban ISA extents.

Recommendations for Rural ISA Mapping
The simple training data collection strategy with a RF classifier utilized in this study was able to improve rural ISA mapping. Compared with the global land-cover products, the mapping results were greatly improved using the medium resolution (30 m) Landsat imagery. This strategy is proposed to achieve rural ISA mapping based on comprehensively collected training data, particularly by collecting samples at every administrative unit by referring to the actual location of villages or towns and creating a large training sample to achieve the required classification accuracy. The suitability of algorithms can vary in terms of the study area and satellite imagery, and the selection of algorithms is not restricted. The RF method could be replaced by other supervised classification algorithms such as maximum likelihood (ML), support vector machines (SVMs), or any other classification method. Although the accuracy of the RF classification was the best throughout our experiments (Table A1), all of the other classification methods were able to estimate the fragmented ISAs in the Liping County. Non-parametric algorithms such as RF and SVM tend to outperform the parametric algorithm like ML in land-cover classification [41]. As SVM is more flexible and complicated, requiring a careful parameter selection, the RF classifier was used to take advantage of its robustness and accuracy for this case study. This strategy can be implemented in other rural regions that have small and fragmented ISAs. In future studies, rural ISA mapping in other areas using this proposed strategy will be explored.

Conclusions
This study highlights the importance of mapping ISAs in rural areas (or less populated areas) that have rarely been focused on. It was found that most of the rural ISAs are not depicted in the global land-cover products because they overlook small and fragmented ISAs. This study demonstrates a simple and feasible training sample strategy with a RF classifier to obtain accurate rural ISA maps. When compared with the global land-cover products, this simple strategy significantly improved the ISA map obtained for the Liping County, a typical rural area in southwest China. Owing to its simplicity, this approach can be applied to any other rural regions. The spatial extent and distribution of the fragmented ISAs in rural areas is an important indicator to monitor the anthropogenic impact on natural resources. This approach provides a valuable means of mapping the fragmented ISAs in rural areas for efficient rural land development and planning.
Author Contributions: J.Z.: data curation, writing-original draft preparation, methodology, software, formal analysis, visualization, investigation. N.T.: conceptualization, methodology, validation, supervision, writing-reviewing and editing, project administration, funding acquisition. All authors have read and agreed to the published version of the manuscript.
Funding: This research were funded by the Kyoto University Foundation and the National Geographic Foundation.