Continuous Monitoring of Nighttime Light Changes Based on Daily NASA’s Black Marble Product Suite

Monitoring nighttime light (NTL) change enables us to quantitatively analyze the patterns of human footprint and socioeconomic features. NASA’s Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band (DNB) atmospheric and Lunar-BRDF-corrected Black Marble product provides 15-arc-second daily global nighttime radiances with high temporal consistency. However, timely and continuous monitoring of NTL changes based on the dense DNB time series is still lacking. In this study, we proposed a Viewing Zenith Angle (VZA) stratified COntinuous monitoring of Land Disturbance (COLD) algorithm (VZA-COLD) to detect NTL change at 15-arc-second resolution. Specifically, we divided the clear observations into four VZA intervals (0-20°, 20°-40°, 40°-60°, 0-60°) to mitigate the temporal variation of the NTL data caused by the combined angular effect of viewing geometry and the various kinds of surface conditions. Single term harmonic models were continuously estimated for new observations from each VZA interval, and by comparing the model predictions with the actual DNB observations, a unified set of NTL changes can be captured continuously among the different VZA intervals. The final NTL change maps were generated after excluding the consistent dark pixels. Results show that the algorithms reduced the DNB data temporal variations caused by disparities among different viewing angles and surface conditions, and successfully detected NTL changes for six globally distributed test sites with an overall accuracy of 99.78% and a user’s accuracy of 68.25%, a producer’s accuracy of 66.89% for the NTL change category.


