Evaluating the Effects of Urbanization Evolution on Air Temperature Trends Using Nightlight Satellite Data

: Confounding factors like urbanization and land-use change could introduce uncertainty to the estimation of global temperature trends related to climate change. In this work, we introduce a new way to investigate the nexus between temporal trends of temperature and urbanization data at the global scale in the period from 1992 to 2013. We analyze air temperature data recorded from more than 5000 weather stations worldwide and nightlight satellite measurements as a proxy for urbanization. By means of a range of statistical methods, our results quantify and outline that the temporal evolution of urbanization affects temperature trends at multiple spatial scales with signiﬁcant differences at regional and continental scales. A statistically signiﬁcant agreement in temperature and nightlight trends is detected, especially in low and middle-income regions, where urbanization is rapidly growing. Conversely, in continents such as Europe and North America, increases in temperature trends are typically detected along with non-signiﬁcant nightlight trends. thermal impact, by relating the temperature variations in recent years to the corresponding variations in nightlight data. We investigate the presence of possible bias in temperature


Introduction
The urban transition leads to alterations in landscape conditions and to important modifications in the urban climate, along with several environmental problems e.g., on water use and quality, on the generation of air pollution, and on the production of solid waste and sewage [1,2]. Because of rapid urbanization growth, even more considerable impacts are expected on a broader scale and especially in developing countries, like higher consumption of energy, goods, services, and resources demand, which have the potential for greater negative impacts on global environments and ecosystems [3][4][5][6][7]. Changes in food supplies, freshwater resources, and increase in extreme weather events (e.g., heatwaves and droughts) are expected to lead to several consequences on human health in terms of e.g., heat stress, cardio-respiratory, and infectious diseases [2,8]. In this regard, we must also consider that 55% of the world population is residing in urban areas in 2018, which is projected to reach 68% by 2050 [9]. In the context of global climate change, it is crucial to better investigate how urban growth affects temperature record trends to consistently attribute the causes of observed warming at wider scales [10][11][12][13][14][15][16].
The impact of urbanization on near surface temperatures has been investigated since the 1980s [17,18]. These studies suggested that a proportion of global warming observed on the last century timescale could be related to local warming induced by urbanization. Rapid urban growth has resulted in the expansion of built-up areas in and around cities, particularly for nations and regions experiencing demographic expansion. This plays a crucial role in the near-surface warming and on temperature Recently, several studies were conducted in order to account and compensate for the effect of urban warming when estimating large scale temperature trends [11,14,15]. Global warming trend analyses are conducted, for example, using adjusted urban temperature data, by comparing observations to the sea surface temperature [11], comparing temperature values in different weather conditions [38], or removing sites with suspected urban warming from global [11] and regional warming analyses [39]. However, very often, the poor spatial coverage of the weather stations network does not guarantee the possibility to compare urban stations with the corresponding rural ones [40]. Most weather station networks are located in heavily anthropized sites [13,15]. Rapid urbanization makes the classification between urban and rural areas difficult in some cases [34].
New techniques and approaches could be helpful in detecting relationships and feedbacks between air temperature and land-use changes. These efforts result in independent and relatively accurate estimates of urban shape and extent. Remote-sensing products, which are primarily MODIS-500 maps [41][42][43][44], Landsat data [45], satellite night time images [3,46], and Synthetic Aperture Radar (SAR) data [47], are widely used to map global urban extent and to identify and geo-localize rural and urban stations [15,39]. In recent studies [46], day-night composites combine nightlights and Landsat to show consistencies in land cover and nightlights brightness. Other relevant works spatially assess the UHI signature on land surface temperature amplitude and its mutual relations with development size, intensity, and in different biomes, for over 3000 cities worldwide combining MODIS and night-time lights products i.e., a nightlight based impervious surface area (ISA) [48]. Nightlight ISA and Landsat ISA for cities in the continental US are first compared to bridge from previous studies. Results highlight significant positive relationships between UHI magnitude, ecological setting and ISA, and that nightlights are good estimators of urban sprawl and are more objective than methods based on population density. Zhou et al. [3] develop a method to map urban extent from nightlights on a global scale and compare it at the pixel level and regionally with five other widely used global urban products, including MODIS [43], GlobCover [49], and GLC2000 [50]. Results show that the nightlight map produced by the author is in high agreement with the other global urban area map products. Bagan and Yamagata [51] also show that the combined use of land cover data and nightlights are both good predictors of population density in Japan. On a local scale, recent studies [52] assess the properties of the temperature station network and the potential urban influence on temperature records by means of land cover, population density, and nightlights at an increasing buffer around the weather station. The three methods are found to be highly correlated, which indicates that nightlights and population density are good proxies of urban areas.
Thus, based on comparison with other existing global urban maps, night-time satellite images are demonstrated to be a good proxy of urban extent and allows for temporal dynamics analyses of urban areas [3,52]. It is shown that satellite nightlights maps (also abbreviated as "nightlights" from now on) allow tracking the spatiotemporal dynamics of the population much better than traditional census and administrative data [53,54]. The fine spatial resolution of nightlights (nearly 1 km at the equator) makes them a valuable proxy for the human presence, which is widely used in many research fields [55]. In the late 1990s, the scientific community started considering nightlights as a proxy for e.g., population density [56]. Many studies followed, focusing mainly on economic activity [31], electric power consumption patterns [32], or the level of poverty estimates [57]. Nightlight imagery are now largely employed to map and measure urban extents [3,56,58] and address environmental topics, e.g., light pollution in Europe [59]. Recently, nightlights data are being employed to relate flood risk to increasing human pressure near rivers on a global scale [60]. The same approach is adopted to study local human exposure to hydrogeological risk, e.g., mapping population exposure to natural hazards [61] and assessing the interaction between water resources and dynamic human systems [62].
In this study, we introduce a new approach to quantify the relationship between urbanization dynamics and thermal impact, by relating the temperature variations in recent years to the corresponding variations in nightlight data. We investigate the presence of possible bias in temperature records induced by land-use change-specifically the growth of urban areas-by means of night time satellite images as a proxy for urbanization.
The work is organized as follows. Nightlights and air temperature data are first repurposed and prepared for the analysis. Then, we assess the relationship between air temperature and nightlight evolution in the last 25 years, by means of trend estimation. We, thus, propose a series of ad hoc statistical techniques in order to shed light on the possible agreement between temperature and nightlight variations in time. Lastly, we discuss our findings and draw some conclusions.

