Mapping Impacts of Human Activities from Nighttime Light on Vegetation Cover Changes in Southeast Asia

: It is commonly believed that the impacts of human activities have decreased the natural vegetation cover, while some promotion of the vegetation growth has also been found. In this study, negative or positive correlations between human impacts and vegetation cover were tested in the Southeast Asia (SEA) region during 2012–2018. The Visible Infrared Imaging Radiometer Suite— Day/Night Band (VIIRS/DNB) nocturnal data were used as a measure of human activities and the moderate resolution imaging spectroradiometer (MODIS)/normalized difference vegetation index (NDVI) diurnal data were used as a measure of vegetation cover. The temporal segmentation method was introduced to calculate features of two sets of time series with spatial resolution of about 500 m, including the overall trend, maximum trend, start date, and change duration. The regions with large variation in human activities ( V ‐ change region) were first extracted by the Gauss ‐ ian fitted method, and 8.64% of the entire SEA (VIIRS overall trend < ‐ 0.2 or >0.4) was set as the target analysis area. According to statistics, the average overall VIIRS trend for the V ‐ change region in SEA was about 2.12, with a slight NDVI increment. The time lag effect was also found between vegeta ‐ tion cover and human impacts change, with an average of 10.26 months. Our results indicated a slight green overall trend in the SEA region over the most recent 7 years. The spatial pattern of our trend analysis results can be useful for vegetation management and regional planning.


