Spatiotemporal Heterogeneity Analysis of Yangtze River Delta Urban Agglomeration: Evidence from Nighttime Light Data (2001–2019)

The long-term changes of the relationship between nighttime light and urbanization related built-up areas are explored using nighttime light data obtained from the Defense Meteorological Satellite Program/Operational Linescan System (DMSP/OLS, data before 2013) and the Suomi National Polar-orbiting Partnership Visible Infrared Imaging Radiometer Suite (NPP/VIIRS, data after 2012) and information of the spatiotemporal heterogeneity of urban evolution. This study assimilates two datasets and diagnoses the spatial heterogeneity in administrative city scale based on built-up area tendencies, temporal heterogeneity in pixel scale based on nighttime light intensity tendencies, and GDP associated spatiotemporal variability over the Yangtze River Delta comparing the first two decades of this century (2001–2010 versus 2011–2019). The analysis reveals the following main results: (1) The built-up areas have generally increased in the second period with the center of fast expansion moving southward, including Suzhou-Wuxi-Changzhou, Hangzhou, Ningbo, Nanjing, and Hefei. (2) Urban development in the original city core has saturated and is spilling over to the suburbs and countryside, leading to nighttime light intensity tendency shift from a “rapid to moderate” and a “moderate to rapid” development (a “hot to cold” and a “cold to hot” spatial clustering distribution). (3) The tendency shifts of built-up area and nighttime light intensity occur most frequently in 2010, after which the urban development is transforming from light intensity growth to built-up area growth, particularly in the developed city cores. The urban agglomeration process with nighttime light intensity reaching saturation prior to the urban development spreading into the surrounding suburbs and countryside, appears to be a suitable model, which provides insights in addressing related environmental problems and contribute to regional sustainable urban planning and management.


Introduction
Urban agglomeration is a highly developed spatial form of integrated cities, which can be defined as a contiguously built-up area, shaped by one core city or by several adjacent cities, sharing industry-, infrastructure-, and housing-land use with high-density levels with embedded open spaces [1][2][3]. It has been witnessed worldwide in recent decades as one of the most significant changes in human society, and as a "dynamic, multidimensional socio-spatial process" [4]. Until now, 55% of the world's population lives in the urban areas, and it is predicted that this number will reach 68% by 2050 [5]. Although urban (according to the 2020 National Statistics Yearbook of China). Suitably employing the tools introduced in the previous section ( Figure 1b) provides a scheme for analyzing urban agglomeration evolution.

Data and Preprocessing
The built-up areas of Suzhou in every four years commencing with 2002 are obtained from the manual interpretation of the Landsat images, which are used as the reference data (30 m × 30 m) to determine the optimal threshold for the DMSP/OLS and NPP/VIIRS nighttime light data based built-up areas. The version 4 cloud-free DMSP/OLS nighttime light data is an annual composite from 2001 to 2013 with a spatial resolution of 30 arc-seconds (1 km), which is downloadable from the Nation Oceanic and Atmospheric Administration's National Geophysical Data Center (NOAA/NGDC). These datasets (1992-2013) were obtained from six satellites (F10: 1992-1994, F12: 1994-1999, F14: 1997-2003, F15: 2000-2007, F16: 2004-2009, and F18: 2010-2013) lacking on-board calibration. To be comparable, the invariant region method [39] based intersatellite calibration from Zhang et al. [32] are used. The NPP/VIIRS nighttime light monthly composite from the NOAA/NGDC is available since 2012 (https://ngdc.noaa.gov/eog/index.html, accessed on 3 July 2020) with a spatial resolution of 15 arc-seconds (0.5 km). Annual means are calculated by averaging data from January to April and August to December (as summer data are missing in the high latitudes of China). Both DMSP/OLS and NPP/VIIRS nighttime light images are reprojected to the Albers Equal Area Conic projection. The NPP/VIIRS nighttime light data was resampled according to the spatial resolution of DMSP/OLS data (1 km × 1 km) by the nearest neighbor method. Geometric correction was applied to DMSP/OLS data by manual interpretation comparing with the NPP/VIIRS data. ArcGIS and ENVI software were used.

Figure 2.
Results from three background noise mask methods for defining the upper/lower limits: (a) an absolute lower limit (0.3 × 10 −9 W·cm −2 ·sr −1 ) and a relative upper limit (the city core); (b) an absolute lower limit (3 × 10 −9 W·cm −2 ·sr −1 ) and a relative upper limit (the city core); and (c) the limits from the NPP/VIIRS annual composite of 2015. Additionally, the corresponding sigmoid regression performance after applying the logarithmic (the middle row, d-f) and IHS transformation (the bottom row, g-i) to the masked NPP/VIIRS data in 2013.