Introduction
Human activities are continuously changing the Earth's natural surfaces and the urban systems, comprising a continuum of socioeconomic and demographic phenomena. Monitoring global human activities and their changes is crucial for understanding global environmental change, sustainability, and socioeconomic status (Leu et al., 2008;Venter et al., 2016;Zhu et al., 2019). Shifts in societies, cultures, economic system structures, policies, technologies, and behaviors are rapidly affecting global ecosystems (Malecki, 1997;Ojima et al., 1994;Steffen et al., 2006). Human footprint expansion and reconstruction actions are mainly driven by socioeconomic changes, such as population growth and the consequential need for natural resources. These drivers have modified the long-term land cover and land use features over large areas (Deng et al., 2009;Turner, 2010). Moreover, disturbance stresses caused by social shocks and behavioral changes (e.g., armed conflicts and gathering events) engendered short-term changes on local scales (Baumann and Kuemmerle, 2016). Data-intensive frameworks for monitoring human-induced land changes at a large scale have become essential to enable a more timely, comprehensive, and deeper understanding of human activity dynamics. This demand calls for the need for reliable, timely, and large-area consistent information. In contrast to the socioeconomic datasets that are at regional or national levels and the field survey measurements that fell short in providing data in the longterm, satellite remote sensing imagery offers a reliable measure of human activity changes at the global scale with much higher temporal and spatial resolutions (Jensen and Cowen, 1999;Xie and Weng, 2016).
To understand human activity trajectories, historically, time series analysis of optical daytime remote sensing data, such as MODIS and Landsat, has been widely applied to capture and categorize humanrelated changes, such as urban conversion and modification , cultivated land changes (Potapov et al., 2022), and forest management (Chen et al., 2019a;Hansen et al., 2013). Advanced time series analysis algorithms, such as Landsat-based detection of Trends in Disturbance and Recovery (LandTrendr) (Kennedy et al., 2010), Breaks For Additive Seasonal and Trend (BFAST) (Verbesselt et al., 2010), Continuous Change Detection and Classification algorithms (CCDC) (Zhu and Woodcock, 2014), and COntinuous monitoring of Land Disturbance (COLD) (Zhu et al., 2020), have been developed to quantitatively retrieve land surface features and changes in urban and non-urban regions based on daytime optical remote sensing datasets. Analysis of the temporal periodicity and variations of the daytime optical observations also improves the extraction of the intra-annual and inter-annual phenology changes (Verbesselt et al., 2010). However, the daytime optical data can only observe the physical land surface changes, and it is difficult to monitor human activity and human behavior changes with limited physical surface changes, such as mobility and power consumption changes. On the other hand, with the direct measurement of both the spatial extent and emission intensity of the artificial light at night, the remote sensing nighttime light (NTL) data provide a unique opportunity to monitor human activity changes (Elvidge et al., 1997;Levin et al., 2020;Zhang and Seto, 2011). Characteristics of artificial nocturnal illumination are often associated with the economic and demographic structures of modern society (Green et al., 2015). The artificial NTL intensity is strongly responsive to the health and development of society (Hölker et al., 2010). Strong correlations have been found between the NTL trends and socioeconomic status , which enable us to estimate the spatial-temporal dynamics of the society based on the NTL changes. Accurate results have been produced by using the NTL data as the major variable for monitoring the urbanization processes (Shi et al., 2014), estimating Gross Domestic Product (GDP) and mapping poverty (Yu et al., 2015), mapping the underserved communities (Machlis et al., 2022;Román et al., 2019), estimating carbon emissions , and analyzing culture behaviors (Román and Stokes, 2015). Human activities are strongly related to NTL changes, however, algorithms that can provide accurate and timely detection of NTL changes based on dense NTL time series (e.g., daily) are very limited. The high temporal variability of the daily NTL satellite observations hampers the accurate detection of NTL changes. The Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band (DNB) data, one of the most valuable remote sensing NTL, provide daily global NTL observations at higher spatial and radiometric resolutions, with significant improvement in the quality, traceability, and consistency over its predecessor: Defense Meteorological Satellite Program's Operational Line Scanner (DMSP/OLS) (Elvidge et al., 2017). However, high uncertainties caused by both the dynamics of the emission sources and the environmental impacts still exist in the VIIRS DNB observations (Coesfeld et al., 2018;Elvidge et al., 2022;Wang et al., 2021) which make time series analysis (e.g., daily) of VIIRS DNB observations extremely challenging. The VIIRS DNB observations are inevitably subject to the extraneous impacts of the angular effect, surface BRDF and albedo, lunar phase, atmospheric effect, cloud and snow contamination, and vegetation canopy . Subsequently, to reduce the data variation, the monthly or the annual average composite images were often recommended in the studies of NTL change (Elvidge et al., 2021;Levin, 2017). For instance, Liu et al. (2019) tracked the seasonal human activity changes during the festival time by using the Seasonal-trend decomposition procedure based on an information ratio index. Zhao et al. (2020) assessed the damages of Hurricane Maria in Puerto Rico by identifying the NTL changes for a particular time based on the seasonal and trend decomposition using Loess (STL). Xie et al. (2019)   Eight VIIRS linear latitude/longitude 10 • ✕10 • tiles were used in this study. The calibration samples were collected from all eight tiles in red squares, the popped-up white squares show the major cities/regions where calibration samples (red dots) are located for each tile. Independent validation samples were collected from the six titles in red stripes (one from each of the six continents) based on the stratified random sampling approach. The row and column numbers represent the "horizontal" and "vertical" number of the tiles, respectively. The background is the Esri ArcMap World image. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the NTL dynamics during the Iraq War with a Sum of City Light (SCL) index that estimated the NTL changes at a regional level. Zheng et al. (2021) mapped urban land changes based on a logistic-harmonic model and provided an in-depth analysis of different urban land change patterns using NTL data. The use of monthly or annual composites of NTL data can substantially reduce the temporal variation (Elvidge et al., 2021;Levin, 2017;Zheng et al., 2022), particularly for developed cities and intra-city areas; however, it will also reduce the temporal frequency of the NTL time series, making it difficult to provide timely information and monitor short-term changes.
The recently released NASA's atmospheric and Lunar-BRDFcorrected Black Marble NTL product (VNP46A2) offers new opportunities for analyzing NTL dynamics by providing daily VIIRS DNB data at 15-arc-second spatial resolution with operational correction for surface reflected lunar radiance (Román et al., 2018). A significant reduction of the temporal variation has been achieved by the daily Lunar-BRDF-corrected Black Marble NTL product . Yet, there are still some remaining factors that could cause large variations in NASA's Black Marble products, such as cloud and snow missed in the Quality Assurance (QA) flag , vegetation phenology, surface albedo (Levin, 2017;Tang et al., 2021), and angular effects from the illuminating artificial lights (Li et al., 2022;Tan et al., 2022), which make their direct usage for time series analysis difficult.
In this study, we proposed a new algorithm for continuous  monitoring of NTL changes based on daily VIIRS DNB observations from NASA's Black Marble standard product suite, which provides accurate and timely detection of human-related changes at night on a daily basis. Specifically, we focused on: (i) filtering out the potential cloud and snow contaminations missed in the Level 3 QA flagging process; (ii) stratifying the data into four viewing zenith angle (VZA) intervals and conducting continuous NTL change detection; (iii) identifying and removing the detected breaks in consistently dark areas.