Introduction
During the last decades of modern industrial societies, the world population has grown rapidly by 1.24% per year, and the rate dropped to 1.10% per year nowadays, still yielding an additional 83 million people annually [1,2]. The increasing population will require more urban or built-up areas for human development, and the global land and environment has been suffering measurable human pressures [3][4][5]. Especially for developing countries, the urbanization becomes the irreversible process during the last decades, and it will occupy a certain amount of vegetation cover, such as croplands, forest, and grassland [3,6,7]. This kind of land-use conversion has resulted in serious ecological problems, such as loss of the ecosystem services [8] and biodiversity reduction [9]. However, some researchers have also found a remarkable growing trend of vegetation cover [10], even in some urban areas [5,11]. This is mainly due to more attention to the environmental and ecological protection and scientific protection schemes, resulting in the vegetation cover increment [11][12][13]. Therefore, the positive or negative correlation between vegetation cover and human impacts is still under debate, and the magnitude of these correlations is ambiguous.
Vegetation dynamics and the influence of human impacts are both long-term and dynamic processes, and satellite-based datasets have provided opportunities for holistic monitoring at a larger spatial and temporal range [5,6,14,15]. Simply, the magnitude of human impacts can be approximately quantified as different land cover/use types from single or integrated diurnal remote sensing data [3,16]. Moreover, the satellite-derived nighttime light (NTL) can provide a uniform and convenient indicator to monitor and analyze the magnitude of human activities [17][18][19][20][21]. The intensity of NTL has proven strong correlations with human activities and socioeconomic parameters at multiple scales [22,23], such as urbanization [18,24], energy consumption [25], regional disparity [26], and declines in nocturnal activities by COVID-19 [19]. The Visible Infrared Imaging Radiometer Suite-Day/Night Band (VIIRS/DNB) nighttime lights data, launched in 2012, can detect nighttime lights with higher spatial and radiometric resolutions, with lighter problems of saturation, blooming, and on-board calibration [22,27]. Similarly, the remote sensing techniques can also be applied to map the large-scale vegetation cover dynamics [28,29]. The vegetation index (VI) from the remote sensing data products can measure the vegetation cover directly, such as Normalized Difference Vegetation Index (NDVI). Among all, Moderate Resolution Imaging Spectroradiometer (MODIS) data have the superior advantages in long-term analysis for temporal consistency and the uniform sensors [30].
According to the existing research, vegetation cover can be influenced by human factors and non-human factors, such as climate and the environment [11,31,32]. In this study, we only focus on the human impacts change that could be detected from NTL data, and areas with little human activity changes are excluded to minimize other influences. To evaluate the correlations of vegetation dynamics and human impact, the core problem is to find out the difference from spatial comparison [9]. However, few researchers have focused on the quantitative assessment of this correlation [16], and most of the existing research studies just directly combine NTL and NDVI datasets, based on negative correlations, for human footprint mapping and related information extraction [33][34][35]. Although some researchers have systematically studied the correlation between NTL changes and NDVI changes, these different analyses are mostly based on short time series or static imageries [11,27]. For these kinds of difference analysis, the derived features or differences are sometimes inaccurate mainly due to the natural vegetation phenology of NDVI data [36] and cloud cover of NTL data [37], along with other atmospheric scattering, geometric errors, and stray light problems. These data problems may result in large interannual variability of change detection when choosing different time nodes, and the timeseries analysis can alleviate these problems [38]. Thus, the time series analysis between NTL data and NDVI data is necessary to reveal the vegetation response to human impacts quantitatively.
The time series analysis has proven successful in vegetation change detection, such as agricultural abandonment [39] and forest disturbance [40]. For NTL data, time series analysis is especially valuable for detecting, monitoring, and estimating the socioeconomic dynamics at various temporal and spatial scales [22,23,41]. For the annual seasonality influence of NDVI and NTL data, there are many existing mature methods or models to acquire the real overall trends, such as the Ordinary Least-Squares model [42], STMAP (Sorted Temporal Mixture Analysis with Post-classification) model [16], and the Theil-Sen median slope method [32,38]. However, for both human impacts and vegetation cover changes, the occurrence date of change may vary greatly in a certain studied time range, and there are sometimes multiple obvious changes with different change directions [40]. Meanwhile, changes of time series are not always gradual and constant, and there are some abrupt changes [43]. Thus, the segmentation of the time series is necessary to identify the temporal pattern [44]. Secondly, time lag effects are common for human impacts on different geographical and environmental processes [45,46], especially for vegetation cover from NDVI datasets [47]. Thus, the temporal segmentation of a time series is also the priority for the time lag estimation by finding the core changes. The existing temporal segmentation algorithms are the LandTrendr (Landsat-based detection of Trends in Disturbance and Recovery) [48], MODTrendr (MODIS-based detection of Trends in Disturbance and Recovery) [40], BFAST (Breaks For Additive and Seasonal Trend) [49], and DBEST (Detecting Breakpoints and Estimating Segments in Trend) [50]. By setting several rules and thresholds, the original time series is divided into piecewise segments with several turning points [39]. Then, different segments are generalized under schemes, and sub trends can be fitted and acquired for each part, respectively [40,48,50]. However, these models all need huge amounts of pre-setting parameters and thresholds. For large scales of time series analysis, an accurate temporal segmentation method should be simplified.
In this study, the human impacts are quantified as the VIIRS/DNB NTL intensity, and vegetation cover is represented by MODIS/NDVI values. Thus, the research aim of this study is to (1) portray the human impact trend by time series analysis of VIIRS/DNB data and identify the VIIRS/DNB change areas during 2012-2018, (2) calculate the related vegetation cover trend by time series of MODIS/NDVI data in VIIRS/DNB change areas, (3) quantify correlation between the human impacts and vegetation cover change by spatial comparison, and (4) evaluate the time lag effect by comparing the start date and duration of the VIIRS/DNB and MODIS/NDVI core change.

Study Area
The study area is the Southeast Asia (SEA) region and contains 11 countries, including Brunei, Indonesia, Cambodia, Laos, Myanmar, Malaysia, Philippines, Singapore, Thailand, Timor-Leste, and Vietnam ( Figure 1, Table 1). The SEA region covers an area of 448.25×10 4 km 2 , with an estimated population of 648.69 million in 2017, which is the one of the most densely populated areas in the world. Meanwhile, the SEA region is also the fastest urbanizing and most rapidly developing area in the world, with an annual average urban population growth of 2.52% and an annual average gross domestic product (GDP) growth rate of 4.59% during 2008-2017. The high-intensive economic activities and rapid human growth exert a great pressure on the vegetation cover. Meanwhile, the SEA region has a tropical climate where vegetation resources are abundant and have little seasonal changes, and it is an ideal place to research the correlation between human influence and vegetation cover change.   The VIIRS/DNB is carried on the Suomi National Polar-Orbiting Partnership (Suomi-NPP), and the National Oceanic and Atmospheric Administration (NOAA) provides the monthly composite products. The products are stored with 15 arc-second (about 500 m at the equator) geographic grids, and they span the globe from 75N to 65S latitude with 6 tiles. The unit of average radiance values for monthly products is nanoWatts/cm 2 /sr, and data accuracy reaches 8 decimal places. The averaging method was employed to generate the monthly products, and some noises can be filtered to exclude in advance, including stray light, lightning, and lunar illumination. Meanwhile, the cloud cover is also removed by VIIRS Cloud Mask (VCM) products. Unfortunately, it is still impossible for several areas to get single good-quality coverage during that month due to consistent cloud cover or solar illumination, which indicates that zero or negative values do not always means no lights observed [37]. The relevant files, with extension "cf_cvg", are integer numbers of cloud-free observation days for monthly composite products, and they are used to crosscheck the data reliability [51]. Finally, 73 monthly images (from April 2012 to April 2018) of "vcm" version of VIIRS/DNB monthly composite are selected, including Tile 3 (75N/60E to 00N/120E) and Tile 6 (00N/60E to 65S/120E).