Data Assimilation
To obtain a consistent long nighttime light series, the approach developed by Zhao et al. [30] (or Gibson et al. [43]) is employed: first, a logarithmic (or inverse hyperbolic Remote Sens. 2021, 13, 1235 5 of 17 sine (IHS)) transformation is performed to adjust the range of the NPP/VIIRS data to be comparable with the DMSP/OLS data.
NPP/V I IRS New = ln(NPP/V I IRS >T + 1) or NPP/V I IRS New = I HS(NPP/V I IRS >T ) Here NPP/VIIRS >T (NPP/VIIRS New ) is the original (adjusted) NPP/VIIRS radiance value; the constant 1 ensure data consistency with transformed NPP/VIIRS nighttime light data greater than 0. T is the threshold for excluding background noise, gas, and fire cases. Then, the overlapping data (2013) from DMSP/OLS and NPP/VIIRS were selected to fit the coefficients a, b, c, and d for a sigmoid function comparison (see Figure 2).

Threshold Setting for the Built-up Area
According to the method proposed by Small et al. [44] and Liu et al. [45], the optimal threshold is determined in 2002,2006,2010,2014, and 2018 when area from DMSP/OLS and NPP/VIIRS nighttime light data closest to the area from the manual interpretation of the Landsat images (see Figure 3a,b). A quadratic regression with R 2 = 0.997 was applied to set the threshold values for the remaining years ( Figure 3d). The quadratic regression-based threshold values of 2013 and 2017 were selected and tested with built-up areas from the manual interpretation of the Landsat images (detailed accuracies see Section 3).

Area Tendency and Light Intensity Tendency
To estimate the long-term urban agglomeration development, a linear regression analysis (against year from yr = 1, 2, . . . n) was employed to determine x, the area tendency and light intensity tendency from the time development on city scale and pixel scale. The results are subjected to an F-test with all significance levels on the p-value < 0.05 level (unless noted otherwise).

Hot and Cold Spot Analysis
The Getis-Ord Gi* statistic has been widely adopted to analyze biological habitat [46], epidemic disease [47], and crime [48] to identify the significant spatial clusters of high values (hot spots) and low values (cold spots), which characterize the intensity of the geographic relation between observed data within a region. This functionality has been included in the automatic output of the hot spot analysis (Getis-Ord Gi*) tool in ArcGIS software [49], and is calculated as follows: with w i,j = 1, distance to x j < 100 km or 3 km 0, distance to x j > 100 km or 3 km and X = ∑ n j=1 x j /n and S = ∑ n j=1 x 2 j /n − X 2 .
Gi * i (Z) is the Getis-Ord Gi* statistical z-score value for an administrative city or a pixel i, which describes the spatial dependency of i on all the surrounding administrative cities or pixels j. x j represents the magnitude of area tendency or light intensity tendency at location j, and X is the averaged magnitude over all the cities or pixels n. w i,j is the spatial weight Remote Sens. 2021, 13, 1235 7 of 17 value between the administrative city or pixel i and j, whose value equals either 1 or 0 defined by the distance between i and j. In this study, distances of 100 km and 3 km were selected for city scale and pixel scale, respectively. Clusters of high values (hot spots) and low values (cold spots) passing the 90% confidence level were categorized from the long-term tendencies of nighttime light intensity and area coverage. For a detailed method description see [50] and https://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-statistics-toolbox/ h-how-hot-spot-analysis-getis-ord-gi-spatial-stati.htm (accessed on 11 January 2021).

Tendency Shift
Pettitt's test is a non-parametric test for detecting the tendency shift in a time series [51].
Here it was applied to the time series of the built-up areas and nighttime light intensities with the null hypothesis that there is no tendency shift or change-point in a time series. The test statistic K t,N and the associated probability P are given. Years in the time series with P(t 0 ) ≤ 0.05 were considered to be statistically significant change-points in a given time series. where

