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

Abstract

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.

1. 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 agglomeration has stimulated rapid economic growth, this process has brought many environmental problems, including green space occupation, water scarcity, water and air pollution, and the urban heat island effects, particularly for developing countries [6,7,8,9,10].

China has been experiencing intensive urbanization since the 1980s, spatial forms of integrated cities are also noticed, for example, the Beijing–Tianjin–Hebei (BTH), the Yangtze River Delta (YRD), and the Pearl River Delta (PRD) regions [11,12,13,14,15]. Spatial and temporal heterogeneities of urbanization were first noticed by the Chinese population geographer Hu Huanyong in 1935 who proposed the “geo-demographic demarcation line” [16] separating the eastern and western parts of the country. The eastern part covers 36% of China, but accounts for 96% of the Chinese population. After 2000, the population in the west grew faster than in the east, but the growth did not change the situation. Urbanization shares a similar geographical distribution as population showing a fast development in the eastern and a slow development in the western regions [17,18].

Understanding the spatiotemporal dynamics of the urban agglomeration is crucial to optimize urban land use patterns, to support regional urban development strategies, and thus to address environmental problems related to fast changes in urban land cover [19]. Nighttime light from the Suomi National Polar orbiting Partnership visible infrared imaging radiometer suite (NPP/VIIRS) offers a new data source for extracting urban information. Efforts have been made to rapidly extract relatively accurate urban areas such as built-up [20] or impervious surface areas [21], and to reveal the responses of nighttime light intensity to socioeconomic development, for example, urban land expansion [22], human population [23], electric power consumption [24], gross domestic product and other socioeconomic variables [25,26,27].

The reported nighttime light based spatiotemporal heterogeneity of urban agglomeration in China is analyzed using data either from a single satellite [28] or from assimilation of more than one satellite [29,30,31]. It considers spatiotemporal heterogeneity by comparing (1) the three major urban agglomeration: Beijing–Tianjin–Hebei, Yangtze River Delta and Pearl River Delta [11,32]; (2) cities from eastern coastal region with cities from western mountainous region [33,34]; and (3) city core with the surrounding suburbs [35]. In applications, NPP/VIIRS nighttime light data are of much better quality than the nighttime light data from the Defense Meteorological Satellite Program Operational Linescan System (DMSP/OLS), but these two datasets suffer from the drawbacks of blooming effects from the near-infrared band [23,36,37]. Furthermore, DMSP/OLS images are inherently blurred, due to the footprint expanding away from the nadir, the limited data storage and the geo-location errors [38].

In summarizing, combined with nighttime light from the DMSP/OLS and NPP/VIIRS, 20 years of nighttime light data are provided as a response to the increasing demand for analyzing the effects of the changes of built-up areas. However, to our knowledge, few studies have been concerned with the data assimilation procedures associated with a time changing threshold to identify the built-up areas from DMSP/OLS (data before 2013) and NPP/VIIRS (data after 2012), although this is necessary due to considerable inconsistencies, different spatiotemporal quality, drawbacks of blooming effects and radiometric heterogeneities; in this paper, a time changing threshold is introduced to identify the built-up areas. Thus, to obtain a consistent analysis, it is urgently needed to adopt such novel approaches in order to calibrate DMSP/OLS (1992–2013) and NPP/VIIRS (2012 to present) nighttime light datasets to monitor urban extension and the kind of its evolution in the last 20 years. After introducing the methodology (Section 2), applications to Yangtze River Delta show the spatiotemporal variabilities of urban agglomeration in terms of urban area cover, nighttime light intensity, tendency shifts, and hotspot regions (Section 3 and 4), which is followed by the conclusion (Section 5).

2. Materials and Methods