Satellite Nightlights Time Series
Annual time series of nightlights are freely provided from the NOAA National Geophysical Data Center (NGDC) as satellite images, collected under the Defense Meteorological Satellite Program (DMSP), Operational Linescan System (OLS) [63]. OLS consists of two sensors operating, respectively, in the near-infrared (0.4 to 1.1 µm) and in the thermal infrared (10.5 to 12.6 µm) spectrums. Satellite images are collected on a yearly basis for a 22-year period from 1992 to 2013. Six satellites are used, with a total of 34 composite images, which generate a product called stable light [62]. Original satellite images are not onboard calibrated. Thus, before running any additional analysis, we inter-calibrate them, according to standard procedures [31,57]. In case of years presenting overlapping datasets, the averaged night time brightness value, extensively used in recent papers [60], is employed. Each pixel in the final product represents the stable light yearly average in a 6-bit format and it is expressed as a Digital Number i.e., a dimensionless numerical integer, which represents the brightness on the Earth's surface in each cell. Sensors could identify cloud-free night time lights from e.g., human settlements, fires, and gas flares. In this analysis, signals pertaining to fires and gas flares are removed, since they are not of interest in studying urbanization and human settlement dynamics. Digital Number values are proportional to radiance and range from 0 (pitch dark areas) to 63 (bright areas). When dealing with satellite imagery, possible discrepancies could be related to, for instance, light saturation and blooming effects, which potentially affects nightlight data. Light saturation and blooming occur primarily in highly built-up areas and developing countries, which are characterized by bright pixels surrounded by broad pitch-dark zones [31]. Saturation and blooming effects may lead to an overestimation of the extent of urban areas [64]. In recent years, machine-learning approaches are being used to map urban areas over broad scales, which account for a blooming effect and identify the boundaries of highly bright settlements [65]. Our analysis is based on calibrated images and adjusted nightlights values, which are not saturated at the highest intensities, in order to reduce saturation and blooming bias [60,62]. Moreover, since our analysis focuses on trends in nightlights (i.e., on the slope of the complete time series) rather than the actual level of urbanization (i.e., nightlights absolute values) [66], the potential blooming effect equally affects all years along the time series within the study period, with no significant impact on the overall trend.
Images are provided as raster products with a spatial resolution of 30 arc seconds (0.00833 • ), i.e., nearly 1 km resolution at the equator, with a spatial extension between 75 • N, 65 • S latitude, 180 • W, and 180 • E longitude.