MODIS NDVI Data
The 463 m resolution 16-day NDVI products (MOD13A1, version 6) were obtained from the Level-1 and Atmosphere Archive & Distribution System (LAADS) Distributed Active Archive Center (DAAC), NASA (https://ladsweb.modaps.eosdis.nasa.gov/archive/allData,last accessed on 5 February 2021). This 16-day composite is retrieved from the daily imagery, and it uses the best available value for each pixel, with low cloud cover, low view angle, and the highest value. The NDVI values range from -1 to 1, and data accuracy can reach 4 decimal places. Larger positive means values (>0.2) probably indicate a larger proportion of vegetation cover, where 0.2-0.4 represents sparsely vegetated or non-vegetated areas, and 0.4-0.8 represents green vegetated areas. To keep more data for trend analysis of relative values, NDVI means values that were in the range from 0.2 to 0.8 were all preserved [43]. The MODIS/NDVI data also have the Reliability layers that record the no data, snow/ice, and cloudy grids with values of -1, 2, and 3, respectively, and they can be applied to extract the reliable grids. To cover the entire SEA region, 19 MODIS tiles are required: h26v06-v07, h27v06-v09, h28v06-v09, h29v06-v09, h30v07-v09, h31v09, and h32v09. Then, the 19 tiles are mosaicked into one imagery for a certain date and converted from sinusoidal projection to the WSG-84 geographical coordinate system, with a pixel size of 15-arc seconds by the nearest resampling method. Thus, the spatial resolution of VIIRS/DNB and MODIS/NDVI data is unified under the same coordinate system. Furthermore, little geometric deviations were removed to ensure that the centroid of corresponding grids for two images were matched properly. Meanwhile, to match the time range of the VIIRS/DNB data, dates from the period 7 April 2012 (2012097) to 24 April 2018 (2018113) are selected for MODIS/NDVI data, and data in this period are all available, including 140 scenes.

Preprocessing
For VIIRS/DNB time series (VTS), the "cf_cvg" file was used to identify the number of clear observations. If the clear observation was zero, the corresponding data value was reassigned to the average value with clear observation (> 0) in 6 neighboring dates (for beginning and ending, the size might be less than 6). If there was no clear observation in 6 neighborhoods, the data value was then set to zero. This ensured the same length of each VTS and the smoothness of the VTS. The same process was applied to the MODIS/NDVI time series (MTS) based on the Reliability layer. Meanwhile, if the values for a certain total VTS or MTS were all < 0, that meant there was no light or no vegetation cover all the time and no need to conduct further analysis for these areas. Then, VTS and MTS were further smoothed using a Savitzky-Golay filter to reduce the remaining outliers with a maximum retention of raw data [52]. Although low-value outliers were removed, there were still several high-value outliers (spike), especially for the VTS, which is ephemeral [48]. Thus, the same 6-size neighborhoods were applied to remove these spikes by defining a standard deviation threshold, which was selected by repeated trials based on several randomly selected samples (Equation (1)). After removal of these outliers, the VTS and MTS can be used for further analysis.
where n is the set of the original 6-size neighborhoods in VTS or MTS and xn is the data value, and are the standard deviation and average value, card is the size of set n. m is the set of effective values in 6-size neighborhoods, and xm is the data value for set m. If the absolute difference between a certain value and average was larger than 5 , the value would be redistributed as the average value of the other effective values.

Segmentation of Time Series
The temporal segmentation process in this study was mainly based on the proposed MODTrendr model [40] and DBEST model [50]. To apply the time series analysis at large scales, we integrated and simplified these two models with less parameters. The major process was as follows: (1) A least squares first-order regression was conducted for the whole time series, where the date was digitized as 1-73 for VTS and 1-140 for MTS. Then, the perpendicular distance was calculated between each point and fitted line in Figure 2(B1).
(2) The point with the maximum distance to the fitted line was set as the potential turning points, and the turning points should also ensure that the period of the segmented time series is larger than 3 months. If the potential point did not satisfy this criterion, the point with the second or third largest distance was selected in Figure 2(B2).
(3) The time series could be divided into two parts by effective turning points. Notably, the valid turning points were the mutual nodes for the connected two sub time series.
(4) Then, the previous steps (1-3) were conducted for each part of the sub time series again in Figure 2(B3), and the iteration process was repeated 3 times for VTS and 4 times for MTS. That meant that the common maximum segments of VTS and MTS were 8 and 16, respectively.
(5) Moreover, the statistical significance test was conducted for each segment based on an F-statistic, with a significant level threshold of 0.05 [40,50]. The local segmentation would continue until statistical significance was detected.
The iteration number decided the number of turning points, and more turning points might reduce the residual sum but would result in over-fitting problems and large calculation complexity. When choosing the optimal turning point number, the Bayesian Information Criterion (BIC) was applied to segmentation results under different iteration numbers to minimize the BIC value [50,53]. After calculation for several samples, 3 and 4 were set as the optimal iteration number for VTS and MTS. , iteration number of 2 as an example) and incorporation (B4-B6). Where "distance_" is the perpendicular distance between each data point and fitted line, "κ" is the slope of the fitted line, and "α" is the included angle between two adjoining fitted lines. The outputs of time series analysis contain four main features, including overall trend, maximum trend, start month, and duration.