Study area and data
We selected eight globally distributed regions with different land cover and land use types, and human activity changes, which corresponded to eight NASA Black Marble Level-3 VIIRS tiles (latitude/ longitude extent of 10 • ✕10 • per tile): h10v04, h11v07, h13v11, h17v08, h18v04, h21v05, h29v05, and h32v12 (Fig. 1). The analysis employed all available 15-arc-second spatial resolution daily DNB atmospheric-and Lunar-BRDF-corrected Black Marble NTL radiance (VNP46A2) collected between 2013 and 2021 for NTL change detection. At the same time, the corresponding Top of Atmosphere (TOA) Black  Comparison of the change detection results with different sets of parameters (i.e., number of consecutive anomaly observations to confirm a change and change probability threshold) based on the VZA-COLD and the original COLD algorithm. The red arrow lines show the VZA-COLD results, and the blue arrow lines show the original COLD results. The shapes of the marker represent the selected number of consecutive anomaly observations, the arrows indicate the change direction of the applied change probability from low to high, and the marker face colors show the corresponding F1 scores. A consecutive anomaly observation of 14 and a change probability of 75% were the optimal parameters selected for the VZA-COLD algorithm. (Change Prob.: change probabilities; Conse.: number of consecutive anomaly observations.) (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
Confusion matrix of the validation tile maps in sample counts, area proportions, and accuracies of the change and no-change stratum.

Map Strata
Change Marble product (VNP46A1) was also used to get the VZA and the QA flag for cloud and snow status. To develop and calibrate the algorithm, we manually interpreted 610 calibration samples that have undergone various kinds of human activity changes, such as construction actions, economic growths, gathering events, armed conflicts, power outages, and new streetlights, as well as stable samples without any NTL change, based on the opportunistic strategy around the major cities and metropolitan areas within the selected eight tiles (Fig. 1). Particularly, to explore the complex angular effects, the calibration samples were selected from areas with typical types of local geometry conditions, such as downtown areas with skyscrapers, single-and multi-story residential regions, and areas with dense vegetation canopy. The collected change samples include both abrupt changes with sudden jumps in the NTL radiances that can transfer to the next stable stage in less than one month and transitory changes that have more gradual shifts (e.g., from positive to negative trends) that can last for more than a year.

Methods
We developed an algorithm named "View Zenith Angle (VZA) stratified COntinuous monitoring of Land Disturbance (COLD)" (hereafter called "VZA-COLD"), which was built based on the Landsatbased COLD algorithm (Zhu et al., 2020). The COLD algorithm detects land surface change by continuously updating a time series model to predict future observations. Anomaly observations are identified if the difference between the predicted and observed values are larger than a threshold and breaks of the time series are confirmed if a number of consecutive anomaly observations are identified. The VZA-COLD algorithm uses daily NTL data as inputs and has three major innovations: (i) cloud/snow buffer; (ii) change detection based on observations from four stratifications of VZA; and (iii) consistent dark pixel removal. The workflow of the VZA-COLD algorithm is illustrated in Fig. 2. First, the remaining cloud and snow observations were screened based on the exclusion of cloud and snow edge pixels that are partially influenced by clouds or snow. Second, DNB observations were stratified into four groups based on the VZA intervals to mitigate the variance caused by the compounded impacts from the different viewing geometry and surface conditions to reduce the overall time-series variation. For each VZA interval, NTL changes were monitored based on the change detection framework similar to the COLD algorithm to obtain the individual sets of estimated time series models and detect a unified set of breakpoints among all the VZA intervals continuously. Third, the consistent dark pixels were filtered with a minimum detectable NTL radiance threshold error caused by the remaining consistent dark pixels with low NTL radiance for a pasture field in Wyoming, PA. The red, green, blue, and grey colors indicate the different VZA intervals, the lines represent the estimated models, the small dots are the DNB observations, and the large magenta dots are the detected changes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) of 1.0 nW * cm − 2 * sr − 1 . The VZA-COLD algorithm was calibrated based on the calibration samples, in which the changes detected within the period of ± six months of the calibration sample change intervals were determined as the right call. Accuracy estimators of omission rates, commission rates, and F1 scores that measure the harmonic mean of the precision and recall (Dice, 1945;Sorensen, 1948) of NTL changes were used to evaluate the algorithm performance. A higher F1 score with balanced omission and commission rates indicates a better performance of the algorithm. The final map accuracy was evaluated based on independent validation samples following the "good practice" protocol (Olofsson et al., 2014).

Remove remaining cloud and snow impacted observations
Considering that cloud and snow edge pixels are very likely to be influenced by thin clouds and snow, a spatial buffer was applied to remove these edge pixels. Confident/probable cloud, cirrus cloud, and snow/ice observations were firstly removed according to the QA flags of the standard NASA's Black Marble product. The cloud/snow edge pixels removal was tested by dilating cloud/snow pixels (at 8-connected directions) from 0 to 11 pixels to find the optimal moving window size based on our calibration samples. Fig. S1 in supplementary material shows the accuracies (Fig. S1a) and the percentages of data removed by the cloud/snow buffer for the images collected at the eight tiles ( Fig. S1b) at different moving window sizes. For the moving window size equal to or less than five pixels, both omission and commission errors dropped gradually along with the increase in the window size (Fig. S1a). A decrease in F1 scores and increases in omission/commission errors were observed when the moving window size was larger than five pixels ( Fig. S1a), which is mostly due to the removal of too many clear observations for change detection (Fig. S1b). Thus, the 5 × 5 pixel moving window was selected as the optimal buffer size for masking potential cloud/snow-influenced pixels.