Air Temperature Datasets
For the purpose of our analysis, nightlights imagery and temperature records should have the same spatiotemporal coverage and resolution. More specifically, (i) the coordinates of air temperature stations should have a spatial resolution of 30 arcsec, (ii) temperature data should be available from 1992 to 2013, and (iii) a global spatial coverage of temperature stations should be guaranteed, including in developing countries. Two temperature datasets are considered in our analysis, which are the Berkeley Earth [67] and the World Meteorological Organization (WMO) [68,69] dataset. More specifically, the Berkeley Earth dataset is used to analyze the temporal trends of temperature records since it is the best compromise between reasonable spatial coverage and temporal availability. The WMO dataset is used as a support to perform a quality control of the geographical position of weather stations.
The Berkeley Earth dataset combines data and metadata from 16 previous existing datasets [67]. The currently-available free dataset, after removing duplicate records, consists of 39000 stations. The dataset consists of three categories of data, as reported in Berkeley Earth [67]: (i) "Source Data", which include raw temperature data as reported from the original agencies, (ii) intermediate data, where all combined data from the original sources are filtered, merged together, quality-checked, and flagged, which result in the "Quality Controlled" dataset, (iii) output data i.e., the Breakpoint Adjusted Monthly Station data, which are post-processed to correct any discontinuity or heterogeneity and local systematic effects (e.g., UHI effect). The complete description of the data quality check processes can be found in Reference [67].
To avoid bias correction that may have been applied to compensate for UHI, we focus on the intermediate "Quality Controlled" product (that will be termed as "Quality Controlled" from now on), which has not been adjusted to resolve temporal discontinuities and long-term in homogeneities.

Data preparation
Data preprocessing entails the following steps: 1.
geo-localization of air temperature stations, 2.
pairing of temperature and nightlights data, 3.
gap-filling procedure for incomplete time series of temperature records.

Geo-Localization of Air Temperature Stations
The spatial resolution of temperature data included in the Berkeley Earth dataset is quite heterogeneous, with some stations localized with a 30 arcsec resolution, and others with a lower resolution. This occurs because this dataset includes station data gathered from different sources. To minimize uncertainties, metadata from the Berkeley Earth dataset are compared to official data provided by the WMO weather reports when available [68,69]. WMO weather reports include a full list of all surface and upper-air stations in use at a 30-arcsec resolution (0.00833 • ) and, thus, provides a precise geo-localization of temperature stations. An iterative statistical-based procedure is performed in order to compare station coordinates in both Berkeley Earth and WMO datasets when available. When inconsistencies in the Berkeley Earth dataset are found, the related coordinates are corrected by using WMO data. Different scenarios could occur. These are described below and we refer to the flowchart in Figure 1 for further details.

1.
If the weather station from the Berkeley Earth dataset is not available in the WMO list and its spatial resolution is coarser than 30 arcsec, the station is removed.

2.
If the weather station from the Berkeley Earth dataset is not available in the WMO list and its spatial resolution is equal or more detailed than the 30 arcsec, the station is included in the final sample.

3.
If the weather station from the Berkeley Earth dataset is included in the WMO list and the spatial resolution provided by the Berkeley Earth dataset is coarser than 30 arcsec, station coordinates are corrected by using those provided by WMO and the station is included in the final sample.

