Night light data were obtained from the China Long Time-Series Night Light Dataset (2000–2020) from the Global Change Science Research Data Publishing System (http://www.geodoi.ac.cn). The spatial resolution of the EANTLI-like night-light dataset from 2000 to 2020 was 1 km.

NDVI data were obtained from the MOD13Q1 product data from the National Aeronautics and Space Administration (NASA) (https://ladsweb.modaps.eosdis.nasa.gov). The resolution was 250 m and the original time resolution was 16 d. Therefore, comprehensive NDVI data covering 2000–2020 were obtained using the maximum synthesis method. Python and MRT were used to preprocess the data after download. The ArcGIS platform was then used for data calculation and creating charts.

The meteorological data used in this study were obtained from the monthly average temperature (Peng, 2019) and precipitation (Peng, 2020) dataset from the National Tibetan Plateau Scientific Data Center (https://data.tpdc.ac.cn/). The selected time spanned 2001–2020 and the spatial resolution was 1 km. This dataset is based on the CRU global 0.5° climate dataset issued by the National Center for Atmospheric Sciences, United Kingdom, and the global high-resolution climate dataset published by World Climate. This dataset was generated in China using a delta spatial downscaling scheme. In addition, 496 independent meteorological observation data points were used for verification, with reliable results (Peng et al., 2017; Peng et al. 2018; Peng et al. 2019; Ding and Peng. 2020). The annual temperature was calculated as the average of the monthly temperature; precipitation was calculated as the monthly cumulative value.

Digital Elevation Model data were obtained from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http://www.resdc.cn). This dataset was generated by collating and splicing provisional data from the most recent SRTM V4.1 data. The resolution was 90 m, the data type was img, and the projection was UTM/WGS84.

The ArcGIS platform was used for projection conversion, clipping, and operations to obtain the data for this study. In addition, the “Slope” function of the “Spatial Analyst” module in the ArcGIS platform was used to process and obtain slope and other data.

The total population (10,000 people) and urban population (10,000 people) were derived from statistical data from the Guizhou Provincial Macroeconomic Database, which can be found on the main board of the Guizhou Provincial Bureau of Statistics (http://hgk.guizhou.gov.cn). In addition, the GDP (100 million yuan), built-up area (km²), number of buses (vehicles), total number of passengers (10,000 persons), and taxis (vehicles) were obtained from the China Statistical Yearbook on the website for the National Bureau of Statistics (http://www.stats.gov.cn/tjsj./ndsj). The data were collected from 2000 to 2020.

Land use and cover change (LUCC) data were obtained from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http://www.resdc.cn). The spatial resolution was 1 km. These data were classified using the LUCC classification system, which is widely used in China. The land use types were divided into six first-class land types: cultivated land, forestland, grassland, water area, construction land, and unused land. The land use data for 2000, 2010, and 2020 were first selected for this study. Then, the scope of this study was determined based on the administrative boundary of the study area. Finally, the ArcGIS platform was used to extract the corresponding research data for analysis and calculation. Spatial data were cropped and calculated on the ArcMap 10.2 platform using the administrative boundaries of the study area.

This study used nighttime light data to reflect human activities. First, the nighttime light index was reduced based on the administrative boundary of the study area. The image was then normalized using the pixel-dichotomy model, calculated as follows: Y i = L i − L min L max − L min , (1) where Y _i represents the nighttime light index (dimensionless) after normalization of the ith pixel, L _i represents the nighttime light index of the ith pixel, L _min represents the minimum nighttime light index over the study period in the study area, and L _max represents the maximum nighttime light index in the study area.

A GeoDetector is a statistical method for analyzing the spatial heterogeneity of factors and explaining the underlying factors. It comprises a risk detector, factor detector, ecological detector, and interaction detector that can be used to quantitatively detect and identify various interactions between spatial distribution attributes and their explanatory factors (Wang and Xu. 2017). It was used to analyze the explanatory power (i.e., q value) of each driving factor on the distribution of a geographical object using quantitative detection of spatial heterogeneity. The formula for the q value is as follows: q = 1 − 1 n σ 2 ∑ h = 1 L n k σ k 2 , (2) where n k represents the number of samples of type h of factor A, n is the total number of samples in the entire region, σ k 2 represents the discrete variance of the h type of factor A, and σ 2 represents the discrete variance of the entire region. The range of q is [0,1]. The greater the value of q, the greater the impact of the factor on the spatial distribution of the landscape.

The least squares method was used in this study to analyze the pixel-by-pixel change rule of images (i.e., annual change rate) based on the MATLAB platform. The calculation formula is as follows, taking the NDVI as an example: θ s l o p e = n × ∑ i = 1 n i × N D V I i − ∑ i = 1 n i × ∑ i = 1 n N D V I i / n × ∑ i = 1 n i 2 − ∑ i = 1 n i 2 , (3) where i represents the serial number of the year, N D V I i represents the NDVI value of the ith year, N is the sample size, and θ s l o p e represents the rate of change in the NDVI during the study period (i.e., the slope of the linear regression equation representing the least squares of interannual change).

The F-test was used to determine the significance of the changing trend, calculated as follows: F = S R S E / n − 2 , (4) where S R represents the sum of the error squares, S E represents the sum of the residual squares, and (n—2) represents the residual degree of freedom.

The correlative analysis method was used to measure the closeness of two or more variables by analyzing them with correlation, calculated as follows: r x y = ∑ i = 1 n x i − x ¯ y i − y ¯ / ∑ i = 1 n x i − x ¯ 2 × ∑ i = 1 n y i − y ¯ 2 , (5) where r x y represents the correlation coefficient between the elements x and y and x ¯ and y ¯ represent the average values of the sample values of the two elements. The range of correlation coefficient r is [− 1, 1]: when r is positive, it is positive; when the r value is negative, it is negative. The significance test was performed after determining the significance degree using the critical value table of the test correlation coefficient.

The land use transfer matrix depicts the land use/cover structure and its type and direction of change. This reflects the dynamic processes that drive land use changes over time. The general formula for the land-use transfer matrix is as follows: s i j = s 11 ⋯ s 1 n ⋮ ⋱ ⋮ s n 1 ⋯ s n n , (6) where S _i is the land-use type area before the transfer, S _j is the land-use type area after the transfer, i = 1,2,3 … n, j = 1,2,3 … n, and n is the number of land-use types.

The global Moran index can be sued to reveal the aggregation of the spatial factor layout as a whole (Zhao et al., 2020), calculated as follows: M o r a n ′ s   I = n ∑ i = 1 n ∑ j = 1 n W i j x i − x ¯ x j − x ¯ ∑ i = 1 n ∑ j = 1 n W i j ∑ i = 1 n x i − x ¯ 2 , (7) where n represents the number of samples, x i and x j represent the observations of attribute x at spatial locations i and j, x ¯ represents the average value of attribute x, and W i j represents the spatial weight matrix. This study adopted an adjacency weight matrix. If the ith region was adjacent to the jth region, W i j was one; otherwise, it was zero. Moran’s I had a range of [− 1,1]. An index > 0 indicates positive spatial autocorrelation (aggregation mode), whereas an index < 0 indicates a negative spatial correlation (dispersion mode).

The local Moran index, LISA, is used for local spatial autocorrelation analysis to reveal the high-high and low-low spatial concentration areas of human activities and to understand the local spatial distribution characteristics: L o c a l   M o r a n ′ s   I = n x − x ¯ ∑ j = 1 n W i j x j − x ¯ ∑ i = 1 n x i − x ¯ 2 , (8) where the definition of each parameter in Eq. 8 is the same as that in Eq. 7.

Spatial correlation was calculated using the ArcMap 10.2 platform. The spatial rate of change and significance test were performed using the MATLAB platform. Linear correlations were analyzed using SPSS 20.0. Spatial autocorrelation analysis was performed using GeoDa.

The analysis revealed that the range and intensity of human activities extracted based on the NLI were positively correlated with the selected indicators; in particular, the urban population, GDP, built-up area (km²), number of buses (vehicles), total number of passengers (10,000 person times), and number of taxis (vehicles) were significantly positively correlated at the 0.05 level (Table 1). In summary, the NLI-derived human activity intensity and range could effectively reflect relevant human activities.

The degree of influence of the factors related to human activities (HA) can be analyzed from space using the GeoDetector model. The analysis results showed that the q value of construction land to the spatial heterogeneity of human activity area calculated based on the NLI was 0.8347 according to land-use type data from the Chinese Academy of Sciences classification system.

Spatial analysis of the extent of the NLI revealed that in 2000, 2010, and 2020, 77.58, 84.98, and 93.76% of the built-up land area, respectively, coincided with the NLI area. However, the range of the nighttime lighting index was significantly larger than that of construction land.

This was significantly higher than that of other land-use types (note: these findings were obtained by analyzing land-use data from the Resources and Environment Center of the Chinese Academy of Sciences). These results showed that the spatial distribution of human activities reflected by the NLI could accurately reflect the spatial status of human activities.

Human activities accounted for 57.78% of the total study area from 2000 to 2020, with a significant annually increasing trend (y = 4,279.2x—2080.5, R² = 0.8564), in particular from 2012 to 2013 (Figures 1K–O). Human activities were more extensive and intensive in the central and southwestern regions than in the eastern region. Human activities in Guizhou Province were relatively intense during this period in regions where the prefectures and municipal governments were located, in particular in the central and southwest regions, such as Guiyang, Zunyi, Anshun, Liupanshui, and Xingyi (Figure 1E). Regarding time, the intensity of human activities increased significantly (y = 0.0006x + 0.0068, R² = 0.9043), particularly in 2012 (Figures 1L–O). Regarding space, 97.56% of the human activity area showed an increasing trend in human activity, 86.05% of the total area of human activity showed a significant increasing trend, 2.44% of the regions showed a weakening trend, and 0.84% of the regions showed a significant weakening trend. The intensity and area of human activities in Guizhou Province increased from 2000 to 2020. Human activities in Guiyang City were the most visible. In addition, human activities increased significantly at the government stations of Zunyi, Anshun, Liupanshui, Xingyi, Tongren, Qiannan, Qiandongnan, Bijie, and other prefectures and cities, and in other county areas. However, there was a notable decrease in human activity in the southwestern part of the study area (Figure 1G).

Human activity accounted for 15.67% of the total study area from 2000 to 2005. The regions with increased human activity accounted for 66.34% of the total human activity area, and 35.42% of the regions showed a significant increase. Reduced human activity areas accounted for 33.66% of the total human activity area, with 12.53% showing a significant reduction. The scope of human activities was limited at this stage of the study. However, human activities increased significantly in Guiyang, Zunyi, Anshun, Liupanshui, Xingyi, and other regions. The increase in human activities was most visible in the trend of cities as the focal point, with a gradual weakening (Figures 1A,F).

From 2006 to 2010, the area of human activities accounted for 18.39% of the total study area. The overall performance in this stage was an increase in human activity, which continued to increase in intensity and scope, with the strongest intensity and scope occurring in the city. However, the southwestern region of the study area (e.g., Weining County) showed a trend of weakening human activity and a narrowing scope. Human activities in the southeastern region of the study area showed a notable trend of belt enhancement due to the construction of roads and other infrastructure (Figures 1B,G).

Human activity accounted for 44.60% of the total study area from 2011 to 2015. The area accounting for 95.49% of the total human activity area showed an increasing trend, with a significant increase in the area accounting for 67.16%. Human activities decreased in the area accounting for 4.51% of the total area of human activities, with 1.10% showing a significant decrease. The area of human activity at this stage was significantly larger than that in the previous stage. Human activities in some counties weakened due to policy adjustments and other reasons. However, overall human activities were expanding in scope and intensity. Due to road construction, human activities were more clearly distributed in strips in the southeast and northeast (Figures 1C,H).

Human activities accounted for 55.03% of the total study area from 2016 to 2020. The area accounting for 66.83% of the total human activity area increased, with 30.25% showing a significant increase. The area accounting for 33.17% of the human activity decreased, with 10.53% showing a significant decrease. Compared to the previous stage, the scope of human activities continued to increase along with the overall intensity of human activities. At this stage, the intensity of human activities in the north and east of the study area increased significantly. In contrast, some urban areas in the middle and southeast of the study area weakened due to policy changes and other human activities. In contrast, human activities in rural residential areas increased significantly, particularly in the southeast and north of the study area. The banding development of human activities caused by infrastructure construction, such as roads, remained visible, in particular in the south and east of the study area (Figures 1D,I).

We analyzed the shifts in human activities, NDVI, temperature, and precipitation gravity centers in four stages in Guizhou Province from 2000 to 2005, 2006–2010, 2011–2015, and 2016–2020, on a 5-year interval basis (Figure 2).

During the study period, the concentrated area of human activity intensity moved 20.71 km to the northeast. According to the specific changes in each stage, the focus of human activities moved 5.8 km southward from 2006 to 2010 compared to 2000–2005; from 2011 to 2015, the focus of human activity intensity moved 8.92 km to the northeast compared to 2006–2010; and from 2016 to 2020, it continued to migrate to the northeast by 15.37 km compared to 2011–2015. The distribution direction of human activities in Guizhou Province did not change, indicating a northeast–southwest direction according to the SDE. The spatial distribution direction of human activities has become increasingly notable over time (Figure 2A).

Temperature and precipitation were relatively stable among the meteorological factors; the SDE and focus of each research period were consistent. The NDVI was also relatively stable, with small variations in the SDE and focus migration. The SDE direction and focus position of the temperature, precipitation, and NDVI were consistent (Figures 2B–D). In summary, except for significant spatial changes in human activities during the study period, the NDVI and temperature were relatively stable and precipitation showed minor spatial changes.

The local Moran’s index (LMI) and local indicators of spatial association (LISA) were used to show the spatial concentration distribution of human activities (Figures 2A–E) and its significance (Figures 3F,G). The global Moran’s index (GMI) was used to reveal the aggregated spatial factor layout (Figures 2K–O). An increase or decrease in the GMI indicated that the factor space was contracting or expanding, respectively.

The GMI was 0.710 from 2000–2005 (Figure 3K), indicating a high degree of aggregation. According to the LMI, human activities were mainly concentrated in the middle and southwest of Guizhou (Figures 3A,F).

The GMI was 0.731 from 2006 to 2010 (Figure 3L), and the spatial concentration of human activities increased. According to the LMI, human activities were mainly concentrated in the middle and southwest of Guizhou, with an increasing trend in the northeast of Guizhou (Figures 3B,G).

The GMI was 0.657 from 2011 to 2015 (Figure 3M), indicating that the aggregation of human activities in this stage was lower than that from 2006 to 2010. However, the scope of human activities expanded significantly during this stage. According to the LMI, high-value areas of human activity were still concentrated in the middle and southwest of Guizhou while low-value areas were in the southeast and north of Guizhou at this stage, indicating that human activities in this area began to increase. However, compared to the original high-value areas of HA, HA in this area remained low (Figures 3C,H).

The GMI was 0.634 from 2016 to 2020 (Figure 3N), which was more extensive than the range of human activities in the previous stage. According to the LMI, there was a greater concentration of low-value areas in the southwest of Guizhou, with the appearance of new high-value areas (Figures 3D,I).

In summary, Moran’s index was 0.670 (Figure 3O), indicating that human activity had relatively notable aggregation from 2000 to 2020. According to the LMI, the human activity intensity and level of urban development were high in the middle, southwest, and northeast of Guizhou. Conversely, the HAI and level of urban development were relatively low in the north, west, and southeast of Guizhou. The reduction in the aggregation of HA, as indicated by the spatial change in the LMI, also indicated that the scope of HA in Guizhou expanded during the study period (Figures 3E,J).

In terms of altitude, Guizhou had low altitudes in the southeast and east, higher altitudes in the middle and west, and the highest altitudes in the southwest (Figure 4B). In terms of slope, 0°–5° slopes accounted for 27.33% of Guizhou, 6°–15° slopes accounted for 48.08%, 16°–25° slopes accounted for 19.34%, 26°–35° slopes accounted for 4.53%, 36°–45° slopes accounted for 0.65%, and > 46° slopes accounted for only 0.06%. The terrain of Guizhou was relatively smooth in the middle and rugged in the west, southwest, and south (Figure 4C).

The annual average temperature in Guizhou was > 7°C. The regions with higher average temperatures were in the east, southeast, and northeast. The regions with lower average temperatures were in the middle. Those with the lowest temperatures were in the west (Figure 4D). Annual precipitation in the study area exceeded 800 mm, with more precipitation in the east than in the west, more in the south than in the north, and gradually decreasing from the southeast to the northwest (Figure 4E). Low temperatures and less precipitation occurred in the middle and west; whereas, higher temperatures and more precipitation occurred in the east, southeast, and northeast.

The spatial distribution effects of elevation, slope, temperature, and precipitation on the spatial heterogeneity of HA were analyzed using the GeoDetectors model. The factors influencing the spatial heterogeneity of human activity intensity in Guizhou Province were the slope, precipitation, altitude, and temperature, ranging from high to low based on the degree of impact. Furthermore, the factor works better in combination with other factors than a single factor (Figure 4F).

Terrain factors significantly impacted human activities in Guizhou Province, with slope factors having a more significant impact than altitude factors. Therefore, human activity was concentrated in low-lying areas. The analysis of concentrated human activity areas (normalized value > 0.1) revealed that human activity intensity areas were mainly distributed in areas with a gradient of 1°–21° (based on 95% of the areas with a normal distribution). Therefore, human activities were concentrated primarily in areas with the gentlest slopes (6°–15°), accounting for 51.47% of the total area of HA. The area with flat slopes (0°–5°) accounted for 39.78%, and the area with 16°–25° slopes accounted for 8.76%. Using the same method, we found that human activities were primarily distributed in areas with an altitude from 489 to 1,982 m, temperatures from 12.36°C to 17.74°C, and precipitation from 1,001.99 to 1,276.99 mm.

In summary, HA were mainly distributed in areas with gentle slopes, higher precipitation, lower altitudes, and higher temperatures in the study area. The factor works better in combination with other factors than a single factor. Slope and precipitation had the greatest impacts.

The ranges of significant positive correlations between precipitation and human activity were slightly higher than that of temperature based on spatial distribution. However, overall, precipitation and temperature significantly impacted human activities, particularly the combination of water and heat (Figures 5A,B; Table 2).

The NDVI value was lower in regions with high HAI and higher in areas with low HAI. The mean NDVI of the study area was 0.823 and that of the human activity area was 0.812. Low NDVI values were found in all areas with high human activity, particularly in urban areas (Figure 6A).

In both Guizhou Province and the human activity area, the rate of change in the NDVI showed a large area occupied by an increase and significant increase. The proportion of NDVI reduction and a significant reduction trend in the human activity area was slightly higher than in Guizhou (Figure 6D; Table 3).

Based on the driving force analysis of the NDVI in Guizhou using the GeoDetector model, we found that human activities were the most important influencing factor of the NDVI distribution in Guizhou, followed by slope, precipitation, altitude, and temperature, which had small effects. Moreover, the effect of each factor acting in combination with other factors was greater than the effect of each factor acting in isolation. The combined effect of human activities and precipitation had the most significant impact on the spatial heterogeneity of the NDVI (Figure 6F).

Regarding the correlation between meteorological factors and NDVI, the difference in the proportion of various correlations between the NVDI and meteorological factors in anthropogenic regions and Guizhou was small. The proportion of areas showing a significant positive correlation between human activities and NDVI was significantly larger than that of meteorological factors. However, we note that the proportion of areas showing a significant negative correlation between human activities and NDVI was also greater than that of the meteorological factors (Figures 6B,C,E; Table 4).

Human activities had a significant influence on changes in the NDVI. Human activities and NDVI in human activity areas mainly showed positive and significant positive correlations. Furthermore, the significantly negative correlation between human activities and NDVI was greater than that between natural factors and NDVI. HA generally had a significant impact on NDVI. Human activities were mostly negatively and significantly negatively correlated with the NDVI in urban areas. In contrast, human activities in mountainous areas at a certain distance from cities were mostly positively and significantly positively correlated with the NDVI (Figure 6).

Analyzing the correlation between the various land-use types and NDVI from the spatial distribution, woodland had the most significant effect on the spatial distribution of the NDVI in Guizhou Province, followed by changing areas, cropland, and built-up land. In contrast, grassland, water bodies, and unused land had less of an influence on the spatial distribution of the NDVI in the study area (Figure 7D).

The land-use information atlas was examined from 2000 to 2020 (Figure 7A). The transferred land-use types represented 31.05% of the total area of Guizhou. The conversion of farmland to woodland from 2010 to 2020 represented 11.46% of the area of Guizhou. Cropland was converted into grassland from 2010 to 2020, representing 4.47% of the area of Guizhou. These changes were mainly due to human activities, such as converting cropland to woodland and grassland.

Land-use transfer occurred more frequently from 2010 to 2020 than from 2000 to 2010. Conversions between cropland and woodland, and between cropland and grassland, were the most frequent and notable. The areas of woodland and grassland converted to cropland were greater than that of cropland converted to woodland and grassland. Overall, human activities were the main drivers of land-use type changes in Guizhou during the study period (Figures 7B,C; Tables 5, 6).

These findings revealed that forestland, followed by changing regions, had a significant impact on the spatial heterogeneity of the NDVI in the study area. The combined effect of forestland and changing regions significantly influenced the spatial heterogeneity of the NDVI. During the study period, the changing area was primarily the transformation of cultivated land, forest, and grassland, with human activities serving as the primary driving forces. Therefore, by encouraging land-use changes, human activities indirectly affected the spatial heterogeneity of vegetation NDVI in the study area.