Background Noise Mask and Data Assimilation
Areas of three background noise masks indicate T 0.3 > T 2015 > T 3 for the Yangtze River delta. Gibson [42] indicates that using the NPP/VIIRS annual composite for 2015 (T 2015 ) as the background noise mask has the advantage of creating a like-with-like series to link to DMSP/OLS (because the NPP/VIIRS annual composite was processed to the same standard as the DMSP/OLS composite). However, when applying the sigmoid regression after the logarithmic and IHS transformation to the masked NPP/VIIRS data in 2013 ( Figure 2), the results indicate that ln(NPP/VIIRS > 0.3 + 1) performs with a higher R 2 value of 0.9697 (see Figure 2). Thus, the logarithmic transformation was used in the following analysis.

Accuracy Test with MODIS Data
Shi et al. [20] reported that for 12 Chinese cities, thresholds vary from 2.9 to 23.7 × 10 −9 W·cm −2 ·sr −1 (or DN from 22 to 61). Thus, a further accuracy test comparing our threshold based built-up area with the urban built-up area from MODIS Land Cover products (MODIS 12Q1, 500 m spatial resolution) in 2005, 2010, and 2015 (see Table 1) was calculated, to test whether the threshold of Suzhou was representative for deriving urban extent for all the cities over the Yangtze River Delta. The results show overall accuracy reaching 97.56%, 97.29%, and 95.49% with Kappa coefficients of 0.66, 0.68, and 0.65, respectively.

Analyses and Discussion
DMSP/OLS data F18 2010-2013 are considered problematic, because it records much higher DN values than earlier DMSP/OLS satellites [52]. Results from splitting the data into two decades (that is, 2001-2010 and 2011-2019) might be affected rather by different sources of nighttime light data than by inherent urban and nighttime light change tendencies. A sensitivity test, to exclude the possible lacking assurance from (1) satellite differences, (2) sensor drawbacks, and (3) data assimilation, is designed by calculating the original (excluding F18 and without data assimilation) area and/or nighttime light intensity tendencies of the DMSP/OLS data in 2001-2009 versus the NPP/VIIRS data in 2012-2019 (also without F18 and data assimilation; Figure 4). The NPP/VIIRS nighttime light intensity tendencies with/without data assimilation are compared (2012-2019) and presented as a frequency (number density) distribution in Figure 4b. Despite regional deviations, the NPP/VIIRS light intensity tendencies with/without data assimilation over the Yangtze river delta were almost linear. That is, the possible drawbacks from satellite, sensor, and data assimilation show less sensitivity in the following tendency analysis. Therefore, two periods

Spatial Heterogeneity in the Administrative City Scale
The built-up areas in each city of Yangtze River Delta from 2001 to 2019 are extracted with the threshold presented in Figure 3d. A linear regression is applied to the built-up area time series of each city and the tendencies with significance p-value < 0.05 are used for a further hot and cold spot analysis ( Figure 5).
In the first period of 2001-2010 (Figure 5a), a rapidly increasing tendency of the built-up area could be diagnosed as a scattered distribution center around the core city of Shanghai, from a fast to a slow development Suzhou > Shanghai > Wuxi > Changzhou with the associated tendency values 126.94 > 66.84 > 42.53 > 32.29 km 2 /year, respectively. While cities far away from Shanghai like Anqing, Chizhou, and Xuanchen in the Southern Anhui province, and Quzhou, Wenzhou, and Taizhou in the Southern Zhejiang province show no significant increasing tendencies of built-up area. In the second period of 2011-2019 (Figure 5b), the rapidly increasing tendencies clustering the surrounding cities into one metropolitan domain, which means that cities close to Shanghai and along the costal line undergo a rapid built-up area expansion, with Ningbo > Suzhou > Hangzhou > Shanghai revealing tendency values 121.81 > 109.05 > 91.86 > 87.69 km 2 /year. Inner land cities with a relatively long distance to Shanghai show lower tendency values of urban spatial expansion, that is Huangshan < Chizhou < Tongling < Huainan < Huaibei with the tendency values 0.78 < 3.26 < 7.33 < 9.01 < 9.36 in Anhui Province. All the findings are in agreement with "the Riverside Development Belt", "the Shanghai-Hangzhou-Jinhua development belt", and "the coastal development belt" development plan (2015-2030) [53].
Hot spot and cold spot distributions for the two periods highlight the significant spatial clusters of high area tendency values (hot spots) and low area tendency values (cold spots). In the first period of 2001-2010 (Figure 5c), the hotspots of built-up area tendencies are clustered around the core city Shanghai (as a concentration distribution), surrounded by Suzhou, Wuxi, Changzhou, Nantong, and Jiaxing, which is consistent with a previous study [54]. By contrast, in the second period of 2011-2019 (Figure 5d), the hotspot areas were no longer clustered with the core city Shanghai but spreading to the southern regions of the Yangtze River Delta. Hangzhou, Shaoxing, and Jinhua from Zhejiang province replacing Wuxi and Changzhou bloom in terms of the built-up area. That is the spatial form of integrated cities from several provinces (or urban agglomeration) was noticed (see also Figure 5b and [55,56]). Cold spots locate on Bengbu, Chizhou, and Huangshan, Anhui province, which may be related to floods, population exodus, and low industrialization in Northern Anhui [57], and mountainous area limitation and natural landscape protection (Huangshan National Forest Park) in Southern Anhui [58].