The study area (Figure 1a), Yangtze River Delta, covers 3.7% of the Chinese land but with 1/6 of Chinese total population and 1/4 of the total gross domestic product (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.

2.1. 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.

2.2. Background Noise Mask

Three background noise mask methods were selected for defining the upper/lower limits to exclude background noise, gas, and fire cases (see Figure 2a–c): (a) an absolute lower limit (0.3 × 10⁻⁹ W·cm⁻²·sr⁻¹) and a relative upper limit (the city core) (referred as T_0.3 hereafter; see [24,30]; (b) an absolute lower limit (3 × 10⁻⁹ W·cm⁻²·sr⁻¹) and a relative upper limit (the city core) (referred as T₃ hereafter; see [40]); and (c) the limits from the NPP/VIIRS annual composite of 2015 (referred as T₂₀₁₅ hereafter; see [41,42]).

2.3. 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 sine (IHS)) transformation is performed to adjust the range of the NPP/VIIRS data to be comparable with the DMSP/OLS data.

$$\begin{array}{c}NPP/VIIR{S}_{New}=ln\left(NPP/VIIR{S}_{>T}+1\right) \\ or NPP/VIIR{S}_{New}=IHS\left(NPP/VIIR{S}_{>T}\right)\\ =ln\left(NPP/VIIR{S}_{>T}+\sqrt{NPP/VIIR{S}_{>T}^2+1}\right)\end{array}$$

(1)

Here $\mathrm{NPP}/{\mathrm{VIIRS}}_{>\mathrm{T}}$ ($\mathrm{NPP}/{\mathrm{VIIRS}}_{\mathrm{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).

$$f\left(x\right)=a+\frac{b}{1+e^{-c\left(\log NPP/VIIR{S}_{New}-d\right)}}$$

(2)

2.4. 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² = 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).

2.5. 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).

$$tendency\left(x\right)=\frac{n\ast {{\displaystyle \sum }}_{yr=1}^{n}yr\ast x_{yr}-\left({{\displaystyle \sum }}_{yr=1}^{n}yr\right)\left({{\displaystyle \sum }}_{yr=1}^n{x}_{yr}\right)}{n\ast {{\displaystyle \sum }}_{yr=1}^{n}y{r}^2-{\left({{\displaystyle \sum }}_{yr=1}^{n}yr\right)}^2}$$

(3)

2.6. 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:

$$G{i}_i^{*}\left(Z\right)=\frac{{{\displaystyle \sum }}_{j=1}^n{w}_{i,j}{x}_j-\overline{X}{{\displaystyle \sum }}_{j=1}^n{w}_{i,j}}{S\sqrt{\frac{\left[n{{\displaystyle \sum }}_{j=1}^n{w}_{i,j}^2-{\left({{\displaystyle \sum }}_{j=1}^n{w}_{i,j}\right)}^2\right]}{n-1}}}$$

(4)

$$with w_{i,j}=\left\{\begin{array}{c}1, distance to x_j<100 \mathrm{km} or 3 \mathrm{km} \\ 0, distance to x_j>100 \mathrm{km} or 3 \mathrm{km}\end{array}\right.$$

and $\overline{X}={{\displaystyle \sum }}_{j=1}^n{x}_j/n$ and $S=\sqrt{{{\displaystyle \sum }}_{j=1}^n{x}_j^2/n-{\left(\overline{X}\right)}^2}$.

$G{i}_i^{*} \left(Z\right)$ 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 $\overline{X}$ is the averaged magnitude over all the cities or pixels n. $w_{i,j}$ is the spatial weight 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).

2.7. 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\left(t_0\right) \le$ 0.05 were considered to be statistically significant change-points in a given time series.

$$K_{t,N}=max\left|{{\displaystyle \sum }}_{i=1}^t{{\displaystyle \sum }}_{j=t+1}^{N}sgn\left(x_i-x_j\right)\right| , 1\le t<N$$

(5)

$$P\left(t_0\right)=2exp\left\{-6K_{t,N}^2/(N^3+N^2)\right\}$$

(6)

where $sgn\left(x_i-x_j\right)=\left\{\begin{array}{c}1, x_i-x_j>0\\ 0, x_i-x_j=0\\ -1, x_i-x_j<0\end{array}\right.$.

3. Results

3.1. Background Noise Mask and Data Assimilation

Areas of three background noise masks indicate T_0.3 > T₂₀₁₅ > T₃ for the Yangtze River delta. Gibson [42] indicates that using the NPP/VIIRS annual composite for 2015 (T₂₀₁₅) 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² value of 0.9697 (see Figure 2). Thus, the logarithmic transformation was used in the following analysis.

3.2. Threshold Setting for the Built-Up Areas

The nighttime light based built-up areas and their thresholds in Suzhou from 2001 to 2019 were extracted with the methods presented in Section 2 based on the manual interpretation of the Landsat images in 2002, 2006, 2010, 2014, and 2018 (see Figure 3). Accuracy test between nighttime light based built-up areas and Landsat images interpreted areas in the years 2013 and 2017 shows overall accuracies of 92.99% and 92.19% with Kappa coefficients of 0.67 and 0.66, respectively. That is, the selected thresholds for built-up regions extraction in the present study have credibility. Therefore, the spatiotemporal evolution of the built-up regions in Yangtze River Delta is extracted based on the threshold from Figure 3d. The following results are noticed (see also Figure 3e):

(1) Nighttime light intensity thresholds for DMSP/OLS indicate a systematic shift from 1 per year between 2000 and 2012 (thresholds change from 44 to 57), while NPP/VIIRS data are much more stable with the threshold value 58 after 2013. (2) Area expansion (averaged nighttime light intensity) increased by four to five times from 5000 km² in 2001 to 24,000 km² in 2019 (from 5.18 in 2001 to 21.2 in 2019). (3) Two stage transforms from moderate to fast developments of the built-up areas and the averaged nighttime light intensity were divided around 2010, and further detailed explanations about the most frequent shifted years with Pettitt’s test see Section 4.3.

3.3. Accuracy Test with MODIS Data

Shi et al. [20] reported that for 12 Chinese cities, thresholds vary from 2.9 to 23.7 $\times$ 10⁻⁹W·cm⁻²·sr⁻¹ (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. Accuracy test comparing nighttime light 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.

	2005	2010	2015
Overall accuracy (%)	97.56	97.29	95.49
Kappa coefficient	0.66	0.68	0.65

) 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.

4. 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).

In built-up area detection (Figure 4a), the possible drawbacks from satellite, sensor, and data assimilation are excluded because, instead of using a constant nighttime light threshold, the optimal thresholds are determined in 2002, 2006, 2010, 2014, and 2018 when area from DMSP/OLS and NPP/VIIRS nighttime light data were closest to the area from the manual interpretation of the Landsat images. Accuracy tests in 2013 and 2017 with built-up areas from the manual interpretation of the Landsat images show overall accuracies of 92.99% and 92.19% with Kappa coefficients of 0.67 and 0.66. In addition, the area tendency ratios of NPP/VIIRS to DMSP/OLS with versus without data assimilation were similar and about 3 (1504/481 km²/year = 3.12 versus 1623/460 km²/year = 3.52).

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 (2001–2010 versus 2011–2019) were analyzed over the Yangtze River Delta, with built-up area tendency on administrative city scale (Section 4.1), nighttime light intensity tendency on pixel scale (Section 4.2), and in spatiotemporal joint analysis combining tendencies of built-up area and nighttime light intensity (Section 4.3).

4.1. 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²/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²/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].