Viewing zenith angle stratification
To reduce the large temporal variances, we stratified the observations by equally dividing the VZA range of 0-60 • into several subsets. We tested the VZA stratification from two to six equally divided groups (based on VZA) and determined the optimal VZA interval number by comparing the cross-group standard deviation and within-group standard deviation of the Black Marble data based on the calibration samples (See Fig. S2 in supplementary material). The three VZA subset strategy (0-20 • , 20 • -40 • , 40 • -60 • ) was selected to maximize the cross-group variance, reduce the within-group variance, and keep the number of VZA subsets small to ensure each VZA interval contains enough observations for time series analysis. This stratification also matches the nearnadir (0-20 • ) and off-nadir (40 • -60 • ) divisions of NASA's Black Marble monthly and annual NTL composite products (VNP46A3/A4) . To maintain the dense temporal resolution of the NTL time series and enable accurate detection of ephemeral changes, a fourth interval that covers the entire VZA range (0-60 • ) was also applied.
Therefore, the DNB observations were grouped into three VZA subset intervals, 0-20 • , 20 • -40 • , and 40 • -60 • , to mitigate the temporal variation for angular affected NTL emission areas, with a fourth group that include DNB observations from all VZA ranges (0 • -60 • ). Fig. 3 shows . This VZA-related DNB radiance disparity is mainly caused by the complex and case-bycase joint angular effects of their viewing angle and local geometry Tan et al., 2022;Wang et al., 2021). For low buildings and open area pixels with streetlights fully shielded, such as rural settlements without adjacent occlusions (Fig. 3a), almost no DNB differences were observed among different VZA intervals. While significant VZArelated disparities were observed in areas with multiple-story buildings and rural settlements surrounded by dense tree canopy ( Fig. 3b-c). In downtown areas with skyscrapers, near-nadir observations have the largest overall radiance and variations (Fig. 3b). In the dense residential areas, however, the off-nadir DNB radiances are dramatically larger than the ones with lower VZA values (Fig. 3c). The inconsistency of the DNB radiances consequently led to inequalities in the sensitivity of monitoring NTL changes, therefore disparities in the occurrences, magnitude, and timing of the changes show up across the VZA. Fig. 4 shows an example with disparate magnitudes and occurrences of DNB radiance changes among the different VZA intervals. The near-nadir time series has a larger magnitude of change in early 2016 than the one with 20 • -40 • VZA interval, while this 2016 change is hardly noticeable from the off-nadir time series with 40 • -60 • VZA interval (Fig. 4).

Continuous monitoring of NTL change
The VZA-COLD algorithm was applied to estimate time series models from DNB observations collected within different VZA intervals while collectively identifying NTL changes across all VZA strata. For DNB observations stratified within each VZA interval, an individual harmonic model (Eq. 1) was estimated to capture the seasonality and trend of the DNB time series, which could greatly reduce the impact of the intraannual (e.g., vegetation phenology and snow) and inter-annual (e.g., gradual economic growth and vegetation long-term growth) dynamics. We tested the models with unimodal, bimodal, and trimodal seasonality (4, 6, and 8 coefficients) based on the Ordinary Least Square (OLS) regression, Least Absolute Shrinkage and Selection Operator (LASSO) regression (Tibshirani, 1996), and robust regression (Hampel et al., 2011;Zhu et al., 2015) to explore the optimal combinations for modeling the daily DNB time series. We observed that models with a single-term harmonics model and based on robust regression had the best results for our calibration samples and were more robust to outliers and less likely to be overfitted. Therefore, the single-term harmonic model (Eq. 1) estimated based on robust regression was selected to predict the overall DNB magnitude, intra-annual seasonality, and interannual trends, which would be used in continuous monitoring of NTL changes.
where, ρ i,x : Predicted DNB value for the ith VZA interval at Julian date x.
x: Julian date. T: Number of days per year (T = 365.25). a 0 : Coefficient for overall value for the DNB. a 1 , b 1 : Coefficients for intra-annual change for the DNB. c 1 : Coefficient for inter-annual change (trend) for the DNB. Continuous change detection was conducted based on the models estimated from each VZA stratification following the COLD algorithm (Zhu et al., 2020), by comparing the actual observations with the model predictions. Breakpoints were identified based on the number of consecutive anomaly observations beyond the applied change probability thresholds. The VZA-COLD algorithm made three major changes compared to the original COLD algorithm. First, due to the high temporal frequency and the large fluctuations observed with the daily DNB time series, VZA-COLD would tolerate one of the observations (except for the first one) not showing up as an anomaly in the consecutive anomaly test. Second, VZA-COLD detected changes based on a set of four time series models estimated from observations of different VZA intervals (instead of observations from different spectral bands), and when a breakpoint was identified by any of the VZA interval models, it would be applied to all the four VZA-stratified models, and thus, dividing their time segments with the same break time. Third, the optimal parameters for the change detection process were different from the default COLD algorithm.
To determine the optimal parameters and evaluate VZA-COLD's performance for the NTL change detection, we compared VZA-COLD with the original COLD algorithm that only has one set of models for all data within the same VZA 0-60 • interval based on different pairs of consecutive anomaly observation numbers and change probabilities. Accuracy estimators of omission rate, commission rate, and F1 score were calculated based on the calibration samples. Fig. 5 illustrates the accuracies of VZA-COLD and the original COLD algorithm when the consecutive anomaly observation number and the change probability vary from 13 to 16 and 60% to 80%, respectively. According to the results, the VZA-COLD algorithm reduced both the omission and the commission rates at a higher change probability, which demonstrates that the VZA-COLD algorithm has mitigated the temporal variability of the NTL time series and detected the NTL changes more accurately than the original COLD. Among all these different pairs of parameters, the consecutive anomaly observation of 14 and the change probability of 75% were selected for the VZA-COLD algorithm to detect NTL changes due to its high F1 score and balanced omission and commission rates. Note that we kept the change detection parameters the same among the different VZA intervals for the VZA-COLD algorithm in this comparison. Considering the different data characteristics between the full VZA interval and the three VZA stratified intervals, we also tested using different parameters for the two different sets of data for the VZA-COLD algorithm. The use of 14 consecutive observations and a change probability of 75% for all VZA intervals still show up as one of the best scenarios (See Table S1 in supplementary material), with top 9 ranked F1 score (< 2% smaller than the best F1 score), smallest omission error, and the smallest number of consecutive observation to confirm a change (faster change detection results). See Section 3 in the Supplementary material for details.