4.
If the weather station from the Berkeley Earth dataset is included in the WMO list and the spatial resolution provided by the Berkeley Earth dataset is equal or more detailed than 30 arcsec, station coordinates from the two datasets are compared. In case of significant differences (i.e., the station is located in two different grid cells, see Figure 2 as an example), the original sources are consulted to precisely locate the station. Station metadata are, thus, corrected and the station is included in the final sample.    Table S1 in the Supplementary Materials reports the final sample of air temperature stations used in this work, which fulfill the geo-localization procedure (i.e., matching between Berkeley Earth and WMO metadata) [68,69].

Pairing of Temperature and Nightlights Data
Once the position of the stations has been checked and corrected if necessary, we define a regular square buffer around the pixel where the station is located, which ranges from 3 to 7 pixels (i.e., from 3 to 7 km at the equator approximately). We compute the average annual Digital Number value (L for lights, from now on) for the year i and the station s (L is ) using the equation below.
where L ik is the value in the kth pixel for the year i and k tot is the total number of pixels in the buffer around station s. For example, k tot = 9 in a 3 × 3 km buffer around the station and up to k tot = 49 in a 7 × 7 km buffer. In this regard, L is indicates that the average is made in 3×3 km buffer, whereas L (7) is refers to a 7 × 7 km buffer. Buffers smaller than 3 km (i.e., 1 × 1 km and 2 × 2 km) are not considered in order to avoid the spatial noise due to a possible inaccurate geo-localization of air temperature stations. Selected buffers allow one to attenuate this spatial noise (i.e., the local noisy effect is reduced by considering neighboring pixels). The choice of analyzing more than one buffer is due to effects of urban warming that could be detected several kilometers away from the instrument site, even if the major impact is evident in the first kilometer [70].

Gap-Filling Procedure for Incomplete Time Series of Temperature Records
After having matched temperature and nightlights data, we turn to fill possible gaps in the annual time series of temperature records. The Berkeley Earth dataset provides monthly average data T ij , where i is the year (from 1992 to 2013) and j is the month (from 1 to 12). Yearly data are derived by averaging available monthly data. Since some stations show gaps in monthly data, we apply a statistical procedure to estimate the annual mean temperature in the presence of limited missing data and then fill the gaps (see Text S1 in the Supplementary Materials).
In detail, for a given station s, we call Nm i (s) the number of monthly data available in year i, i = 1992 . . . 2013 (0 ≤ Nm i (s) ≤ 12). We also call Ny j (s) the number of available records for month j, j = 1 . . . 12 (0 ≤ Ny j (s) ≤ 22). The distribution of Nm i and Ny j values is reported in Figure 3. By performing a sensitivity analysis (see a detailed description of the method and outcomes in S1 Text in the Supplementary Materials), the acceptable number of missing values to perform the gap-filling procedure is identified. The thresholds Nm* = 9 and Ny* = 18 are selected: if Nm i (s) < Nm* (for any i) or Ny j (s) < Ny* (for any j) the air temperature station s is discarded from our analysis. Otherwise, missing data are reconstructed as follows: suppose the temperature observation is missing in station s for month j* in year i*. The reconstructed value is the average of the temperature values available for the same month in other years using the equation below.
where Ny j* (s) is the number of available records for month j* and T kj* (s) are available temperature data for month j* in year k, k = 1, . . . ,Ny j* (s), respectively. Table S2 in the Supplementary Materials provides an example of application of the performed gap-filling procedure.
Overall, 5530 air temperature stations are selected in this study, meeting both geo-localization and gap-filling requirements (Figure 3c and Table S3 in the Supplementary Materials). The performed data preparation procedure produces a loss of information, from the original 39,000 stations, which show different behaviors across continents. In the case of North America, the loss of information does not significantly influence the analysis, in view of both the high station density and spatial coverage. Conversely, the consistent number of gaps in the temperature datasets in Africa and South America, along with the scarce network coverage, more significantly impacts the consistency of the selected dataset.

Methods
In order to investigate the nexus between temperature variations and urbanization trends, as derived from nightlights, a linear regression analysis estimating trends of temperature and nightlights versus time is first performed. A statistical analysis is then proposed to measure the degree of agreement between trends in temperature and nightlights. Results are analyzed on both global and regional (i.e., continental) scales. For the sake of clarity, examples of application are included throughout the text.

