Mapping Changing Population Distribution on the Qinghai–Tibet Plateau since 2000 with Multi-Temporal Remote Sensing and Point-of-Interest Data

: Advanced developments have been achieved in urban human population estimation, however, there is still a considerable research gap for the mapping of remote rural populations. In this study, based on demographic data at the town-level, multi-temporal high-resolution remote sensing data, and local population-sensitive point-of-interest (POI) data, we tailored a random forest-based dasymetric approach to map population distribution on the Qinghai–Tibet Plateau (QTP) for 2000, 2010, and 2016 with a spatial resolution of 1000 m. We then analyzed the temporal and spatial change of this distribution. The results showed that the QTP has a sparse population distribution overall; in large areas of the northern QTP, the population density is zero, accounting for about 14% of the total area of the QTP. About half of the QTP showed a rapid increase in population density between 2000 and 2016, mainly located in the eastern and southern parts of Qinghai Province and the central-eastern parts of the Tibet Autonomous Region. Regarding the relative importance of variables in explaining population density, the variables “Distance to Temples” is the most important, followed by “Density of Villages” and “Elevation”. Furthermore, our new products exhibited higher accuracy compared with ﬁve recently released gridded population density datasets, namely WorldPop, Gridded Population of the World version 4, and three national gridded population datasets for China. Both the root-mean-square error ( RMSE ) and mean absolute error ( MAE ) for our products were about half of those of the compared products except for WorldPop. This study provides a reference for using ﬁne-scale demographic count and local population-sensitive POIs to model changing population distribution in remote rural areas.


Introduction
The spatial distribution of the human population and its change are critical components in studies focusing on the relationships between humans and the environment [1,2]. The enormous

Study Area
The QTP, consisting of Qinghai Province and the Tibet Autonomous Region, is the largest and highest mountainous region in the world, with an average elevation of over 4000 m above sea level ( Figure 1). The QTP is characterized by low temperature and strong solar radiation and hosts the largest permafrost area in a high-altitude region worldwide. This complex topography and uniquely harsh climate of the QTP play important roles in the uneven human population distribution [30]. The QTP is a typical remote rural region in China, more than two-thirds of the QTP is covered by grassland. Compared with other areas of China, the level of urbanization of the QTP is low. Evidence indicates that humans had occupied this eastern plateau since the Middle Pleistocene epoch [31], however, the north part of the QTP is considered as one of the largest terrestrial wilderness areas in the world [32].
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 18 largest permafrost area in a high-altitude region worldwide. This complex topography and uniquely harsh climate of the QTP play important roles in the uneven human population distribution [30]. The QTP is a typical remote rural region in China, more than two-thirds of the QTP is covered by grassland. Compared with other areas of China, the level of urbanization of the QTP is low. Evidence indicates that humans had occupied this eastern plateau since the Middle Pleistocene epoch [31], however, the north part of the QTP is considered as one of the largest terrestrial wilderness areas in the world [32].

Data Source and Pre-Processing
Demographic data for 2000, 2010, and 2016 and contemporary high-resolution covariate datasets were collected for the estimation of the population density of the QTP. Due to the consideration of previous studies [7,33,34] and the data availability, a wide range of factors known to correlate with population density were selected as predictive covariates, including land use/land cover, topography, vegetation productivity, climate pattern, intensity of nighttime lights, road networks, and the locations of POI data (Table 1). In order to make the data from the contemporary covariate datasets match the ancillary input data used for population density estimation, we resampled all datasets to a spatial resolution of 1000 m using the bilinear algorithm, and all datasets were projected in an Albers equal-area projection.

Data Source and Pre-Processing
Demographic data for 2000, 2010, and 2016 and contemporary high-resolution covariate datasets were collected for the estimation of the population density of the QTP. Due to the consideration of previous studies [7,33,34] and the data availability, a wide range of factors known to correlate with population density were selected as predictive covariates, including land use/land cover, topography, vegetation productivity, climate pattern, intensity of nighttime lights, road networks, and the locations of POI data (Table 1). In order to make the data from the contemporary covariate datasets match the ancillary input data used for population density estimation, we resampled all datasets to a spatial resolution of 1000 m using the bilinear algorithm, and all datasets were projected in an Albers equal-area projection. Table 1. Summary information of the spatial covariate datasets and derived covariates used for the estimation of human population density on the Qinghai-Tibet Plateau (QTP).