Incorporation and Feature Extraction
The segmentation would result in excess problems and the merge or incorporation should be applied for all the segments to simply the results. We used the sorted angle (included angles between adjoining fitted lines, [48]) method to composite each sub time series by a forward stepwise order in Figure 2 (B4). Here, for the integration of both longterm sub time series (VTS>24, MTS>48, about 2 years), an incorporated process happened only when neighboring segments had the same change direction (sign of the change value, increment or decrement). For the composited new time series, the simple least square linear regression and significance test (0.05) was conducted again to get the new trend, until there were no more segment pairs that were satisfied with incorporation rules. In addition, once there was a new calculated trend, this iteration process started from the beginning to find the global minimum angle.
Once this merge process was finished, there were several new vertexes for segmented fitted lines, and the "Overall Trend" was fitted based on these vertexes. For middle vertexes, there might be two values for adjoining parts at one date, and the average value was reassigned. Meanwhile, for merge complement status, the segment with maximum absolute change value and same change direction with "Overall Trend" is set as the core change period in Figure 2(B6). To evaluate the time lag effect, 3 indexes of this core change period were calculated: "Maximum trend", defined as the absolute change value of fitted line for the core change period; "Start date", defined as occurrence date for the core change period; and "Duration", defined as the duration time of the core change period in Figure  2(B5). That means, for a certain VTS or MTS, 4 main selected features are extracted to describe the features.

Time Lag Calculation and Validation
Considering the large areas in SEA, areas with large variation in human activities (Vchange region) were first extracted based on the VTS overall trend, which could reduce the calculating amount and identify the target analysis area. A simple Gaussian fitted method was applied to identify the threshold of this kind of large variation of VTS [54]. Then, the vegetation cover trend was synchronously calculated in these V-change areas. Thus, the time lag effect (Lag_effect) and duration difference (Duration_diff) could be simply calculated by unifying the time range of these two time series (Equation (2)). The length was 143 for 16 day MTS, and it was unified into a time range of 70, which was the length of the VTS.
The validation for time series analysis had no mature scheme, and we conducted two parts according to the common validation schemes. For the NTL data, the other-source reference data were not common, and we randomly selected 300 grids with a large VIIRS/DNB trend, and 100 grids with a nearly zero VIIRS/DNB trend. The validation just focused on the segmentation and incorporation results for the time series itself by manual interpretation of the increasing and decreasing trends. Meanwhile, for vegetation cover, there were several other high spatial resolution diurnal data that could help validate the results. The randomly selected 300 grids were converted into polygon shapefiles and imported into Google Earth Engine. The NDVI from Landsat imageries in these polygons were calculated and the trend direction was compared with our MODIS/NDVI results.