Consistent dark pixel removal
The change probability in VZA-COLD was calculated based on the normalized Root Mean Square Error (RMSE) values from the robust regression, which could be still sensitive to pixels that are consistently dark due to their relatively low temporal variations and result in very small RMSE values in model estimation. Therefore, a slight change in the DNB values caused by outliers could lead to a substantial increase in the final normalized change probability, and results in commission errors. To mitigate this issue, the detected breakpoints with low overall values of DNB and small change magnitudes are considered as low confidence changes that are more likely to be commission errors caused by the outliers. Changes detected over the consistent dark area were identified based on the predicted overall DNB values before and after the change, and the corresponding change magnitudes. A threshold of 1.0 nW * cm − 2 * sr − 1 , which is two times the breakthrough value of the NTL detection limit (L min = 0.5 nW * cm − 2 * sr − 1 ) defined in the daily Black Marble product (Román et al., 2018), was applied to exclude the low confidence changes in consistent dark areas. The detected breakpoints, with before-break overall value, after-break overall value, and change magnitude value all <1.0 nW * cm − 2 * sr − 1 , would be identified as consistent dark pixels and would be removed from the final change detection results. Commissions caused by the scattering light (See Fig. S3 in supplementary material) and salt-and-pepper noise (See Fig. S4 in supplementary material) were substantially removed by this approach.

Accuracy assessment
A total of six VIIRS tiles, including h10v04, h13v11, h17v08, h18v04, h29v05, and h32v12, were selected as the validation tiles to cover all six non-polar continents (Fig. 1). The stratified random sampling strategy (Olofsson et al., 2014;Stehman et al., 2012) was applied for selecting the independent validation samples. Here, the sample population was the total number of NTL pixels, excluding the ocean areas, multiplied by the number of years. Based on the annual change maps from 2013 to 2021, the pixels were assigned to either the change stratum or the no-change stratum for each calendar year. Each sample represents not only a location on the ground but also a place in time. The sample sizes were calculated with the pixel count area estimation of the two strata based on the annual change maps following Olofsson et al. (2014), with target user's accuracies of 70% for the change stratum and 90% for the no-change stratum, and a target standard deviation of 0.01 for the over accuracy. A minimum number of 200 samples per stratum was applied to make sure enough samples can be collected for the minority stratum. Therefore, a total of 901 no-change samples and 200 mapped change samples were selected. Manual interpretation of the validation samples was conducted based on original NTL imagery, high-resolution images from Google Earth and PlanetScope data, International Space Station (ISS) nighttime photos (https://eol.jsc. nasa.gov/), and other available socioeconomic data. A reference sample will be labeled as correct if the mapped category (e.g., change or stable) is the same as it is interpreted in the selected calendar year. Based on the manual interpretation, 8 of the samples were excluded due to the low confidence in interpretation. Thus, a total number of 1,093 samples, of which 900 from the no-change stratum and 193 from the change stratum were selected (Table 1). The accuracy of the change map was then evaluated following the unbiased stratified sampling estimators of overall accuracy, user's accuracy, and producer's accuracy adjusted by the estimated area proportions of the strata (Olofsson et al., 2014;