Demographic Data and Administrative Unit Boundaries
Demographic counts for 2000 and 2010 on the QTP both at the town-level and county-level were collected from the database of the fifth and sixth census of China, respectively. As the time interval of the national census in China is ten years, the demographic counts in 2016 at the town-level were obtained from the China County Statistical Yearbook (town volume (2017)) in 2017.
The administrative unit boundaries both at the town-level and county-level on the QTP for 2000 were obtained from the Institute of Geographical Sciences and Natural Resources Research [35]. The administrative unit boundary data for 2010 were retrieved from the vector map in 91 Weitu Software. Based on the information on the changes of administrative units in China during 2010-2016 (available at http://www.xzqh.org), the administrative unit boundary data in 2010 were altered to conform with the administrative unit boundaries in 2016. Then, the tabular census data for 2000, 2010, and 2016 were linked to the corresponding administrative units (Figure 2). Due to the lack of town-level administrative unit boundaries for 2000 in 16 counties on the QTP, the missing town-level administrative units were merged and then replaced by county-level administrative boundaries according to previous studies [24,35].  Figure 2). Due to the lack of town-level administrative unit boundaries for 2000 in 16 counties on the QTP, the missing town-level administrative units were merged and then replaced by county-level administrative boundaries according to previous studies [24,35].

Land Use/Land Cover Data
Land use/land cover data for the QTP were obtained from China's land-use/cover datasets (CLUDs). The datasets have been identified as one of the highest accurate land use/land cover products at the national level for China [36]. These datasets have a spatial resolution of 30 m and provide a long time series from the late 1980s to 2015 at five-year intervals [36,37]. This study used the datasets for 2000, 2010, and 2015. The CLUDs are divided into six land use categories with 25 land use sub-categories. The six categories are cultivated land, forest land, grassland, water body, builtup land, and bare land. Among the six categories and 25 sub-categories of land use, we found a small area of mining land on the QTP, mostly in Qinghai Province, which was classified as built-up land. To improve the accuracy of population mapping, according to the sub-categories of land use, this study initially divided the category of built-up land into built-up areas and mining areas. We then converted each category of land use/land cover data into proportion data at a spatial resolution of 1000 m using the aggregate tool in the ArcGIS software. The proportion data were rescaled by 0 to 100% for each category (Figure 3). Compared with the usual binary classifications of land use/land cover at a resolution of 1000 m, the proportion data with continuous values reserves more detailed spatial information.

Land Use/Land Cover Data
Land use/land cover data for the QTP were obtained from China's land-use/cover datasets (CLUDs). The datasets have been identified as one of the highest accurate land use/land cover products at the national level for China [36]. These datasets have a spatial resolution of 30 m and provide a long time series from the late 1980s to 2015 at five-year intervals [36,37]. This study used the datasets for 2000, 2010, and 2015. The CLUDs are divided into six land use categories with 25 land use sub-categories. The six categories are cultivated land, forest land, grassland, water body, built-up land, and bare land. Among the six categories and 25 sub-categories of land use, we found a small area of mining land on the QTP, mostly in Qinghai Province, which was classified as built-up land. To improve the accuracy of population mapping, according to the sub-categories of land use, this study initially divided the category of built-up land into built-up areas and mining areas. We then converted each category of land use/land cover data into proportion data at a spatial resolution of 1000 m using the aggregate tool in the ArcGIS software. The proportion data were rescaled by 0 to 100% for each category (Figure 3). Compared with the usual binary classifications of land use/land cover at a resolution of 1000 m, the proportion data with continuous values reserves more detailed spatial information.
To improve the accuracy of population mapping, according to the sub-categories of land use, this study initially divided the category of built-up land into built-up areas and mining areas. We then converted each category of land use/land cover data into proportion data at a spatial resolution of 1000 m using the aggregate tool in the ArcGIS software. The proportion data were rescaled by 0 to 100% for each category (Figure 3). Compared with the usual binary classifications of land use/land cover at a resolution of 1000 m, the proportion data with continuous values reserves more detailed spatial information.

Geospatial Vector Data
POI data have been increasingly used as ancillary data to improve human population mapping [18][19][20]. The data of rural village POI for 2000 and road networks were obtained from the Institute of Geographical Sciences and Natural Resources Research (IGSNRR). The other POI data including government agencies, tourist destinations, hospitals and clinics, education facilities, and temples, were collected from Baidu Maps (http://map.baidu.com). The POI records from each category were then processed to two types of raster-based derived layers, including a layer with the density of POI and a layer with the Euclidean distance to the nearest geospatial record [7,33,38]. According to existing studies, the layer with the density of POI was generated using the kernel density estimation function in the ArcGIS software [38].

Nighttime Light and Additional Ancillary Data
Nighttime light (NTL) data has been widely used for mapping population density estimation [7,39,40]. The NTL data used in this study included DMSP/OLS NTL data and NPP/VIIRS NTL data. DMSP/OLS data were downloaded from the National Centers for Environmental Information of the National Oceanic and Atmospheric Administration (NOAA) (available at https://ngdc.noaa.gov/eog/dmsp/). The NPP/VIIRS data were downloaded from NOAA/NGDC (available at https://www.ngdc.noaa.gov/eog/viirs/). The DMSP/OLS data were available from 1992 to 2013, while the NPP/VIIRS data were available since 2012. DMSP/OLS data for 2000 and 2010 and NPP/VIIRS data for 2016 were used.
Topographic data were obtained from the United States Geological Survey's global digital elevation model (DEM). Based on the DEM data, the elevation and slope were extracted. The normalized difference vegetation index (NDVI) datasets from MODIS-Terra (MOD13A2) for the period 2000-2016 were obtained from the National Aeronautics and Space Administration (NASA). To reduce contamination by cloud, snow, and ice cover, we smoothed the annual NDVI cycle using an adaptive Savitzky-Golay approach [41,42]. We then calculated the average NDVI for the growing season for the period 2000-2016. Daily precipitation and temperature records for the QTP and its surrounding areas during 2000-2016 were collected at 315 meteorological stations from the China National Stations' Fundamental Elements Datasets V3.0. The gridded datasets of growing-season precipitation and temperature at a resolution of 1000 m were interpolated using the ANUSPLIN version 4.2 software.

Other Gridded Population Distribution Datasets
Five global and national gridded population distribution datasets for 2010 with different spatial resolutions were used in this study. Two global datasets were employed, namely the WorldPop dataset (henceforth referred to as WorldPop) and the Gridded Population of the World v4 dataset (hereafter referred to as GPW v4) [7,12]. There are three national datasets for mainland China [11,43,44].
Main covariate datasets for mapping these three datasets are land cover, nighttime light, and elevation data. These five gridded population distribution datasets were created based on the county-level census counts. It should be noted that the Tibet Autonomous Region was excluded in the gridded population distribution dataset from Gaughan et al. (2016, CnPop3) [11]. The aggregate tool in the ArcGIS software was used to upscale the high spatial resolution (e.g.,~100 m) in the data to 1000 m. Summary information of the five population gridded datasets is provided in Table 2. 2.3. Methods

Random Forest-Based Dasymetric Population Mapping Approach
The RF-based dasymetric population mapping approach is a combination of an RF regression model and a dasymetric population mapping method, and has been successfully used as an underlying model for the estimation of human population density both at national and global scales, due to its excellent performance and generalization capabilities [7,11,38]. In brief, the major processes of this hybrid population mapping approach can be divided into two steps: (1) an RF regression model is initially used to generate a population density weighting layer, which (2) is subsequently applied to dasymetrically redistribute the original census counts within administrative units [7,11].
The general process for RF-based dasymetric population mapping is described in Stevens et al. (2015) [7] and Gaughan et al. (2016) [11]. Several studies have noted that the higher accuracy of predicted population values at the administrative-unit level leads to a higher accuracy of the gridded datasets of population density [14,45]. The town-level (Jiedao/Xiangzhen in Chinese, equivalent to fourth-level administrative units) is the finest administrative unit for which demographic data are publicly released in China; however, it is very difficult to assess the accuracy of the gridded datasets using finer-scale administrative units for census counts. Therefore, in this study, the hybrid population mapping approach was modified to map multi-temporal high-resolution human population density in the QTP. Four steps are outlined in Figure 4. These steps are as follows: (i) fitting the RF model to generate a population density weighting layer using town-level census counts and their corresponding predictive covariates ( Figure 5); (ii) dasymetric redistributing of the county-level census counts (see Figure 2a) proportional to the predicted population density weighting layer at the grid-cell level (namely, county-based QTPpop); (iii) comparing the accuracy of the county-based QTPpop with currently available gridded population datasets; and (iv) when the assessment shows that the county-based QTPpop is robust, dasymetric redistributing the town-level census counts (see Figure 2b) proportional to the predicted population density weighting layer at the grid-cell level (namely, town-based QTPpop). According to previous studies [11,29,46], the RF model was used for each census year (2000, 2010, and 2016) to generate multi-temporal gridded datasets of the population density on the QTP, which were subsequently used to detect the changes in population density. The RF model was fitted in the R statistical software using the randomForest package. More details on the parameterization of the RF model can be found in Stevens et al. (2015) [7].

Technical Validation and Inter-Comparison
(1) Out-of-bag (OOB) error estimation The out-of-bag (OOB) error estimation, which is an internal cross-validation component of the RF model, provides an accurate and unbiased estimate of the overall accuracy of the RF model [47,48]. During the RF model fitting, one-third of the input data are held in reserve from the bagging iteration process and utilized to obtain an OOB error. The prediction error of the RF model for population density mapping was estimated by averaging the mean squared errors (MSE). Moreover, the OOB error was utilized to evaluate the importance of each model covariate by randomly permutating a given covariate's OOB data with random noise and computing the percentage increase in the mean squared error (%IncMSE) across all trees of the RF model. Higher values of %IncMSE represent greater variable importance [11].
(2) External Accuracy Assessment The mapping accuracy is compatible among population density datasets within a given year, nevertheless, the comparisons cannot be reliably performed for datasets from different years [49].

Technical Validation and Inter-Comparison
(1) Out-of-bag (OOB) error estimation The out-of-bag (OOB) error estimation, which is an internal cross-validation component of the RF model, provides an accurate and unbiased estimate of the overall accuracy of the RF model [47,48]. During the RF model fitting, one-third of the input data are held in reserve from the bagging iteration process and utilized to obtain an OOB error. The prediction error of the RF model for population density mapping was estimated by averaging the mean squared errors (MSE). Moreover, the OOB error was utilized to evaluate the importance of each model covariate by randomly permutating a given covariate's OOB data with random noise and computing the percentage increase in the mean squared error (%IncMSE) across all trees of the RF model. Higher values of %IncMSE represent greater variable importance [11].
(2) External Accuracy Assessment The mapping accuracy is compatible among population density datasets within a given year, nevertheless, the comparisons cannot be reliably performed for datasets from different years [49].

Technical Validation and Inter-Comparison
(1) Out-of-bag (OOB) error estimation The out-of-bag (OOB) error estimation, which is an internal cross-validation component of the RF model, provides an accurate and unbiased estimate of the overall accuracy of the RF model [47,48]. During the RF model fitting, one-third of the input data are held in reserve from the bagging iteration process and utilized to obtain an OOB error. The prediction error of the RF model for population density mapping was estimated by averaging the mean squared errors (MSE). Moreover, the OOB error was utilized to evaluate the importance of each model covariate by randomly permutating a given covariate's OOB data with random noise and computing the percentage increase in the mean squared error (%IncMSE) across all trees of the RF model. Higher values of %IncMSE represent greater variable importance [11].
(2) External Accuracy Assessment The mapping accuracy is compatible among population density datasets within a given year, nevertheless, the comparisons cannot be reliably performed for datasets from different years [49]. Thus, this study assessed the accuracy of our population density product on the QTP for 2010 by comparing it with five other gridded population datasets (WorldPop, GPW v4, CnPop1, CnPop2, CnPop3) for the same year. The mapping accuracy was assessed via the following steps. First, the county-based QTPpop for 2010 was compared to the town-level census counts by summing gridded population estimates within each town administrative unit. Then, scatter plots of the town-level census count and the town-level estimated population were constructed. The coefficient of determination (R 2 ), root-mean-square error (RMSE), and mean absolute error (MAE) were selected as the statistical measures of mapping accuracy. The same steps were used to assess the accuracy of the five gridded population products.

Accuracy Assessment and Inter-Comparison
The RF model used to predict human population density in the QTP has a consistently high The county-based QTPpop for 2010 was highly consistent with the original census counts at the town-level, with an R 2 of 0.97 ( Figure 6). This value of R 2 was higher than those of the five publicly released gridded population density datasets; the R 2 value for Worldpop was 0.90, while those of the other four datasets were lower than 0.90, with the lowest value being obtained for GPW v4 (R 2 = 0.40). Spatially, the county-based QTPpop better mirrored the characteristics of the population pattern compared to the five publicly available products ( Figure 6). First, compared with Cnpop1, CnPop3, and GPW v4, in county-based QTPpop, the phenomenon in which the population density maps displayed abrupt changes along the administrative unit boundaries was mitigated. Second, compared with Cnpop2, in county-based QTPpop, according to the original census counts at the town-level (Figure 2b), the large underestimation of population density in the areas with low population density in the northwestern QTP was mitigated. The county-based QTPpop was more consistent with the WorldPop product than the four other publicly available datasets in terms of the spatial pattern of population distribution on the QTP.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 18 Thus, this study assessed the accuracy of our population density product on the QTP for 2010 by comparing it with five other gridded population datasets (WorldPop, GPW v4, CnPop1, CnPop2, CnPop3) for the same year. The mapping accuracy was assessed via the following steps. First, the county-based QTPpop for 2010 was compared to the town-level census counts by summing gridded population estimates within each town administrative unit. Then, scatter plots of the town-level census count and the town-level estimated population were constructed. The coefficient of determination (R 2 ), root-mean-square error (RMSE), and mean absolute error (MAE) were selected as the statistical measures of mapping accuracy. The same steps were used to assess the accuracy of the five gridded population products.

Accuracy Assessment and Inter-Comparison
The RF model used to predict human population density in the QTP has a consistently high The county-based QTPpop for 2010 was highly consistent with the original census counts at the town-level, with an R 2 of 0.97 ( Figure 6). This value of R 2 was higher than those of the five publicly released gridded population density datasets; the R 2 value for Worldpop was 0.90, while those of the other four datasets were lower than 0.90, with the lowest value being obtained for GPW v4 (R 2 = 0.40). Spatially, the county-based QTPpop better mirrored the characteristics of the population pattern compared to the five publicly available products ( Figure 6). First, compared with Cnpop1, CnPop3, and GPW v4, in county-based QTPpop, the phenomenon in which the population density maps displayed abrupt changes along the administrative unit boundaries was mitigated. Second, compared with Cnpop2, in county-based QTPpop, according to the original census counts at the town-level (Figure 2b), the large underestimation of population density in the areas with low population density in the northwestern QTP was mitigated. The county-based QTPpop was more consistent with the WorldPop product than the four other publicly available datasets in terms of the spatial pattern of population distribution on the QTP. Furthermore, the RMSE and MAE of the county-based QTPpop were 4027.23 and 2024.67, respectively, both of which were lower than those of the five publicly available gridded population datasets (Figure 7). Specifically, the RMSE of the county-based QTPpop was about half that of Cnpop1, Cnpop2, WorldPop, and CnPop3, and about one-fifth that of GPW v4. The MAE of the county-based QTPpop was approximately half that of the five publicly available gridded population products except for WorldPop. Internal measures of the accuracy of the RF model and an inter-comparison of external accuracy assessments among the five recently released gridded population datasets indicated that the county-based QTPpop yielded a higher mapping accuracy compared to the five other products. One important explanation is that the use of fine-scale administrative units for census counts improved the model's ability for population mapping. As above-mentioned, the accuracy of the population density dataset based on county-level census counts was higher than those of the five publicly available population density products. Meanwhile, the town-level was the finest statistical unit in publicly released demographic data for China, and the town-based QTPpop took full advantage of the spatially detailed town population counts. It is reasonable to expect that the town-based QTPpop can yield higher accuracy than the county-based QTPpop. More interpretation is given in Section 4.1. Thus, we used the town-based QTPpop in subsequent analyses.

Importance of Variables in Explaining Population Density
The relative importance of variables in explaining population density on the QTP is shown in Figure 8. Of the 24 variables used in the RF model, the 10 that contributed most to the estimation of population density can be subdivided into three categories: (1) six POI-related variables; (2) two land-use related variables; and (3) variables related to topographical features and vegetation productivity. The variable "Distance to Temples" has the largest contribution to the prediction of population density. This variable is closely related to the religious beliefs of inhabitants across the QTP. Aside from this variable, the two variables that contributed the most to the prediction of population density were "Density of Villages" and "Elevation". The fact that elevation is among the top three variables is understandable since most inhabitants on the QTP live in areas with relatively low elevation such as valleys. The "Density of Villages" and "Distance to Villages" variables were good measures of population density in the QTP. Meanwhile, the importance of the variable "Percentage of Cultivated Areas" was greater than that of "Percentage of Built-up Areas", which may be related to the low level of urbanization across the QTP and the fact that the population density of farmers depends on the distribution of cultivated areas. Surprisingly, the importance of the "Night-lights' intensity" variable was outside the top 10, likely due to the low level of urbanization in the study area.

Spatio-Temporal Change of Human Population Density on the QTP
The town-based QTPpop provided an accurate description of the population density on the QTP in 2000, 2010, and 2016 ( Figure 9). Overall, the QTP has a sparse population distribution. In about 14% of the QTP, the population density is zero; these areas are mainly distributed around the Kunlun Mountains and Altun Mountains, including the north of the Tibet Autonomous Region and western Kekexili, Qinghai. Areas with relatively high population density are primarily located in the Yellow River-Huangshui River Valley and the Yarlung Zangbo and Lhasa River Valley.
From 2000 to 2016, the extreme sparsely populated areas in the QTP gradually decreased while the relatively densely populated area rapidly increased (Figure 9). The proportion of areas with a population density of less than 0.5 persons/km 2 decreased from 41.62% in 2000 to 40.23% in 2010 to 36.52% in 2016. In contrast, the proportion of areas with a population density between 5 and 10 persons/km 2 increased from 4.82% in 2000 to 6.25% in 2010 to 6.70% in 2016, while the proportion of areas with a population density of above 10 persons/km 2 increased from 4.26% in 2000 to 4.91% in 2010 to 5.35% in 2016.
The areas in which the population growth between 2000 and 2016 exceeded 0.1 persons/km 2 , accounting for 42.81% of the total area of the QTP, mainly distributed in the eastern and southern parts of Qinghai Province and the central-eastern parts of the Tibet Autonomous Region (Figure 10c). The proportion of areas with population growth between 2 and 5 persons/km 2 and above 5 persons/km 2 were 5.14% and 2.57%, respectively. In contrast, the areas with a population decrease above 0.1 persons/km 2 accounted for 15.71% of the total area, and were mainly distributed in agricultural and pastoral areas around relatively densely populated areas. Additionally, in both the periods of 2000-2010 and 2010-2016, the areas with increased population density tended to be areas with higher population density, while the areas with unchanged or lower population density were mainly located in more sparsely populated areas (Figure 10a,b).    (Figure 10c). The proportion of areas with population growth between 2 and 5 persons/km 2 and above 5 persons/km 2 were 5.14% and 2.57%, respectively. In contrast, the areas with a population decrease above 0.1 persons/km 2 accounted for 15.71% of the total area, and were mainly distributed in agricultural and pastoral areas around relatively densely populated areas. Additionally, in both the periods of 2000-2010 and 2010-2016, the areas with increased population density tended to be areas with higher population density, while the areas with unchanged or lower population density were mainly located in more sparsely populated areas (Figure 10a,b).

Advancement Using Town-Level Population Counts in Population Mapping for the QTP
In this study, a detailed multi-temporal human population distribution dataset was generated for the QTP based on town-level population counts. A comparison with five existing national and global gridded population products for 2010 indicated that our new population products more accurately characterized the human population distribution in the QTP (Figures 6 and 7). In this work, the higher accuracy of the population density weighting layer was initially predicted using the RF model based on the town-level population counts. Previous studies have noted that the choice of finer census units generally leads to a larger range of measured population densities and better reflects geographic variation in population, which can minimize the span of the model's spatial-scale conversion [12,25,45]. Meanwhile, the dasymetric redistribution of lower population counts in finer administrative units using the predicted population density weighting layer can reduce the error in the final human population density datasets [45]. Therefore, compared with the county-based QTPpop, the town-based QTPpop had a higher accuracy.
The size of administrative units at the same hierarchy of China exhibited considerable differences over the QTP [16], with the size of administrative units of lower hierarchy being much larger than that of higher hierarchy (Figure 2). This phenomenon may introduce uncertainties in human population mapping when the population counts with the higher level of administrative units are used. In the following, Golmud County (located in western Qinghai Province) is taken as an example. This county consists of two unconnected parts, namely Tanggulashan Town and other areas. Tanggulashan Town has a low population density, while other areas in Golmud County have a higher population density ( Figure 2). Some of the currently released population products analyzed in the present stud greatly overestimate the population in Tanggulashan Town, which were mapped based on the county-level census (Figure 6c,g). Similar phenomena also occurs in the northwestern QTP ( Figure 6).

Effectiveness of Local Population-Sensitive Points-of-Interest (POIs) for Improving Population Mapping in Remote Rural Areas
Detailed covariates have been shown to improve the accuracy of models of population density [7]. This study used population-sensitive POI data as input variables in the RF regression model, which improved the accuracy of the estimation of the human population density in the QTP. It was found that of the 10 variables that were most important for explaining population density, six were POI-related variables such as "Distance to Temples", "Density of Villages", and "Distance to Villages". This result indicates that, for the mapping of the population density of the QTP, it is necessary to introduce local social and cultural variables as well as physical variables (e.g., topographical features) [35,50]. In contrast, although the intensity of nighttime lights has been considered as an excellent predictor of urban population density [16,17], the importance of the variable "Night-lights' intensity" may significantly weaken in low-and middle-income regions, especially in remote rural areas such as the QTP.
The town-based QTPpop product that was constructed in this study consisted of temporally comparable population datasets based on demographic data for 2000, 2010, and 2016 and a number of covariates that are time-invariant and temporally explicit. Compared with national population datasets (e.g., Xu et al. (2017) [43], Gaughan et al. (2016) [11], and Tan et al. (2018) [16]), we not only used variables related to built-up areas and night-light intensity to reconstruct the changes in population distribution over time, but also employed population-sensitive POIs variables (e.g., distances to temples and settlements) as predictors of population distribution in remote rural areas.

Population Growth on the QTP
The results of this study suggest that between 2000 and 2016, the QTP experienced a reduction in areas with extremely low population density. This is contrary to the results of a study that found that, in mainland China, between 2000 and 2016, there was an overall expansion of areas with low population density [16]. One potential explanation for this is that the level of urbanization on the QTP is much lower than that in Central and Eastern China. The results of the present study indicate that in rural areas of the QTP, about half of the total area showed an increase in population density, which may increase human pressures on this world's largest pastoral alpine ecosystem [51]. Up to now, more than 150 nature reserves have been created on the QTP, which account for one-third of its total area [52]. Furthermore, since 2004, an ecological restoration project called the Grazing Withdrawal Program (GWP) has been implemented to restore the grassland ecosystems of the QTP. The widespread population growth on the QTP may offset the effectiveness of its nature reserves and the GWP. Therefore, more attention should be paid to coordinating population growth and the preservation of the fragile ecological environment of the QTP, especially in nature reserves.

Limitations and the Future Research
This work provides a case study of population density mapping in remote rural areas based on a RF model. Even though the accuracy of our newly products were higher than five publicly released national or global gridded population products, recent studies noted that deep learning could achieve a higher accuracy of population mapping than the shallow machine learning [53]. Meanwhile, with the improvement in the geospatial big data of POIs and the location-based social media data, more attention should be paid to improving the accuracy of human population mapping in remote rural areas in future work [1], especially for temporally comparable population datasets.

Conclusions
In this study, the temporal and spatial change of human population density on the QTP from 2000 to 2016 was estimated based on demographic data at the town-level, remotely sensed data, and local population-sensitive POI data. The results indicated that our new population map exhibited higher accuracy compared with five recently released gridded population density datasets. Both the root-mean-square error and the mean absolute error for our map were about half of those of the compared products (with the exception of the WorldPop). The higher accuracy of our new population map largely depended on taking advantage of the finer scale demographic count and local population-sensitive POIs. POI-related variables (e.g., "Distance to Temples" and "Density of Villages") were more important for explaining population density on the QTP than the variable "Night-lights' intensity". Furthermore, in large areas of the northern QTP, accounting for about 14% of the total area of the QTP, population density was zero. About half of the QTP showed a rapid increase in population density from 2000 to 2016, mainly located in the eastern and southern parts of Qinghai Province and the central-eastern parts of the Tibet Autonomous Region.