4.2. 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 (2011–2019 in Figure 6b), rapid (moderate) increasing tendencies of the nighttime light intensity are concentrated as a cluster around Suzhou surrounded by Wuxi, Jiaxing, Changzhou, and Shanghai, with averaged nighttime light intensity tendency values of 2.50 > 2.15 >1.87 > 1.84 > 1.61 (1.04, 1.06, 1.78, 1.48, 0.94) per year. In contrast, inner cities in Zhejiang province, Quzhou, Huzhou, Jinhua, and Zhoushan have been undergoing a “moderate to rapid” transition increasing from 0.58, 0.92, 0.93, 1.23–1.80, 1.85, 1.80, and 2.08 per year in terms of the averaged nighttime light intensity tendencies (2001–2010 versus 2011–2019). Similarly, inner cities in Anhui province, Chuzhou, Anqin, Lu’an, and Bengbu have increased from 0.60, 0.68, 0.68, 0.69–1.74, 1.64, 1.84, and 1.62 per year.

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.

4.3. Spatial-Temporal Joint Analysis

Two periods’ hot spot and cold spot analyses (2001–2010 versus 2011–2019) over the Yangtze River Delta with built-up area tendency in administrative city scale (Section 4.1) and nighttime light intensity tendency in pixel scale (Section 4.2) list cities with high and low tendency values. To further diagnose (1) the tendency shift in the time series and (2) the total gross domestic product (GDP) associated spatiotemporal variability in the 41 administrative cities over Yangtze River Delta. Pettitt’s test [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.

4.3.1. 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 Figure 5d and Figure 6d) are highlighted in bold (Figure 7a,b) emphasizing both a southward and an inner cities’ development.

4.3.2. 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:

(1) The fast-developing cities show an increase (decrease) of the built-up area (nighttime light intensity) tendency from the first to the second period, with 0.14–0.47 to 0.57–0.85 (0.40–0.77 to 0.08–0.47). (2) The moderate- and slow-developing cities increase both in the built-up area and the nighttime light intensity tendency, with 0.03–0.08 to 0.29–0.51 and 0–0.06 to 0.08–0.19 in terms of built-up area tendency, 0.25–0.34 to 0.41–0.75 and 0.11–0.21 to 0.40–0.67 in terms of nighttime time light intensity tendency. (3) Comparing 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.

5. 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.