Trend Analysis
where b T identifies the slope of the temperature regression line, a T is the intercept, and t is time. The corresponding linear regression model to fit L values versus time is shown below.
where b L and a L indicate the slope and intercept of the nightlights regression line, respectively. Regression coefficients are estimated via the ordinary least squares method. The slope of the regression line represents the percentage of variation of the temperature (Figure 4a) or the nightlights (Figure 4b) per year.   In order to test the significance of trends, we compute the p-values-p T and p L for temperature and nightlights, respectively-corresponding to the empirically determined slope values on a two-tailed Student's t distribution with n-2 degrees of freedom, where n = 22 is the sample size, e.g., the length of the observation period in years. The null hypothesis of the test is that there is no trend and we adopt a significance level α = 0.1, i.e., a 0.05 significance level on each tail of the distribution. In the following, we allocate positive significant trends in the class c = 1 (p-value ≥ 0.95), negative significant trends in the class c = 4 (p-value ≤ 0.05), while positive and negative non-significant trends are placed in the classes c = 2 (0.5 < p-value < 0.95) and c = 3 (0.05 < p-value < 0.5), respectively. Figure 4 shows the application of the linear regression model to a specific station in a 3 × 3 km spatial buffer. In this case, increasing variations in time for both T and L are found, with p-values of the slope regression lines equal to 0.9846 and 0.9998 for temperature and nightlights, respectively. Since both p T and p L are in class 1 (p value ≥ 0.95), this means that significantly positive temperature and nightlight variations occur in the considered buffer.

Statistical Indicators to Measure the Agreement between Temperature and Nightlights Trends
In order to measure the degree of agreement between the trends of temperature and nightlights, we define an algorithm whose main steps are described below. In the following paragraphs, the results obtained for Asia are used to provide an example of application of the method.
With respect to the considered geographic area (Asia in this case), we compute the percentage of stations with significant (or non-significant) increasing (or decreasing) trends based on related p values. This is performed separately on T and L. We denote these percentages as w Tc and w Lc , with c = 1, . . . , 4 by using: w Lc = n Lc n TOT (6) where n Tc and n Lc indicate the number of p T and p L values in the cth class and n TOT the total sample size in the study area. For Asia, n TOT = 1153 (whereas, for the whole globe, n TOT = 5530). Table 1 shows the percentage of stations in each class of significance c for Asia. For example, 40.8% of stations show a significant and increasing temperature trend and nearly 56% of stations are located in areas of significant increasing luminosity.    We define two variables, denoted as V T and V L , and we assign values to each station and class of significance c: (i) 1 for c = 1, (ii) 0.5 for c = 2 (iii) −0.5 for c = 3, and (iv) −1 for c = 4. Specifically: The product of V T and V L define the final score assigned to the station. We then compute the expected values of the two variables V T and V L based on the probability of occurrence in class c and value assigned to the stations, which is shown in the equations below.
Likewise, variances are computed using the equation below.
In the case of Asia, the expected values and corresponding variances are We then define a concordance index CI (which is formally a covariance), which allows one to assess the degree of agreement between the two considered variables V T and V L .
where the product between V T i and V L i represents the contribution of each station to the concordance index. A decreasing index entails a decreasing agreement between p T and p L , which means increasing temperature T (c = 1,2) and decreasing nightlights trends L (c = 3,4) or vice versa. An increasing index is instead associated with increasing agreement. In the case of Asia, CI = 0.22. If T and L were statistically independent and, therefore, not related, the mean and variance of CI would be computed by using the equation below.
Therefore, Equation (14) provides a reference value to be compared with the CI given by Equation (13) to assess the degree of agreement between T and L. Values of CI can be standardized to make norm-referenced interpretations. We compute a standardized score z using the equation below.
The larger the z value, the more concordant are the variations of temperature and nightlights. Assuming CI is approximately Gaussian-distributed and borrowing some limit values commonly adopted in z-scores interpretation, we have that: if z < −2, T and L variations are strongly discordant. If −2 < z < −1, T and L variations are discordant. If −1 < z < 0, the discordance is weak. If 0 < z < 1, the concordance is weak. If 1 < z < 2, T and L variations are concordant and, if z >2, T and L variations are strongly concordant.
In the case of Asia, we obtain E[CI] = 0.18 σ[CI] = 0.019 and z = 2.1, which has strongly concordant variations in temperature and nightlights.

