The Visible Infrared Imaging Radiometer Suite (VIIRS) Day Night Band (DNB) on board the Suomi NPP satellite now provides almost a decade of daily observations of night light. In comparison to the Defense Meteorological Satellite Program (DMSP) Operational Line Scanner (OLS) night light imagery, VIIRS DNB provides greater dynamic range (14 vs 6 bit), higher spatial resolution (∼0.7 vs 5 km), on board calibration and greater sensitivity for low light imaging (Elvidge et al., 2013). All of these features allow VIIRS to detect a much greater diversity of night lights than DMSP-OLS was able to resolve. In addition, VIIRS imagery is available as individual swaths and daily, monthly and annual composites. The greater temporal frequency of sampling, without the degree of temporal averaging of annual composites, makes it necessary to consider the distinction between apparent temporal changes of night light related to the imaging process and actual changes in the underlying sources of the night light being imaged.

Temporal changes in imaged night light arise from a variety of factors related to ambient phenomena (e.g. stray light, lunar cycle, aurora and lightning) (Elvidge et al., 2017; Ji et al., 2018), atmospheric effects (Fu et al., 2018; Román et al., 2018), view geometry (Li et al., 2019), overpass time (Li et al., 2020), background reflectance (Levin 2017; Levin and Zhang 2017), (Chen et al., 2019) and instrument calibration/drift (Zeng et al., 2018), as well as actual changes in terrestrial light sources. Actual changes in light sources also arise from a variety of factors such as disruptions of electricity supply (Kohiyama et al., 2004; Cao et al., 2013; Mann et al., 2016), conflict (Li et al., 2013; Li and Li 2014; Levin et al., 2017; Li et al., 2018), cultural and religious activity (Roman and Stokes 2015), lighted infrastructure development (Kuechly et al., 2012; Hale et al., 2013; Small and Elvidge 2013; Levin et al., 2014), and gas flaring (Elvidge et al., 2009; Zhang et al., 2015; Elvidge et al., 2016; Coesfeld et al., 2018). While the vast majority of research applications focus on changes in actual light sources, the potential contribution of non-source phenomena are often not considered in analyses of night light change. The most common approach to night light change detection involves direct attribution of observed changes to the phenomenon of interest. Implicit in this approach is the assumption that other forms of actual and apparent change in the light source are negligible or non-existent. An alternative approach is to characterize the spatiotemporal variability prior to deductive attribution of causation so that the attribution can be made in the context of the full range of spatial and temporal variation.

The primary objective of this study is to characterize VIIRS night light variability over a range of spatial and temporal scales to provide a context for interpretation of night light changes observed on both subannual and interannual time scales. Specifically, to introduce a robust methodology for characterization of VIIRS’ spatiotemporal variability that can be used to distinguish multiple sources of apparent and actual change in night light. The strategy implemented in this study is based on the combined use of low luminance thresholds and Empirical Orthogonal Function (EOF) analysis of night light time series on both subannual and interannual time scales.