Results
Quantitative evaluation of the VZA-COLD algorithm performance was estimated based on the 1,093 stratified random samples collected from the validation tiles and following the "good practice" recommendation (Olofsson et al., 2014) and the manual interpretation (See Table S2 in supplementary material). A confidence interval of 95% was applied to present the uncertainties of the unbiased estimators. The overall accuracy of the map is 99.71%±0.16%, with a user's accuracy of 87.18%±2.40% and a producer's accuracy of 68.88%±15.16% for the NTL change category (Table 1). For the 170 samples that were correctly detected by VZA-COLD as NTL changes (Table S2), 35% of them were caused by new constructions, 10% were caused by reconstruction actions, 5% related to mining, 14% related to agriculture processes, 10% related to roads, and 19% are NTL changes without any physical land surface changes (hereafter will be called "other"). The near-nadir interval (0-20 • ) is particularly good at detecting new construction, road, and other change types, while the all data interval (0-60 • ) detected the most changes for reconstruction, mining, and agriculture change types (See Fig. S5a in supplymentary material). The near-nadir interval (0-20 • ) and all data interval (0-60 • ) observations are particularly good at detecting NTL changes that occurred in non-residential areas, and there are no major differences in detecting NTL change in residential areas among the different VZA intervals (See Fig. S5b in supplementary  material). The increasing NTL change samples are more common (approximately three times) than the decreasing ones. Despite correctly identified results, omission errors caused by the undetected reconstruction changes in non-residential areas and commission errors of false detections still existed. For example, Fig. 6a shows an omission error in a factory area with new construction. Due to the high temporal variations of the NTL time series and the subtle increases in the NTL radiances, the algorithm missed the change in 2020. Fig. 6b shows an example of the delayed detection of a transition change of a new factory around 2017 that caused both an omission error and a commission error. Fig. 6c shows an example of a commission error of a residential pixel in New York, the higher NTL radiances affected by snow albedo caused commission errors during wintertime. And Fig. 6d shows an example of another typical type of commission error in a dark pasture area with slightly higher overall NTL radiances than the applied consistent dark pixel removal threshold.
We applied the VZA-COLD algorithm to all the selected tiles ( Fig. 1) to examine its change detection performance for various kinds of human-related NTL changes in different regions. The annual and day-ofyear NTL change maps from 2013 to 2021 were created for the study area. To demonstrate the algorithm's capability in monitoring NTL changes, we investigated a range of urban and peri-urban regions with human activity changes corresponding to different land use, demographic, and socioeconomic typologies . The results, shown in Figs. 7-13, illustrated how the algorithm can accurately capture NTL changes caused by the major types of both short-term and long-term transitions. These factors include but are not limited to, urbanization processes in sub-urban areas, non-residential constructions of public facilities, land cultivation of a new modern agriculture field, redevelopment of a pre-existing urban area triggered by the economic growth, de-electrification derived by the renovation of lighting technologies and environmental policies, power grid loss caused by natural hazards, and armed conflicts. Meanwhile, the identified change areas covered a wide range of human footprints with different land cover and land use types, including the highly populated urban areas, urban green space, suburban and rural areas, agricultural fields, roads, and barren land regions with human activities.
Urbanization with the construction of residential and non-residential developments is one of the most prevalent human-driven land cover and land use changes. Fig. 7 shows the urban expansion process of a new residential community in the suburban area of Melbourne, Australia that converted the agriculture fields into impervious surfaces. A gradual NTL change was captured by the algorithm at this site, which was consistent with the built-up period of this new settlement from 2017 to 2021 (Fig. 7c). Fig. 8 shows the construction actions of a newly built international airport in the suburban areas of Beijing, China. Multiple NTL changes were identified for the new international airport, which aligned well with the different construction stages (Fig. 8c). The timings of these three identified changes agreed with the start of land clearance in February 2016, the finishing time of the major construction in 2018, and the opening of the airport in September 2019. In addition to urban developments, the land cultivation engendered by the food consumption in agricultural land can also be captured by the DNB time series. Fig. 9 shows a new modern organic vegetable greenhouse built with LED plant light at night in the low light area of Canada. According to the time series result of the selected pixel (Fig. 9c), a dramatic increase in the NTL was captured in 2017 after the greenhouse was put into use.
Redevelopment driven by the population and economic growth, (de)-electrification caused by new technologies and environmental policies, natural hazards, and armed conflicts can also temporarily or permanently change the intensities of the artificial NTL over pre-existing developments without land conversions. Fig. 10 shows the redevelopment of urban areas in Abidjan, Ivory Coast. Foreign investments promoted both GDP and population density of the urban environment of Abidjan (Ramiaramanana et al., 2021), which was detected in 2014. The large-scale renovation of LED streetlights in the suburban areas of Milan that were planned by environmental policies encouraged by the International Registered Exhibition in 2015 can also be detected from the annual NTL change map (Fig. 11), in which a significant drop in the NTL radiances was captured in 2014 after the new energy-saving LED streetlights were installed. These newly renovated LED streetlights are more energy efficient with less upward emissions and are typically emitting more blue wavelengths light that fall out of the VIIRS DNB spectral range. In September 2017, Puerto Rico was hit by a powerful hurricane, and the abrupt power outage and gradual restoration were also successfully identified by the VZA-COLD algorithm (Fig. 12). In urban areas with high dynamics of NTL caused by armed conflicts such as the Syrian Civil War in Raqqa, the algorithm identified multiple NTL changes between the stable periods with relatively short durations (Fig. 13).