Results
The main outcomes are shown and discussed in the following paragraphs. Hereinafter, we show the results for only the 3 × 3 km buffer. Analyses on larger spatial buffers lead to similar results compared to the 3 × 3 km buffer (see Figures from S1 to S8 in the Supplementary Materials).
Regression analyses on mean annual T and L are carried out for each station as outlined in the Methods section. Results are reported in Figure 5a-i, where the slopes of T and L regression lines are reported. Sectors 1 and 3 represent positive and negative concordant trends, while sectors 2 and 4 refer to discordant trends, e.g., a rise in temperature in correspondence of a decreasing nightlights trend (sector 4) and vice versa (sector 2). Interesting differences emerge at the continental scale. A positive agreement is detected in more than 50%, 69%, and 43% of stations in Asia, Africa, and South America. In those regions, most stations are located in more and more anthropized areas and experienced an increase of temperature in the period from 1992 to 2013. A unique pattern can be noticed in South America, where 30% of stations are located in areas with increasing nightlights but decreasing temperature trends. In Europe, concordant and discordant patterns are almost balanced. In Oceania, most stations are experiencing warming i.e., more than 70% out of the total, with both increasing (37%) and decreasing (35%) nightlight trends. More than 62% of North American stations show negative luminosity tendency (mostly due to measures to reduce light pollution) in more than 49% of the cases in conjunction with warming trends. This tendency is reflected in the Northern Hemisphere and worldwide due to the large number of stations located in the USA. Figure 6 summarizes the outcomes from our algorithm to assess the concordance between trends of different variables. As the first outcome, the majority of air temperature stations is experiencing warming trends: p T for more than 70% of the selected stations is in class c equal to 1 and 2 (w T1 , w T2 ), with the sole exception of South America (Figure 6f), where negative and positive temperature trends are almost balanced. The global distribution of L is clearly bimodal, with two peaks in class 1 and class 4 (the latter are mostly related to stations in North America).
The computation of the concordance index CI (Equation (13)) reveals an appreciable degree of concordance between T and L trends in Asia and Europe and a weak positive concordance in Africa, South America, and Oceania. The discordance detected in North America moves the CI toward a negative value on a global scale (Figure 6a,h). The standardized z scores, thus, move in the direction of an overall tendency toward positive concordance, except for North America.        (13)). E[CI]: expected mean of CI (Equation (14)). σ [CI]: standard deviation of CI