Temporal Heterogeneity in the Pixel Scale
As nighttime light intensity is closely related to electricity consumption, it has been used as an indicator of the growth of built-up areas. The linear nighttime light intensity tendencies of significance (p-value < 0.05) over Yangtze River Delta were calculated for a further hot and cold spot analysis ( Figure 6) and show the following results: From the first to the second period, spatial heterogeneity of a "rapid to moderate" and a "moderate to rapid" increasing tendency could be diagnosed. In the period of 2001-2010 in Figure 6a  Hot spot and cold spot distributions for the two periods highlight the significant spatial clusters of high nighttime light intensity tendency values (hot spots) and low nighttime light intensity tendency values (cold spots). Heterogeneity of a "hot to cold" and a "cold to hot" spatial distribution could be diagnosed. In the first period of 2001-2010 (Figure 6c), the statistically significant hotspots of nighttime light intensity tendency are mainly concentrated around Suzhou with a "Z"-shape distribution as mentioned in previous studies [56,59], surrounded by cold spots in its suburbs and the countryside. In contrast, for the period of 2011-2019 (Figure 6d), the statistically significant cold spots of the nighttime light intensity tendency occupied the "Z"-shape, including Shanghai, Suzhou, Wuxi, Changzhou, Hangzhou, and Ningbo.  [51] is applied to the time series of the built-up areas and nighttime light intensities, and then 41 administrative cities are categorized based on local total gross domestic product, into fast-developing cities, moderate-developing cities and slow-developing cities. Jenks Natural Breaks Classification Method [60] is applied to make sure within category are more similar than the between category.

Tendency Shift Detection
In the urbanization process of southward development and several inner cities' booming, 22 (29) cities undergo statistically significant tendency shifts in terms of the annual built-up area (nighttime light intensity) (Figure 7). Years with statistically significant tendency shifts occur most frequently in 2010: in total, 19/41 cities show the tendency shift for built-up area (8/41 cities) or nighttime light intensity (13/41 cities). In the Yangtze River Delta, the cities with statistically significant tendency shift belong to the hot spot cities (see Figures 5d and 6d) are highlighted in bold (Figure 7a,b) emphasizing both a southward and an inner cities' development.

GDP Associated Spatiotemporal Variability
Gross domestic product (GDP) is an indicator of economic activity, which closely relates to the development level of a region and has a strong correlation with nighttime light intensity [61]. Thus, fast-, moderate-, and slow-developing cities were categorized according to the total GDP (2001-2018) of the 41 cities in Yangtze River Delta (Figure 8). The normalized statistics provide the following information: fast-, moderate-, and slow-developing cities in the first period, the increases of nighttime light intensity tendency all exceed the increases of built-up area. In the second period, the increases of the built-up area tendency in moderate-and slow-developing cities are still less than the increases of the nighttime light intensity tendency. However, the increases of the built-up area tendency in fast-developing cities overwhelm the increases of the nighttime light intensity tendency due to the nighttime light saturation in the city cores and a fast expansion to suburbs and countryside. In summary, in the first period of their evolution, fast-developing cities were dominated by increasing nighttime light intensity, while in their second period this increase was shifted to an increase of the built-up area. The increases of the built-up area tendencies were positively related to the GDP, showing fast > moderate > slow city development, and the growth of built-up areas lagged behind the growth of light intensity. The "fast to slow" transition of the nighttime light intensity tendency indicates whether the city core was fully developed. In this sense, moderate-and slow-developing cities show potential for a "fast to slow" conversion of the nighttime light intensity tendency in the future. In the second period, the increases of the built-up area tendency in moderate-and slow-developing cities are still less than the increases of the nighttime light intensity tendency. However, the increases of the built-up area tendency in fast-developing cities overwhelm the increases of the nighttime light intensity tendency due to the nighttime light saturation in the city cores and a fast expansion to suburbs and countryside.
In summary, in the first period of their evolution, fast-developing cities were dominated by increasing nighttime light intensity, while in their second period this increase was shifted to an increase of the built-up area. The increases of the built-up area tendencies were positively related to the GDP, showing fast > moderate > slow city development, and the growth of built-up areas lagged behind the growth of light intensity. The "fast to slow" transition of the nighttime light intensity tendency indicates whether the city core was fully developed. In this sense, moderate-and slow-developing cities show potential for a "fast to slow" conversion of the nighttime light intensity tendency in the future.