Discussion and conclusions
The high temporal variation of the VIIRS DNB observations Wang et al., 2021) has made it extremely challenging to monitor and detect NTL changes at the temporal scales necessary to capture the various stages of urban developments (Boucher and Seto, 2009). With the major environmental effects of atmosphere and moonlight removed, the atmospheric-and Lunar-BRDF-corrected NASA's Black Marble product provides new opportunities to use extremely dense (e.g., daily) NTL time series (Román et al., 2018). Using the daily Black Marble data as inputs, we developed new methods to filter, fit, and further reduced the high temporal uncertainties of the daily Black Marble data. As the viewing angle impacts dominate the remaining external uncertainties of the data Wang et al., 2021), we mitigated and simplified the complex effects of the wide ranges of sensor viewing angles and different surface geometry conditions by stratifying the DNB time series based on four sets of VZA intervals. Compared with correcting the viewing angle effect (Tan et al., 2022), the stratification strategy avoided the need for collecting local 3-dimensional structure data and the results would not be influenced by the potential bias introduced in the correction approaches. The combined knowledge derived from VZA interval subsets, and all data models reduced the variance across VZAs without decreasing the temporal frequency of the data.
Among all the four VZA intervals, the all data interval (0-60 • ) identified the most changes (34% of all changes) for all the study tiles, which is not surprising as a large proportion of the NTL changes are either having homogeneous NTL emissions among different viewing angles or showing abrupt and short-term NTL radiance changes, and the higher the data density the more transient changes could be captured by the all data interval models. On the other hand, 66% of the changes were detected by the other three VZA subset intervals, with the near-nadir interval (0-20 • ) capturing the second largest amount of changes (23%), and the other two intervals share similar and slightly less amount of detected changes (22% for the 20 • -40 • interval and 21% for the 40 • -60 • interval). This suggests that even though the all data interval can provide much denser (three-fold dense) NTL observations, the large data variation will still impact their use of change detection. The three stratified intervals showed better performances for two-thirds of the change detected, which further demonstrated the benefits of the stratified models that are useful for pixels with strong angular effects and changes that can only be observed under specific viewing angles.
In addition to the detected NTL change patterns (change time and change location), other change metrics, such as the different VZA intervals that captured NTL changes, change magnitude (NTL differences between model predictions and the actual observations), change direction (increase or decrease of NTL between the model predictions and the actual observations), and coefficients of the time series models can also provide additional information about the change events (Kyba et al., 2021;Wang et al., 2022). For example, Fig. 14 illustrates the maps of the NTL changes identified at different VZA intervals (Fig. 14a), the corresponding change magnitudes and change directions (Fig. 14b), and a high-resolution image (Fig. 14c), for the area around the Gaza Strip that experienced air strikes in 2014. With smaller temporal variations compared to the all data interval, the 0-20 • and 20 • -40 • intervals detected more change in the densely populated areas with multi-story buildings in this region. For the other areas with lower and sparser developments, most of the changes were detected by the all data VZA interval (0-60 • ). For the predicted change magnitudes (Fig. 14b), the whole region of the Gaza Strip experienced substantial decreases in the NTL radiances during this war with clusters of more severe declines around Gaza City, Nuseirat Camp, Khan Yunis, and Rafah. The patterns of the predicted decreasing NTL magnitudes agreed with the key areas that suffered more from the air attacks (Fig. 14c).
The NTL data enables us to better understand human activity patterns and urban infrastructural transitions from a unique perspective, which is important for linking development patterns to ecological, economic, and social health (Elmqvist et al., 2013;Ranis et al., 2000;Stokes and Seto, 2019). While daytime optical remote sensing has tracked urban land cover change for decades (Miller and Small, 2003;Weng, 2012;Zhu et al., 2019), there is a need for new data and analytical tools to capture the land use patterns and tie them to socioeconomic factors. By providing the direct measure of the artificial light at night, the NTL data looks at the human footprint influenced by human activity and behavior aspects rather than the surface physical patterns captured by the daytime optical observations. Due to the different aspects observed by nighttime and daytime data, disparities show up in the location, time, and types of the detected changes. Figs. 15-16 show the comparison of change maps for Pyeongchang, South Korea and Puerto Rico, generated based on the NTL data using the VZA-COLD algorithm (hereafter will be called "nighttime change map") and the Harmonized Landsat Sentinel-2 (HLS) surface reflectance data (hereafter will be called "daytime change map") (Claverie et al., 2018) using the COLD algorithm, respectively. For the Pyeongchang site, although both daytime and nighttime change maps captured changes triggered by the urbanization process for the 2018 Winter Olympics, these changes have different spatial patterns, timings, and meanings (Fig. 15). The daytime change map shows places where the land surface has been physically changed due to the newly built infrastructures, and with the higher spatial resolution data used (30-m), we are able to see the detailed spatial pattern of the change object (e.g., building boundaries). On the other hand, the changes detected in the nighttime change map cover a larger area during the 2018 Winter Olympics, which are corresponding well to this large gathering event. The daytime change map shows the timing of the clear land actions (Fig. 15e) while the NTL change map shows the timing of the installation of electrification infrastructure of the building with a lag of around one year (Fig. 15d). NTL changes can also provide additional information about human activity and infrastructure usage patterns after the new developments were built. For example, after the 2018 Winter Olympics, the built infrastructure was still there but much less used which can be detected from the NTL trajectories by the VZA-COLD algorithm (Fig. 15d). When it comes to the monitoring of natural hazard impacts in Puerto Rico, the daytime and nighttime change maps show even larger differences in their spatial distribution and change types (Fig. 16). For the daytime detection results, changes were detected over the developed (P1 in Fig. 16d) and natural surfaces (P2 in Fig. 16f) that were damaged by the hurricane. The NTL changes showed more on the human development area and captured the power grid variations, power outages, and the corresponding recoveries over the settlements with or without building damages (P1 in Fig. 16c, P3 in Fig. 16g). Generally, the NTL data and daytime data capture the changes from different perspectives and scales, in which the daytime HLS data can provide detailed spatial information of the physical land surface changes while the NTL data can provide human activity and behavioral information at a relatively coarser resolution. By combing the nighttime and daytime changes, we can better differentiate the human activity changes from the physical land surface changes and have a more comprehensive picture of the coupling humanenvironment systems (Chen et al., 2019b;Zhang et al., 2015).
VZA-COLD also has some limitations. For example, high latitude areas with large coverage of winter snow/ice could still lead to commission errors (e.g., Fig. 6c) if the cloud/snow QA and their buffers did not exclude winter snow/ice well (Frey et al., 2020). This error could be reduced in the upcoming Collection 2 Black Marble products after reprocessing all the data with an improved snow flag derived from the VIIRS snow product with a higher spatial resolution. Second, the applied 1.0 nW * cm − 2 * sr − 1 threshold for consistent dark pixel removal can still miss some consistent dark pixels and result in commission errors (e.g., Fig. 6d). Third, VZA-COLD might miss some of the transition changes with subtle shifts in the NTL radiances (e.g., Fig. 6a) or might be able to detect the change but with some delays (e.g., Fig. 6b). Fourth, as VZA-COLD will model intra-annual changes, these annually repeated NTL changes caused by festivals and holidays will be likely modeled as stable cyclical patterns instead of NTL changes. Therefore, more work is needed to detect these transient and repeated intra-annual changes. Finally, VZA-COLD needs 14 consecutive anomaly observations to confirm a change, and even for the all data interval without any cloud/ snow coverage, it will take at least half a month to fully confirm an NTL change, which may limit its capability for near real-time applications.
In conclusion, we developed a VZA stratified COLD algorithm to continuously monitor NTL changes based on the daily atmospheric-and Lunar-BRDF-corrected Black Marble product. This algorithm enabled the time series analysis of the daily DNB data by reducing the viewing angle introduced variations without decreasing the temporal frequency of the data. VZA-COLD can accurately detect different kinds of humanrelated NTL changes, either with or without land cover conversion, and provide concurrent change information regarding urbanization processes, construction actions, and land cultivation. The results indicate that this method could be applied for operational mapping of global NTL changes at a spatial resolution of 15-arc-second with a daily updating frequency.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The NASA Black Marble operational change product is generated based on this algorithm and will be released to the public in 2023.