Spatial Comparison and Pattern Analysis
Meanwhile, the Global Moran's I index (−1 to 1) was introduced to evaluate the spatial distribution pattern of the above-mentioned time series features, where larger positive Moran's I values indicated a more positive spatial distribution pattern and smaller negative Moran's I values indicated a more negative spatial distribution pattern [55]. The Moran's I index was first applied to the univariate overall trend, maximum trend, time lag effect, and duration difference. When evaluating the pattern of the relationship between the VIIRS and MODIS trend, a simple combination was conducted based on the quantization of the VTS and MTS trend level (Equation (3)). This grading could also help display and map the VIIRS and MODIS trend relation directly.
where a larger VII_MOD_level value indicates a greener trend (large MOD_index value) with less human impact increase (small VIIRS_index value).
All the above-mentioned processes were conducted by Python 3.6 and related libraries. Considering the similarity between the spatial distribution of the VIIRS/DNB and road network (Figure 1), we also conducted the pattern analysis with the distance to the present road networks by ArcMap 10.2. The road network data were acquired from the Open-StreetMap (https://download.geofabrik.de/asia.html ,last accessed on 5 February 2021) by 2018.04. The major roads, including the motorway, trunk road, and primary road, were selected to evaluate the road influence on the VIIRS/DNB and MODIS/NDVI trends.

Results
We conducted the statistics for the percent of grid counts under different overall VIIRS/DNB trend ranges (Figure 3). Ten of 11 SEA countries exhibited a typical Gaussianshaped distribution, except for Singapore (SGP) with a small land area, which showed consistence with the existing research results [3,56]. Thus, we set the range (-0.2 to 0.4, peak areas of Gaussian curves) as the no VIIRS/DNB change area, and we set other regions (>0.4 or <-0.2) with large variation in human activities (V-change) as the target area for further analysis (Figure 3). About 1.85×10 6 (8.64% of entire region) grids were extracted as V-change regions, with a maximum trend of statistical significance (0.05).

Human Activities Trend from VIIRS/DNB Time Series
The spatial distribution of the V-change regions mainly concentrated inside or around the big cities ( Figure 4A), with a Moran's I value of 0.535 (p<0.01), which meant a strong positive spatial autocorrelation ( Table 2). The magnitude of the overall VIIRS/DNB trend increment presented a radial shape descending from the city center. Meanwhile, there were still some VIIRS trend decreases around the big cities, such as the capital of Indonesia (IDN) and Brunei (BRN), mainly because of the urban renewal or construction of the ecological land according to the Google Earth historic imageries. Interestingly, a VIIRS/DNB increment could also be clearly observed along road networks ( Figure 4A), which meant the construction of road or improvements of road lighting. According to the national statistics in Table 3, we could see that Singapore (SGP), BRN, Vietnam (VNM) and Malaysia (MYS) had the larger average overall VIIRS/DNB trend. These countries also had large GDP growth rate or urban population growth rate (Table 1). Apart from SGP, the large negative VIIRS/DNB trend also occurred in IDN, BRN, MYS, and VNM, which was mainly due to the balance between rural and urban population mobility [57].    The maximum VIIRS/DNB trend showed a similar spatial pattern with the overall VIIRS/DNB trend, with a Moran's I of 0.420, and the scatter plot showed a strong linear relation between them ( Figure 4B, Table 2). The average duration for the VIIRS/DNB core change period was about 48.04 months, where 39.95% had a duration over 60 months, indicating a continuous process of increased human activities ( Figure 4C). Meanwhile, the occurrence of VIIRS/DNB change for about 74.73% of grids started at the benchmark month (2012.04), and the rest showed a relatively average distribution of the start month ranging from 1 to 65 ( Figure 4C). Although the spatial distribution of this duration did not show an obvious pattern, but areas with longer duration tended to be located far away from the city center ( Figure 4D).

Vegetation Cover Trend in V-Change Regions
The same features of MODIS/NDVI time series were extracted based on V-change regions, where the maximum trend for all target grids shows statistical significance (0.05). According to statistics, the average value of the overall VIIRS/DNB trend for the SEA region was about 2.12, with a related MODIS/NDVI trend of 2.27×10 −2 (Table 3). That meant that although the human activities intensified in the SEA region, the vegetation cover still increased. Among all countries, Myanmar (MMR), VNM, and Thailand (THA) contributed the largest rate of vegetation cover increment. Only Cambodia (KHM) and SGP showed a decrease of vegetation cover both for total and positive VIIRS/DNB trend areas (Table 3). For negative VIIRS/DNB trend areas, the corresponding MODIS trend for all countries showed an increment, indicating the promotion of vegetation cover in areas with declining human activity (Table 3). For positive VIIRS/DNB trend areas, the average VIIRS/DNB trend for the SEA region rises to 2.38, while the related MODIS trend dropped to 2.15×10 −2 . Thus, we set these two values as the threshold, combining with the abovementioned V-change threshold, and each trend could be divided into three parts: VIIRS/DNB-high increase (>2.38), slight increase (0.40-2.38), decrease (<-0.20); MODIS/NDVI-high increase (>2.15×10 −2 ), slight increase (0-2.15×10 −2 ), decrease (<0). Thus, the composite of VIIRS/DNB trend and MODIS/NDVI trend could be divided into nine categories ( Figure 5B). For the spatial distribution of these nine categories, the "green" trends were observed to concentrate in THA, Laos (LAO), KHM, and IDN, while the "brown" trends were mainly in THA, MYS, and IDN ( Figure 5A). When away from the center of the residential areas, the human impacts showed a decreasing pattern with an increase of the vegetation cover change, both for positive and negative VIIRS trends ( Figure 5A). Meanwhile, the "green" trend and human activities weakening also showed a certain spatial clustering, with a Moran's I of 0.289 (Table 2). We could also observe that for each country, at least 50% of grids underwent the process of vegetation cover increase or stabilize (green and purple part in Figure 5B), where more than 30% underwent a high increased (>2.15×10 −2 ). Amongst them, MMR and Timor-Leste (TLS) had the highest proportion of grids to undergo the green process. For the whole region, most areas (≈41.42%) were under the slight VIIRS/DNB increase and high MODIS/NDVI increase, and even 6.32% of the area went "green" with heavy human activities increase ( Figure 5B). Only about 1.46% went "brown" due to heavy human activities increase, and the rest (30.76%) of the "brown" areas had little human influence, where similar statistical conditions were also found for each country except SGP ( Figure 5B).

Time Lag Effect of Vegetation Cover Response to Human Impacts
The time lag effect was obvious between the VIIRS/DNB and MODIS/NDVI trend in the SEA region, with an average of about 10.26 months (Table 4), comparing with 2-year time lag effect of non-human impacts [45]. This time lag effect could measure the feedback time of the vegetation cover change to the human activities. Meanwhile, the average duration of the MODIS trend was about 9 months shorter than that of the VIIRS/DNB trend for SEA region ( Table 4). The spatial distribution of the time lag effect also showed a less regular pattern, with a Moran's I value close to 0, but areas closer to cities were more likely to have a smaller time lag effect ( Figure 6A, Table 2). That meant more populated area tended to have faster responses to the human activities and ecological degradation. Amongst them, 37.10% of the grids had no time lag effect, and the start month of the MODIS trend for 44.79% of the grids were posterior to that of the VIIRS/DNB trend (Figure 6C). The spatial distribution of duration difference between the start month of VIIRS/DNB and MODIS/NDVI maximum trend showed a uniform pattern due to the limit of the study period ( Figure 6B). We also conducted the bivariate kernel density analysis between the time lag effect and duration difference. The plot showed two straight lines with high kernel density: a line with zero lag effect and a line with a slope of negative one ( Figure 6D). The sum value of the diagonal line was about zero, which represented the same end time (end of research time, 2018.04) and a large probability of these two trends ongoing. While the vertical line represented that although the start time of the VIIRS/DNB and MODIS/NDVI change was likely to be the same, the duration difference was various, and the MODIS/NDVI change was likely to have longer durations ( Figure 6D).

Validation and Pattern Analysis
For NTL data, all 400 sites showed satisfied and reasonable results between manual interpretation and time series analysis. When comparing the MODIS NDVI trend of Landsat NDVI change, 300 samples also showed a full accordance in change directions and similar change magnitude, indicating the reliability of our proposed segmentation and incorporation method. Apart from samples, we also compared our results with other annual land cover products [58], and the difference by annual products showed a good accordance with our trend analysis. For pattern analysis, the V-change grids indeed showed a close relation with major road networks in SEA regions ( Figure 7A,B). The average value of the VIIRS/DNB overall trend to the major roads exhibited a decrease when the distance to the major roads increased ( Figure 7C). That meant areas closer to the major roads tended to have larger changes of human activities. For NDVI, the change rule was just the opposite, which meant that most of the green trend areas were far from the major roads ( Figure 7C). This simple analysis also proved the correlation between the human impacts and vegetation cover change, where human activities along major roads were intensive and major roads could stimulate human activities change.

Mapping Human Activities from Nighttime Lights
Mapping human activities is always a hot study issue, and there are several perspectives to map the human activities. Among all, the human mobility dataset, such as call detail records and social media data, can be the most accurate data source for the intuitive utility of geo-locations [59,60]. However, these kinds of data are hard to acquire for large scales due to some privacy policy, and huge amounts of data also make the processing extremely difficult. Alternatively, the open-access remote sensing dataset has naturally become the optimal choice to map long-term and large-scale human activities. Meanwhile, as a comprehensive problem, the human impact mapping should contain multiple factors, including population, infrastructure, and land cover. Thus, many human activities products are mostly based on multi-source diurnal data, such as terrestrial human footprint [3], Global Human Footprint [61], and Global Urban Footprint [62]. However, this integration process will bring some system errors, and unification of the temporal range and spatial resolution for all input data is sometimes tough. The NTL data can provide a reliable single indicator to represent the intensity of the human activities from the perspective of nocturnal activities [18,22], although the NTL intensity is not totally consistent with real human activities in some rural or urban-fringe areas [35,63]. Meanwhile, the trend analysis of NTL data can also provide some potential for further analysis. For example, the positive and negative VIIRS/DNB trend can be seen as an indicator to evaluate the magnitude of urban-rural population mobility or land-use conversion [64].

Human Impacts on Vegetation Cover Changes
Understanding vegetation cover dynamics is essential for land surface systems, environment protection, and human daily life, and it is influenced by multiple human and non-human factors [11,31]. Human impacts, mainly the urban process, have directly converted vegetation land into urban land and inhibited vegetation growth [6]. Meanwhile, the "green" trends were also found due to promoting growth and policy changes [5,10,11]. NTL data can only detect partial human activities change, while other human activity forms cannot be detected, such as cut down and plantations [65,66]. Thus, the human impacts from NTL data in this study could be concluded as one kind of representation index for the urbanization process, not for all the human activities. According to our study results, the green trends were found in the areas of the SEA region with increases in human activity (Table 3), indicating more attention to the human living environment in urban areas. This result was consistent with several recent mainstream findings, although the climatic factors are not considered in this study [10][11][12][13]. Thus, the simple assumption of the negative correlations between VIIRS/DNB nocturnal and MODIS/NDVI diurnal data seemed inappropriate. This result also reminded us of careful correlation analysis before data integration.
When comparing different combination of VIIRS trend and MODIS trend, the grid density heat map revealed a rhombus pattern for most countries and the whole SEA region, where the center has the highest grid density (Figure 8). For different countries, IDN, MYS, Philippines (PHL), THA, and VNM showed a more similar pattern with the whole region, while other countries showed a constrictive pattern due to little human activities change or small land area. For the whole region, the density distribution showed a symmetry both on the x-axis (VIIRS/DNB trend) and y-axis (MODIS trend), where negative VIIRS/DNB trend areas represented a tighter distribution ( Figure 8). Meanwhile, this rhombus pattern had four obvious boundaries, and the trend value of both VIIRS/DNB and MODIS/NDVI were basically inside these boundaries ( Figure 8). These boundaries could be seen as a reasonable and common balance between human activities and related vegetation cover change.

Influential Factors of Spatial Difference
Human impacts and vegetation cover change did not show an obvious spatial pattern ( Figure 5A), but generally, "brown" trend areas were closer to the city center, where the city functions were intensive and there was little land for vegetation. For areas far from the center, there was enough land for ecological use even with heavy human impacts. For different countries, the "green" and "brown" trend difference can be attributed to multiple factors, such as economic structure, population changes, policy making, and climatic changes [10,11,[67][68][69]. For example, THA had the highest green trend due to dependence on the agricultural economy, especially in urban and peri-urban areas [70]. The largest "brown areas" in THA and MYS were mainly due to a high urbanization rate, population growth rate, and economic development (Table 1). For the time lag effect, we could see that SGP, MMR, and BRN had the least time lag effect (Table 3), while for KHM, LAO, and THA, the vegetation response time was larger than 12 months. Meanwhile, the duration difference was obvious in BRN and SGP, with a larger or nearly the same duration for the MODIS/NDVI trend than the VIIRS/DNB trend. A smaller time lag effect and larger duration for a MODIS/NDVI increase might indicate a better ecological and environmental adjusting ability to the human activities. Apart from the above-mentioned factors, different protection schemes or policies of vegetation and concepts of urban planning can also result in this difference [68,69]. For example, in SGP, the rooftop garden is popular and can make urban areas greener [71], and BRN and SGP always have the high-level rate of urban forest cover with little construction in recent years.

Evaluation and Validation of Time Series Analysis
For a certain time series, the trend, occurrence date, and duration are the most important and core features. The temporal segmentation is the efficient feature extraction method for time series analysis, but it also brings some user-defined parameters or thresholds. Different artificial settings and choices will greatly influence the results, and lacking efficient reference data will magnify this problem. The proposed temporal segmentation method in this study needed less input parameters by omitting some steps of the DBEST and MODTrendr algorithms. By comparison with the same samples, it improved the calculation efficiency by reducing about 2/3 of the calculation time, also with good performance. There were also several other means to enhance the accuracy of trend analysis; for example, the object-based time series analysis was recently developed to take full advantage of imagery spatial and texture information [36,72]. However, for medium and low spatial resolution data, texture information and geographical object size are not enough for image segmentation. Thus, the pixel-based time series analysis was still the optimal scheme for VIIRS/DNB and MODIS/NDVI data. The validation of the time series analysis results was also a hard problem for time series analysis, and few studies gave the efficient validation schemes. The simple subtraction validation was sometimes contrary to the original idea of time series analysis, and the efficient validation could only be based on the qualitative assessment of the preselected reference data [73]. Thus, although our results showed a good performance with manual interpretation, our results were still the experimental consequence without rigorous validation, and they should be checked by more field data.

Conclusions
The negative and positive correlations between human impacts and vegetation cover are being widely discussed: the occupation of some ecological green land versus promotion of the vegetation cover by the ecological awareness and technology. Thus, Southeast Asia (SEA), with a rapid population and GDP increment, was selected to test this assumption over the last 7 years (2012-2018). The VIIRS/DNB NTL time series represented the change of human activities, and the MODIS/NDVI time series represented the vegetation cover change. The modified temporal segmentation was introduced to analyze and extract the features of these two sets of time series, including the overall trend and maximum trend, start date, and duration for the core change period. The Gaussian fitted model was first used to extract regions with large variation (VIIRS overall trend <-0.2 or >0.4) in human activities (V-change region) as the target analysis areas, about 1.85×10 6 grids and 8.64% of the entire SEA. The spatial comparisons revealed that the average VIIRS/DNB overall trend was about 2.12 in the SEA V-change region, with an average MODIS/NDVI trend of 2.27×10 −2 , indicating a slight green trend under the human impacts increment. When away from the center of residential areas, the human impacts showed a decreasing pattern with an increase of vegetation cover change. Meanwhile, the time lag effect was also discovered between the VIIRS/DNB and MODIS/NDVI trend, with an average of about 10.26 months, although the average duration of the MODIS trend was about 9 months shorter.
The results reveal the correlation between the human impacts and vegetation cover, and they can provide guidelines for the integration of similar datasets. However, there are still some points that deserve improving for this study. The results in this study only focused on the absolute trend value in a relatively short period. The DMSP (1992-2013) NTL data can be integrated to conduct a longer trend analysis, and other kinds of data can also be combined to improve the accuracy of the human impacts estimation, such as social media datasets. Meanwhile, the driving forces behind these changes should be evaluated, such as climatic changes, socioeconomic influences, and the vegetation management policy. Furthermore, although the vegetation cover changes were limited in areas with great change in human activities, the influence of non-human factors should also be removed to reduce the error. Moreover, the temporal segmentation, incorporation, and validation method at finer scales could also be promoted.