Discussion and Conclusions
In this work, we introduce a new way of quantifying the nexus between the evolution in time of thermal impacts and urbanization at multiple scales. Our results confirm the overall tendency of urbanization trends to affect temperature changes. A significant increment of temperature in concomitance with increasing luminosity variations are detected worldwide, and regional-scale results are generally in agreement with this overall trend.
We quantitatively assess the degree of agreement across the four different classes of trend significance using a statistical indicator, which computes the percentage of stations included in each class of significance, based on temperature and nightlights trends. The majority of stations are experiencing warming trends since p T values are mainly included in the first two classes (w T1 , w T2 ).
Recent studies at a very local scale show that the size of the source radius (i.e., buffer) matters when assessing urban influences on station temperature trends, and the near-station surroundings have been demonstrated to be more influential for temperature values [52]. Nevertheless, analyses with larger spatial buffers tend to confirm the patterns detected in the 3 × 3 km buffer (Figures S1-S8 in the Supplementary Materials). Cities expand toward neighboring non-urbanized areas, which cause their edges to become brighter [71]. This shows that (i) anthropogenic pressure still grows towards new and peripheric areas, and (ii) most parts of the suburbs continue to experience warming relative to nearby rural sites [70].
As confirmed by recent findings at the regional scale in Italy, with urban dispersion, peri-urban areas become progressively more and more vulnerable to climate change [72]. Besides the background climate, urban landscape morphology likely plays a major role in shaping the spatial variability of urbanization-induced warming such as for some Chinese regions [27]. From an urban extent point of view, results on a number of European cities confirm that UHI in an urban sprawl, stretched and polycentric cities is less intense [36]. Not only urban expansion, but variation in density within cities is another relevant issue to be further addressed in the future [73].
The analyses on the agreement of temperatures and luminosity variations on a continental scale could provide interesting insights into the sociodemographic and economic features of single continents. For instance, Africa and Asia reveal significant increasing temperature trends along with nightlight increments in recent years, which reflect the fast-evolving and uncontrolled urbanization in these areas. In such areas, recent studies highlight that urban land expansion is explained by different factors. For example, China's urban land expansion is driven more by economic rise than population while, in India, population growth drives urban land sprawl [74]. Most urban growth around the world, and especially in middle-income and low-income continents such as Africa and Asia, is taking place on prime agricultural land [75]. By being less urbanized, region-specific projections reveal that these areas are expected to experience relatively faster rates of urbanization over the next 30 years [76,77].
Increasing temperatures along with slightly positive or negative nightlight trends mainly correspond to developed continents. It is worth noting that developed countries have already experienced major urbanization expansion in past decades. In such context, the heating effect due to urban sprawl tends to flatten with the increase of urban areas beyond a certain size because the increase in the density of the urban built environment and activities is not infinite [18,20]. As cities grow, temperatures in pre-existing urban areas in some cases are pretty stable, while temperatures only rise in newly urbanized areas [78]. This is confirmed by recent studies showing that, when temperatures rise in pre-existing urban areas, the air temperature change is smaller than newly urbanized ones [79]. Moreover, based on UN projections [9] and other studies, the urban population in some cities across Europe is expected to decline by 2050 [9,76,77].
The significant decrement of night time luminosity detected in North America could be the result of policy-driven initiatives as light pollution abatement programs promoted in these last years in western countries such as in Canada and in at least 18 U.S. states [55,71,80,81]. The same initiative has been undertaken in many countries in Northern Europe such as the United Kingdom and the Netherlands [59,62]. This makes the interpretation of urbanization dynamics more complex on temperature trends. In locations such as North America and parts of Europe, the mere use of nightlights as a proxy of urbanization may lead to the misleading conclusion that T increases independently of the entity of urbanization, which suggests a minimum UHI effect in regional warming. Nevertheless, this is not confirmed by regional analyses worldwide, where the increment of urbanization seems to be of value in explaining temperature variations [28]. Therefore, in such cases, contingent regionally-based studies are crucial in solving the issue [27,28]. Regionally-based studies could reveal more about the importance of scale-dependency of the nexus between climate patterns and urbanization [20,27]. Studies on cities across North America reveal that the local background climate greatly contributes to UHI [82] and that temperature trends differ between northern and southern latitudes because of biomes [20]. Relevant urbanization patterns exist in Europe that can be captured by nightlights and recent studies confirm that there is a strong variation of urbanization dynamics within single countries and regions, which is in line with our findings [66]. Subnational levels of analysis in Europe reveal that changes in nightlights luminosity are influenced by a number of drivers e.g., the socioeconomic context and national and regional policies [66]. Preliminary work reveal that, when analyses on Central-Northern and Southern countries are performed separately, the Mediterranean area countries reveal a clear tendency toward a significant positive concordance, while other factors, such as light pollution abatement programs, could be responsible for the discordance detected in Northern Europe [83]. On the other hand, cooler Northern European cities appear to be more vulnerable to UHI effects than Southern cities, which appear to be well-adapted [84]. These studies, along with our findings, stress the importance of mutual relations and feedbacks between regional warming and urbanization dynamics. Decreasing temperature trends detected in South America could be related to mesoscale effects (i.e., related to meteorological phenomena from approximately 10 to 1000 km in a horizontal extent such as between a microscale and synoptic scale). Recent findings outline that the intensification of the South Pacific Anticyclone during these last years, which is a consequence of global warming, could contribute to the coastal cooling and warming in the continental Chile and Andes [85,86]. Other previous studies show that, in the past, there has been net cooling over Argentina (about −0.04 • C/decade) in contrast to other land areas. This occurs mainly in the areas where precipitation has increased the most likely due to an increase in the moisture transport from the Amazon to Northern Argentina by the low-level jet [87]. Nightlights trends in South America are growing with a slower rate than in Asia and Africa. South America is more urbanized than these regions since the level of urbanization was around 81%, which matches that of North America (82%) as well as many European countries (68%) [9]. Consequently, the rate of urbanization in South America is quite slow compared to other developing regions, as confirmed by previous studies and projections [76].
Although satellite data provide information that is not directly related to the quantitative rise of temperature records, in this work, we prove how these data can support both global and local analyses of urban and global warming-related issues. Our analysis performed over a wide range of spatial scales presents some uncertainty associated with the reliability of urban temperature records, which provides the ground for future discussion on the effect of urban heating on climate data. Further advances in this direction could benefit from the perspectives offered by new approaches and techniques, and from further studies showcasing the multi-scale impacts of urbanization on the climate. Merging high-resolution data as nightlights, which was made available by new advances in remote-sensing, statistical models, and concepts, could represent the opportunity for an unconventional strategy to analyze such issues.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4433/10/3/117/s1, Text S1: Reconstruction procedure in annual average temperature data. Table S1: Matching between the final sample of Berkeley Earth and official WMO stations used in this work (5530 stations). Stations are provided with related ID codes, coordinates, and precision. In the case of WMO stations, the precision of the coordinates is not reported since the minimum resolution is always equal or higher than 30 arcsec. When no matching with the WMO station is available, WMO coordinates are denoted as "NaN". Table S2: Table. Example of application of a reconstruction procedure of annual average temperature. In the selected station Torino Caselle (Berkeley ID 155990) two gaps in November and December 2013 (i = 22) are found. Table S3: Active and selected air temperature stations in the study period from 1992 to 2013 for "Quality Controlled" data derived from the Berkeley Earth dataset. Number of active (i.e., having at least one year of data) stations from 1992 to 2013 and available stations after the application of thresholds for the reconstruction of mean annual temperature from the mean monthly data and spatial localization. The selected thresholds are Nm i ≥ 9 and Ny j ≥ 18. The percentages refer to the total sample of active and selected stations respectively. Figure S1: Slope of T (b T ) and L (b L ) regression trend lines at a global and continental scale, 4 × 4 km buffer. Sectors 1 and 3 correspond to agreeing trends while sectors 2 and 4 refer to discordant trends. The number of stations included in each sector is in bold as well as the number of stations with L systematically equal to zero (b L = 0) on the horizontal axis. Figure S2: P values density plots at a global and continental scale, 4x4 km buffer. The color scale represents the data density. Figure S3: Slope of T (b T ) and L (b L ) regression trend lines on a global and continental scale, 5 × 5 km buffer. Sectors 1 and 3 correspond to agreeing trends while sectors 2 and 4 refer to discordant trends. The number of stations included in each sector is in bold as well as the number of stations with L systematically equal to zero (b L = 0) on the horizontal axis. Figure  S4: P values density plots at global and continental scales, 5 × 5 km buffer. The color scale represents the data density. Figure S5: Slope of T (b T ) and L (b L ) regression trend lines at a global and continental scale. 6 × 6 km buffer. Sectors 1 and 3 correspond to agreeing trends, while sectors 2 and 4 refer to discordant trends. The number of stations included in each sector is in bold, as well as the number of stations with L systematically equal to zero (b L = 0) on the horizontal axis. Figure S6: P values density plots at a global and continental scale. 6 × 6 km buffer. The color scale represents the data density. Figure S7: Slope of T (b T ) and L (b L ) regression trend lines at a global and continental scale. 7 × 7 km buffer. Sectors 1 and 3 correspond to agreeing trends, while sectors 2 and 4 refer to discordant trends. The number of stations included in each sector is in bold, as well as the number of stations with L systematically equal to zero (b L = 0) on the horizontal axis., Figure S8: P values density plots at a global and continental scale, 7 × 7 km buffer. The color scale represents the data density.