Acknowledgment
This work was supported by NASA's Terra, Aqua, Suomi-NPP, and NOAA-20 program grant 80NSSC22K0199, NASA's Remote Sensing Theory for Earth Science program grant 80NSSC20K1748, and the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via 2021-2011000005. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied, of ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for governmental purposes notwithstanding any copyright annotation therein. Fig. S1. Analysis of the optimal cloud/snow buffer moving window size. (a) Accuracies for cloud/snow buffer at different moving window sizes. (b) The percent of data masked by the cloud/snow buffer for images collected at the eight titles at different moving window sizes.
Considering the data density of the all data interval (0-60°) is always higher than the other subsets and considering that the VZA stratification can reduce the temporal variation, here, we utilized rules that: (i) the number of consecutive anomaly observation for the all data interval should always be equal or larger than the ones applied to the VZA subset intervals; (ii) the change probability for the all data interval should always be equal or smaller than the ones applied to the VZA subset intervals. Based on the two rules, we run all different possible scenarios (a total of 588 scenarios) for the two parameters and the table below shows the top 9 scenarios (ranked by F1 score) evaluated based on the calibration data. Compared to the VZA-COLD with the same parameters among all VZA intervals (first row), the VZA-COLD with different sets of parameters for the all data interval and VZA subset intervals only improved the optimal F1 scores by less than 2%, but at the sacrifice of a higher omission error (omission error is usually considered more serious than commission error for change detection). Moreover, the smaller the number of consecutive observations to confirm a change, the faster a change can be confirmed. Thus, we decided to keep using the same number of consecutive anomaly observations of 14 and the change probability of 75% among all VZA intervals for the VZA-COLD. Fig. S3. Consistent dark pixel removal for a scattering light pixel. This pixel is located near a desert highway with improvements in electrification around 2017 in northwestern Saudi Arabia. (a) DNB time series at a dark barren pixel with a commission error before the consistent dark pixel removal was applied. The red, green, blue, and grey colors indicate the different VZA intervals, lines represent the estimated models, small dots are the DNB observations, and the large magenta dot is the detected change (red cross in Fig. S3e) that would be excluded by the dark pixel removal process (red cross in Fig. S3f).