Conclusions
The global urbanization process, where human population, material demands of production, human consumption, and urban waste discharge expand and become more intense, has recently emerged as a sustainability challenge in an increasingly urbanized world, particularly for more rapidly changing developing countries [62]. Nighttime light signals derived from the Defense Meteorological Satellite Program's Operational Line scan System (DMSP/OLS) first provide striking remotely sensed data to analyze spatiotemporal changes in global urbanization processes [63]. Then, the Suomi National Polar orbiting Partnership Visible Infrared Imaging Radiometer Suite (NPP/VIIRS) nighttime light data of much higher quality than the nighttime light data from the DMSP/OLS are used to quantify the relationship between night lighted areas and urbanization variables, including urban expansion [22], population growth [64], electric power consumption [24], gross domestic product [42], and socioeconomic parameters [23,25]. However, to our knowledge, few studies have been concerned with the data assimilation from DMSP/OLS (data before 2013) to NPP/VIIRS (data after 2012) associated with a time changing threshold to define the built-up areas, but which is necessary due to considerable inconsistencies, different spatiotemporal quality, drawbacks of blooming effects and radiometric heterogeneities.
Thus, research on long-term changes of urbanization over rapidly changing Yangtze River Delta may provide insights in diagnosing how long-time spatiotemporal heterogeneities of the urban evolution so as to provide insights in predicting a further urban agglomeration on a regional and global scale. This study compares dynamics of urbanization in 2001-2010 and 2011-2019 in the city unit scale and pixel unit scale, and the main conclusions are as follows.
(1) Comparing the two periods of area expansion indicates a general increase of the built-up areas in the second period of 2011-2019, and its geographical development gradually moves southward in the Yangtze River Delta. The hierarchy of the urbanization process of focusing on the core city is less obvious in the second period. Suzhou-Wuxi-Changzhou, Hangzhou, Ningbo, Nanjing, and Hefei contribute a fast development in the urban agglomeration, which is in agreement with the development plans for the Yangtze River Delta (2015-2030) to construct a networked of spatial pattern distribution characterized by "one mega, five circles and four belts" [53].
(2) The development transitions of a "rapid to moderate" and a "moderate to rapid" and the "hot to cold" and a "cold to hot" distribution indicate a continuous urban expansion. That is, regional development is first associated with a rapid increasing light intensity in the core cities, until the saturation of the light intensity. Then the rapidly increasing light intensity spills over, and thus accelerates the development of the surrounding suburbs and countryside in terms of both built-up area and light intensity. Likewise, the development of the surrounding suburbs and countryside appear to repeat this kind of development process "rapid to moderate" or "hot to cold", which spills over (once more) when the nighttime light intensity saturated.
(3) The tendency shift of the built-up areas and nighttime light intensities for the 41 administrative cities over Yangtze River Delta occurs mainly in the year 2010, which indicates a decade as a characteristic time scale of urban evolution. Thereafter, the urban development characteristic has changed from a centralized pattern to a scattered distribution and undergoes an urban agglomeration process. In the process, the development time sequence is first the light intensity growth and then built-up area growth. In a developed (developing) city core, light intensity tendency is less (greater) than the built-up area tendency.
In general, urbanization is one of the extreme processes in land use and land cover change. Urbanization induced change of in surface parameters, including increased atmospheric greenhouse gas concentration, aerosol emissions, and land use and land cover change, is a major component for the formation and evolution of urban climates [65][66][67]. It accounts for effects like the urban heat island [68,69], affecting not only local and regional climates, but also water resources, air quality, human health, and biodiversity and ecosystem functioning [70]. Thus, research on the long-time evolution process of global/regional urbanization may provide insights in addressing the environmental problems arising from long-term and rapid urbanization and strengthening global or regional sustainable urban planning and management.