The Visible Infrared Imaging Radiometer Suite (VIIRS) sensor was launched on board the NASA-NOAA Suomi satellite in 2011. The day/night band (DNB) of the sensor collects low light imagery in a 3,000 km swath at a fixed resolution of 742 m with an equator overpass time of ∼1 AM local time. Individual VIIRS acquisitions are often composited to exclude clouds and intermittent sources like fires. In comparison to DMSP, VIIRS provides higher dynamic range, on-board calibration, and multiple optical bands that can be used to distinguish different light sources. More detailed descriptions of the data, products, and applications of VIIRS imagery are given by (Elvidge et al., 2013) and (Miller et al., 2013). The VIIRS monthly mean night light composites and cloud free coverages used in this study were produced by the Earth Observation Group at the Colorado School of Mines (https://payneinstitute.mines.edu/eog/). All analyses in this study use the stray-light-corrected monthly mean radiance product. Because VIIRS radiances typically span four orders of magnitude, all analyses are performed using Log₁₀ (radiance).

This analysis is based on a combination of temporal moments, spatial correlation and Empirical Orthogonal Function (EOF) analysis. Computation of the temporal mean (μ) and temporal standard deviation (σ) of monthly mean radiance and number of cloud free acquisitions per month allows for the combined use of moment composite maps and moment spaces to illustrate the relationships between monthly means and standard deviations. Spatial correlation matrices computed for all pairs of monthly mean radiance images provide a statistical measure of similarity of spatial distributions of night light brightness. Temporal moments and spatial correlations provide complementary aggregate metrics of variability. While these metrics offer the benefit of intuitive interpretability, a more comprehensive depiction of this variability is given by a combined spatiotemporal analysis. Empirical Orthogonal Function analysis, originally developed for statistical weather prediction (Lorenz 1956), is now a standard tool for analysis of spatiotemporal patterns and processes. Overviews of the use of EOF analysis in oceanography and meteorology are given by (Preisendorfer 1988; Bretherton et al., 1992; von_Storch and Zwiers 1999). The combined use of EOF analysis with temporal feature spaces and temporal mixture models, with application to DMSP night light time series is described in detail by (Small 2012) and (Small and Elvidge 2013).

EOF analysis uses the principal component transform to represent spatiotemporal patterns as orthogonal modes of variance. Rotating the spatiotemporal coordinate system to align with orthogonal dimensions of uncorrelated variance allows any location-specific pixel time series P _xt in an N image time series to be represented as a linear combination of temporal patterns, F, and their location-specific components, C, as: P xt = ∑ i = 1 N C i x F i t (1) where C _ix is the spatial Principal Component (PC) and F _it is the corresponding temporal Empirical Orthogonal Function (EOF) and i is the dimension. EOFs are the eigenvectors of the spatiotemporal covariability (either covariance or correlation) matrix that represent uncorrelated temporal patterns of variability within the data. The PCs are the corresponding spatial weights that represent the relative contribution of each temporal EOF to the pixel time series P _xt at each location x. The relative contribution of each EOF to the total spatiotemporal variance of the observations is given by the eigenvalues of the covariance matrix. N is the number of discrete dimensions represented by the time series of observations. Principal Components are uncorrelated but not necessarily independent–unless the data are jointly normally distributed. In systems in which the same deterministic processes are manifest at many locations, but stochastic processes are uncorrelated, the variance of the deterministic processes may be represented in the low order PC/EOF dimensions while the stochastic variance may be relegated to the higher order dimensions (Preisendorfer 1988). If a clear distinction can be made between a small number of physically meaningful EOFs (or PCs) distinct from a continuum of uninterpretable EOFs (or PCs), this can provide a statistical basis for attribution of deterministic and stochastic components of an image time series. However, the transformation is purely statistical so there is no guarantee that the attribution will be physically meaningful, or even exist at all.

Computation of the temporal means and standard deviations for each pixel time series provides the bivariate distribution of temporal moments necessary to assess the relationship between annual mean brightness and monthly variability of brightness. Computation of pairwise spatial correlations for each same month pair of images using Pearson’s product-moment correlation coefficient provides the 2019–2020 correlations with which to identify the cross year and cross season lowest correlation pairs. Temporal EOFs and corresponding spatial PCs are computed using the eigenvectors of the covariance matrix of the time series of monthly mean radiance. The mean is removed before the data are projected onto the transpose of the eigenvectors so the first PC generally represents the spatial variation in overall brightness.

Temporal moment maps illustrate two of the most pervasive sources of subannual variability in VIIRS night light. Figure 1A shows the latitude-dependent variability in solar illumination (and lack thereof) for Europe, southwest Asia and north Africa (henceforth Eurasia) and Figure 1B the effect of summer monsoon cloud cover for southeastern Asia. In both images, the red channel represents the temporal mean of 24 monthly radiance composites for 2019 + 2020 while the cyan channel represents the temporal standard deviation. The inset maps show the corresponding quantities for the monthly number of cloud free observations contributing to each monthly composite. Both maps for both regions are strongly bimodal, indicating strong regional distinctions between areas with low mean brightness and high temporal variability and areas with high mean brightness and low temporal variability. Areas with high mean brightness and high variability are limited to anthropogenic night lights associated with human settlements and gas flaring. Both high latitude and monsoon summer gaps in coverage result in at least 1 month without sufficient data to detect even bright lights. Examples are shown in the inset in Figure 1A. Offshore, lights from fishing boats also results in high variability and low mean brightness as a result of fleets’ seasonal mobility. In both Eurasia and southeastern Asia the temporal moment spaces (inset) show the bimodal distributions explicitly as two nearly orthogonal limbs corresponding to the aforementioned geographic partitions. The bright + variable night lights appear as diagonal spurs extending away from the high variability limb of each distribution. The high variability, low brightness limb corresponds to large areas of highly variable background luminance. The lower, horizontal limb of each distribution represents the more stable anthropogenic night lights corresponding to human settlements and other lighted development. Note that in both distributions, the range of standard deviation diminishes monotonically with increasing mean brightness on the lower limb. Similar distributions are observed elsewhere on Earth (Small 2021).

Spatial correlation matrices of VIIRS monthly mean radiance composites for 2019 + 2020 show correlations between 0.7 and 0.9 for all pairs of months. Heteroskedasticity of the mode of the distributions results from month to month variability in background luminance which reduces the spatial correlation for large numbers of low luminance pixels. Tri-temporal color composites of mean radiance for the months with the three lowest spatial correlations are shown in Figure 2. Areas with similar brightness in all 3 months appear shades of gray. Color implies change. The upper composite shows the anomalously low brightness for large areas of Russia, the Alps and Anatolian plateau in December 2020 as shades of red/orange and a narrow high latitude zone of anomalously high brightness in May 2020 as green along the top of the panel. The lower panel illustrates the effect of seasonal differences. This composite shows the periurban peripheries of lighted developments on the Russian steppe, Alps and Anatolian plateau with higher brightness in January 2020, consistent with the presence of high albedo snow in winter and low albedo vegetation in summer. This is consistent with the seasonality observed by (Levin 2017) and (Levin and Zhang 2017), although it should be noted that this variability is limited to periurban areas with lower luminance–and presumably much larger areas of illuminated snow and vegetation than the brighter urban cores. As with the moment spaces in Figure 1, the inset scatterplots in both panels show strong heteroskedasticity with much greater dispersion at lower brightness levels, diminishing monotonically as brightness increases. It is also noteworthy that the considerable skewness magnitude and polarity of the lower tails of all these bivariate distributions varies among all pairs while the upper tail of each distribution remains symmetric about the 1:1 line.

The seasonal periodicity indicated by the monthly spatial correlation matrix suggests that interannual changes should be most apparent when comparing the same months of different years. Figure 3 shows scatterplots comparing monthly mean radiance for 2019 and 2020 for the Eurasian region shown in Figures 1, 2. As with the subannual comparisons in Figure 2, the interannual bivariate distributions are all strongly heteroskedastic with skewed lower tails and symmetric upper tails. There is no obvious consistency in the month-to-month skewness variations, while the upper tails are consistently symmetric for radiances > 10^0.5 nW/cm²/sr. Similar distributions are observed elsewhere on Earth (Small 2021).

The spatial principal components (PCs) of the 24 monthly mean radiance composites from 2019 + 2020 show the spatial distribution of the dominant modes of spatiotemporal variability for Eurasia driven by the contrast between large numbers of low luminance pixels and much smaller numbers of much brighter pixels (Figure 4). Specifically, a strong contrast between the bright stable anthropogenic lights associated with settlements and the high variability of background luminance in non-lighted areas. Again, the strong contrast between higher latitude and mountainous areas where background reflectance changes seasonally, and lower latitude deserts with negligible seasonality of background reflectance is apparent. This contrast is reflected in the bimodal structure of the temporal feature space of PCs one and 2. The structure of the feature space is similar to that of the moment space as the PC1 represents overall brightness varying in space while PC2 represents the dominant mode of temporal variability.

In order to suppress the effects of variance heteroskedasticity on the PC rotation, the 24 months image time series was rotated with all pixel time series with radiances less than 10^0.5 nW/cm²/sr masked. Limiting the rotation parameters to those time series with larger radiances effectively focuses the EOF analysis on the variance of anthropogenic light sources. Figure 5 shows a tri-temporal composite of the three low order PCs from the masked time series. As with the other composites, the effect of background reflectance of land cover is clearly apparent–even though it is limited to the larger periurban areas at the scale of the figure. However, it is also clear that there is much greater variability in most of the urban areas, indicated by a greater variety of colors resulting from varying combinations of the low order EOFs. The inset EOFs show a clear distinction among spatial variations in overall brightness (PC1 – red), seasonal variations (PC2 – green) and interannual variations (PC3 – blue). Note that all 3 PCs have both positive and negative values, implying that each corresponding EOF occurs in both polarities. Hence, the upward trend of EOF3 would correspond to decreases in brightness in pixels with negative PC3 values. The overall variance partition derived from the eigenvalues of the covariance matrix are 88, 2 and 2% for dimensions 1, two and three respectively. This suggests that almost 90% of the spatiotemporal variance is associated with spatial variations in average brightness while 8% of variance is associated with the remaining higher (>3) order modes representing stochastic variability and only 2% each for seasonal and interannual temporal variability. With background luminance suppressed, most anthropogenic night light detected by VIIRS appears to be remarkably stable in at monthly time scales.

The spectral feature space of the three low order dimensions of the low-luminance masked image time series shown in Figure 6 is consistent with both the temporal moment spaces (Figure 1) and the bivariate brightness distributions (Figure 3). As overall brightness (PC1) increases, both seasonal (PC2) and interannual (PC3) variance decrease monotonically. The inset time series taken from the upper edge of the envelope of PC1/PC2 distribution show considerable month-to-month variability superimposed on the seasonal cycles, consistent with a variance partition between deterministic (PC ≤ 3) and stochastic (PC > 3) variance. The PC3/PC2 space is effectively a phase plane in which different combinations of positive and negative PC weights for EOFs two and three represent different combinations of seasonality and interannual change between 2019 and 2020. Note that the example time series 2 and 4 are not actually monotonic trends as EOF three depicts, but abrupt increases and decreases. This is consistent with the way that abrupt changes are often represented in low order EOFs (Small 2012).

Using the spatial PC maps for dimensions 2 and 3 as anomaly detectors makes it possible to immediately identify locations exhibiting strong seasonality or changes between 2019 and 2020. Figures 7–10 show contrasting examples of seasonal changes in 2020 as tri-temporal composites with example time series for individual pixels. Figures 11–13 show examples of interannual changes between 2014, 2017 and 2020 as tri-temporal composites from april of each year. April was chosen to minimize the effects of snow cover at higher latitudes and elevations and the effect of monsoon cloud cover at lower latitudes.