Monitoring the Impact of Land Cover Change on Surface Urban Heat Island through Google Earth Engine : Proposal of a Global Methodology , First Applications and Problems

All over the world, the rapid urbanization process is challenging the sustainable development of our cities. In 2015, the United Nation highlighted in Goal 11 of the SDGs (Sustainable Development Goals) the importance to “Make cities inclusive, safe, resilient and sustainable”. In order to monitor progress regarding SDG 11, there is a need for proper indicators, representing different aspects of city conditions, obviously including the Land Cover (LC) changes and the urban climate with its most distinct feature, the Urban Heat Island (UHI). One of the aspects of UHI is the Surface Urban Heat Island (SUHI), which has been investigated through airborne and satellite remote sensing over many years. The purpose of this work is to show the present potential of Google Earth Engine (GEE) to process the huge and continuously increasing free satellite Earth Observation (EO) Big Data for long-term and wide spatio-temporal monitoring of SUHI and its connection with LC changes. A large-scale spatio-temporal procedure was implemented under GEE, also benefiting from the already established Climate Engine (CE) tool to extract the Land Surface Temperature (LST) from Landsat imagery and the simple indicator Detrended Rate Matrix was introduced to globally represent the net effect of LC changes on SUHI. The implemented procedure was successfully applied to six metropolitan areas in the U.S., and a general increasing of SUHI due to urban growth was clearly highlighted. As a matter of fact, GEE indeed allowed us to process more than 6000 Landsat images acquired over the period 1992–2011, performing a long-term and wide spatio-temporal study on SUHI vs. LC change monitoring. The present feasibility of the proposed procedure and the encouraging obtained results, although preliminary and requiring further investigations (calibration problems related to LST determination from Landsat imagery were evidenced), pave the way for a possible global service on SUHI monitoring, able to supply valuable indications to address an increasingly sustainable urban planning of our cities.


Introduction
Cities all over the world are experiencing very fast growth and changes due to the current rapid urbanization.This is the reason why the United Nations included Goal 11: Make cities inclusive, safe, resilient and sustainable within the 17 Sustainable Development Goals (SDGs) issued in 2015.A key issue after the definition of the SDGs is the monitoring of their progress; with regard to SDG 11, there is a need for well globally-defined, easy computable and comparable indicators, representing different aspects of city conditions, obviously including the Land Cover (LC) changes and the urban climate with its most distinct feature, the Urban Heat Island (UHI) [1,2].
Specifically, the term UHI refers to the mesoscale phenomenon associated with higher atmospheric and surface temperatures occurring in urban environments compared to in the surrounding rural areas due to urbanization [3].It is characterized by a large expanse of non-evaporating impervious materials covering a majority of urban areas with a consequent increase in sensible heat flux at the expense of latent heat flux.Therefore, the UHI and its spatio-temporal evolution are connected to the LC changes, mainly to the transformation from non-urban to urban LC.
Globally, urban footprints only account for approximately 2% of the planetary surface [4].It was well documented that with the rapid urbanization and population growth, the increase in built-up land and the natural land replacement with artificial buildings has altered the surface energy budgets and the hydrological cycle and may affect local and regional climate by changing the physical properties of the Earth's surface [5].Over the past decades, worldwide urbanization and climate change have been two interconnected processes describing the role of human activities in altering and modifying the climate system, which, in turn, may result in uncertain feedback for human welfare.
The UHI effects are exacerbated by the anthropogenic heat generated by traffic, industry and domestic heating, influencing the local climate through the compact mass of buildings that affect the balances of momentum, heat and moisture.The higher temperatures in UHIs increase air conditioning demands [6,7], with a retrofit effect of increasing the outdoor temperature even more, and may modify precipitation patterns.As a result, the magnitude and pattern of UHI effects have been major concerns of many urban climatology studies [8].
The UHI has been a well-known phenomenon since the XIX Century, and it has been studied for many decades; an excellent review of the state-of-the-art (measurements, models, result interpretations) is given in [9], which addresses the four different aspects of UHI, including the Surface Urban Heat Island (SUHI) besides others (canopy urban heat island, boundary urban heat island and subsurface urban heat island).
Traditionally, the UHI is investigated using ground-based in situ measurements of air temperature acquired by automatic weather stations (e.g., see [10]).Air temperature measurements allow one indeed to characterize the fine spatial and temporal variations of the UHI, and they are suitable for studying the impacts of local indicators of urbanization (e.g., sky view factor, floor area ratio) and climate factors (e.g., wind, cloud) contributing to the phenomenon [11][12][13][14].Nonetheless, the urban meteorological networks are not always as complete as could be desired, and frequently, their stations are not uniformly distributed within the territory of the investigated cities.Consequently, large zones may remain without coverage, and it is not possible to analyze the different spatial pattern.
Another possible approach is to investigate the UHI phenomenon through numerical simulations.For instance, ongoing analyses about the UHI at the mesoscale level investigate the effectiveness of superficial properties by using the Weather Research and Forecasting (WRF) mesoscale simulation [15,16].Anyway, additional information and/or local measurements are still required for estimating the parameters of the simulation models [17].
On the other hand, remote sensors on board satellites are able to acquire thermal infrared data from which it is possible to retrieve the urban Land Surface Temperature (LST), allowing thus the study of the SUHI [18].The major advantage of the SUHI is that it can be calculated easily for a large number of cities across large spatial domains.This allows to investigate the role of urbanization in influencing the SUHI over large regions [19][20][21][22][23][24][25][26].A wide literature about the determination of the urban LST is therefore available, based on remote sensing data and its connection with SUHI.
At the same time, though, it is not easy to handle and process the huge amount of Earth Observation (EO) data actually available because of the continuously increasing number of sensors, spatial and temporal resolutions and this will be even more evident in the near future with the development of the micro-satellite constellations.As a matter of fact, even today, much of the archived EO imagery is underutilized despite modern computing and analysis infrastructures, mainly due to processing limitations when faced with a huge amount of data on standard computers and of the high technical expertise needed to use traditional supercomputers or large-scale commodity cloud computing resources [27].In the present era of Geo Big Data, new computing instruments to support remote sensing investigations are therefore necessary, in order to fully exploit and make available the information content of the acquired data.
This work is precisely included in this background: the aim is to show the present potential of Google Earth Engine (GEE) to process the huge and continuously increasing free EO Big Data (i.e., Landsat and Sentinels' imagery) for long-term spatio-temporal monitoring SUHI and its connection with the LC changes.
GEE is indeed the computing platform recently released by Google "for petabyte-scale scientific analysis and visualization of geospatial datasets".Through a dedicated High Performance Computing (HPC) infrastructure, it enables researchers to easily and quickly access more than thirty years of free public data archives, including historical imagery and scientific datasets, for developing global and large-scale remote sensing applications.In this way, many of the limitations related to the data downloading, storage and processing, which usually occur when such a large amount of Geo Big Data is analyzed, are effortlessly overcome [28,29].These features have attracted the attention of researchers from different disciplines, which have integrated GEE into several third-party applications [28].Among these, Climate Engine (CE) [30] can play an important role in the field of large-scale climate monitoring: it is a free web application powered by GEE that enables users to process, visualize, download and share several global and regional climate and remote sensing datasets and products in real time.
Therefore, by continuing the work started in [31], a general methodology for the large-scale monitoring of the SUHI was implemented through the joint use of GEE with CE and tested on six different U.S. cities.The work is organized as follows.In Section 2, an overview of the EO Big Data context is given, paying specific attention to the most promising analysis platforms and underlining the reasons why GEE was selected among them.In Section 3, the proposed methodology is described in detail.A large-scale spatio-temporal procedure was implemented thanks to GEE, also benefiting from the already established CE tool, to extract the LST from Landsat imagery, and a simple indicator was introduced to represent the SUHI evolution and its dependence on LC changes (Figure 1).In Section 4, the results obtained for the investigated six US cities are presented and analyzed: more details are outlined and discussed for the metropolitan area of Atlanta, then the main methodological conclusions derived with regard to Atlanta are applied to the other five metropolitan areas.Finally, in Section 5, conclusions are drawn and future prospects are outlined.

Geo Big Data Analysis Platforms
The latest generations of EO satellites are acquiring increasingly significant volumes of data with a full global coverage, higher spatial resolution and impressive temporal resolution.A new era of EO based on the Big Data paradigm is coming, where innovative tools able to support remote sensing investigations will be increasingly necessary in order to exploit all the potentially available geospatial information contained therein, to effectively contribute to study phenomena at the global scale and to monitor them in (near) real time [32,33].Indeed, it is no longer technically feasible or financially affordable to transfer the data and use traditional local processing methods to address this data scaling challenge.The size of the data and complexities in preparation, handling, storage and analysis are the key aspects that should be taken into account [27].This is the reason why, in the last few years, both satellite EO technology and innovative computing infrastructures have advanced significantly.Such solutions, developed by national space agencies, open source communities and private companies have a great potential to streamline data distribution and processing capabilities, while simultaneously lowering the technical barriers for users to exploit the EO Big Data and to develop innovative services.
In this context, ESA started the EO Exploitation Platforms initiative in 2014, a set of research and development activities that in the first phase (up to 2017) aimed to create an ecosystem of interconnected Thematic Exploitation Platforms (TEPs) on European footing, addressing the most important topics on remote sensing (i.e., coastal, forestry, hydrology, geohazards, polar, urban themes and food security).Briefly, TEPs are a collaborative, virtual work environment providing access to EO data and tools, processors and IT resources required using one coherent interface [34].In 2017, the European Commission in collaboration with ESA and EUMETSAT launched an initiative within the Copernicus program to develop data and information access services that facilitate access to EO Big Data.
Earth Servers 1 and 2 are instead two EO Big Data cubes founded by the Horizon 2020 framework [35,36].Their service interface is rigorously based on the Open Geospatial Consortium "Big Geo Data" standards, web coverage service and web coverage processing service [37].
Considering the high potential of EO Big Data also in commercial applications, Amazon included the full archive of Copernicus Sentinel-2 data in the Amazon Web Service platform through a collaboration with the start-up Sinergise [38,39].They developed the Sentinel-Hub, a service-oriented satellite imagery infrastructure that takes care of all the complexity of handling of the satellite imagery archive and makes it available for end-users via easy-to-integrate web services [40].
Since 2015, the innovative GEE processing platform dramatically transformed and pushed one step ahead the EO satellite data user community.GEE consists of a multi-petabyte analysis-ready data catalog with a high-performance, intrinsically parallel computation service.It is accessed and controlled through an Internet-accessible application programming interface and an associated web-based interactive development environment that enables rapid prototyping and visualization of results [28].Following the user demands, GEE has created a technological solution that removes the burden of data preparation, yields rapid results and maintains an active and engaged global community of contributors.Extensive research and development activity has been performed using GEE, and new applications and services can be deployed to deliver great impact to important environmental, economic and social challenges, including at the local, regional and global scales [41][42][43].Such applications highlight the value of exploiting the EO data, pointing out that the main challenge is in providing the proper connection between data, applications and users.Nowadays, GEE appears as the Geo Big Data analysis platform as more functional from an operational perspective, powerful and at the same time easy to use, and this is the reason why it was selected to carry out the investigations herein presented and discussed.

Data and Methodology
As mentioned before, in this work, a general methodology able to investigate the temporal variations of the SUHI as a whole and its connection with LC changes was implemented through the GEE platform.The main novelty of this approach lies indeed in its easy applicability to the analysis of a huge amount of data, thanks to the Big Data processing capabilities of GEE.
Specifically, a large-scale spatio-temporal analysis was carried out in order to investigate the possible connections between the Land Surface Temperature (LST) trends and Land Use/Land Cover (LULC) changes.
The study was performed through the joint use of GEE and CE, focusing the analysis on six different U.S. Metropolitan Areas (MAs) (Figure 2), characterized by different climate conditions and the availability of a long time series of EO data: Atlanta (Georgia), Boston (Massachusetts), Chicago (Illinois), Houston (Texas), Phoenix (Arizona) and San Francisco (California).In the last few decades, all these cities have experienced a significant urban expansion, which has profoundly modified the main features of their territories.Residential suburbs, commercial and industrial areas have replaced forests and/or agricultural hinterlands surrounding the traditional downtown urban centers: soils that were once permeable and wet have been transformed into waterproof and dry surfaces.Such growth in the peri-urban fringe has resulted in a loss of cultivated areas and natural open spaces: the low values of albedo, vegetative cover and moisture availability in combination with the presence of high levels of anthropogenic heating can thus strengthen the UHI intensity.The selected U.S. MAs therefore represent valuable test sites to assess the effectiveness of the proposed methodology, hereafter described.

Land Surface Temperature
The CE web application was used to compute the annual median of the LST over the Landsat 4/5/7 top of atmosphere reflectance data, for every year of the two decades comprised between 1992 and 2011 in each of the Regions Of Interest (ROIs) corresponding to the six above-mentioned U.S. MAs.Each ROI was selected in order to assure enough surface variability and, at the same time, to consider surface portions remained unaltered and purely rural or natural.
Specifically, more than 6000 Landsat images were processed within CE.In general, their seasonal distribution is well balanced between the cold and the warm period, as reported in Table 1.Furthermore, the use of a robust estimator such as the annual median guarantees a good reliability of the methodology, especially at the beginning of the analysis period, where the number of available images is lower and sometimes slight inhomogeneities are present.and cold (c) period (September, October, November, December, January, February).The orbits of Landsat satellites are solar synchronous, and passages are between 10:00 and 10:30 mean local time at the Equator, in order to provide maximum illumination with minimum water vapor (haze and cloud build-up) [44].Hence, for every MA, a set of 20 thermal images (one per year) was retrieved from CE: the value stored in each pixel is the LST annual median computed over all the Landsat images available in CE for the considered year.
As an example, Figure 3 gives an overall look at the investigated long-term phenomenon, as detected by considering the largest available set of Landsat imagery.In particular, it shows two of the LST annual median maps obtained for the Atlanta MA, respectively for the years 1994 (left panel) and 2010 (right); these two years were selected (as close as possible to the start and the end of the investigated period 1992-2011) since for there is a good balance between the number of images acquired during the cold and the warm seasons (see Table 1), and in this way, the seasonal effects are filtered out.Just a simple visual comparison of these two panels highlights a global increase of the LST.
Thus, with the aim of seeking the simplest suitable model for the temporal variation of the LST, an analysis of the surface temperature was performed at the pixel scale over the 20 years, on the basis of the 20 mentioned computed LST annual median maps.The analysis showed that a linear model can be considered as a reasonable approximation of the phenomenon.Thus, for every pixel of each ROI, the parameters of a simple linear model: describing the LST trend as a function of time were robustly estimated; here, T is the LST, T 0 is the estimated LST referring to t 0 (the analysis starting year, here 1992), t is the current time and r is the rate of the temperature variation per year.In particular, the original LST images retrieved from CE were undersampled with a factor of 10, in order to speed up the estimation of the parameters of the linear model (Equation (1)), computed through the robust regression method [45,46].Therefore, since the LST images have a Ground Sampling Distance (GSD) of 30 m, as the raw Landsat images, the linear model parameters were estimated at a 300-m spatial resolution.Figure 4 reports the maps of the T 0 (left panel) and r (right panel) parameters obtained for the Atlanta MA.Hence, it is worth trying to investigate the possible relation of such observed LST trends, and their relative growth rates, with the corresponding LC transformations.

National Land Cover Database
The USGS National Land Cover Database (NLCD) [47] was used to assess the LC changes that occurred in the six investigated MAs during the twenty years under consideration.The NLCD is a 30-m Landsat-based land cover database spanning four periods (1992, 2001, 2006 and 2011).It was directly retrieved from GEE on the ROIs corresponding to the six MAs for the years 1992 and 2011.
Specifically, the database provides an LC classification for functional classes whose identification codes suffered some modifications during the considered interval 1992-2011: indeed, the NLCD classification scheme used for the NLCD1992 product differs slightly from that adopted in the more recent products (i.e., NLCD2001, NLCD2006, NLCD2011).In order to overcome this problem and to compare LC homogeneous data, the various LC classes were grouped into main three classes related to urbanized, cultivated and forest/shrubland areas as follows: Moreover, after grouping into the chosen three main classes, an undersampling with the same factor of 10 as before was performed on the LC database; in this way, the spatial resolution of the derived LC database for the investigated MAs was equal (300 m GSD) to the one of the derived parameter maps, allowing a one-to-one pixel correspondence.This is a simplification, but it is convenient to highlight the main influences of the long-term relevant LC changes (from rural-cultivated and forest/shrubland to urban areas and vice versa) on the variations experienced by the LST over the considered period.

Correlation between Land Surface Temperature Data and Land Use Data
A spatial analysis was performed to find a possible correlation between the observed LST trends and the corresponding changes in the LC data.In other words, the specific aim was to verify whether the ROI pixels that experienced the same LC changes during the twenty years showed similar LST growth rates or not.Operatively, on the basis of the adopted linear temperature growth model (Equation ( 1)), a (3 × 3) matrix (Rate Matrix) was computed for each MA investigated for the rate (r) parameter.In particular, each cell (i, j) of this matrix (Figure 5): • aggregates all the ROI pixels characterized by the LC class i in the initial year (1992) and the LC class j in the final year (2011); • stores in three different layers the hereafter described statistical parameters (mean, median, numerosity, Permanence/Variation Index), computed considering the values of r of all the ROI pixels that underwent the LC change from class i to class j.
A schematic representation of the obtained Rate Matrix is shown in Figure 5.

Numerosity
This layer quantifies the number of ROI pixels effectively involved in the change from the LC class i to the LC class j.Therefore, the values on the main diagonal indicate the number of stable pixels, which did not change their original LC and which constitutes, as expected, the largest portion.Then, there are the pixels that from purely rural classes (Cultivated or Forest/Shrubland) turned into urban ones (lower triangle) or vice versa, (upper triangle); finally, there are also pixels changing within rural classes (from Cultivated to Forest/Shrubland or vice versa).Recalling that the pixels have a 300-m GSD, it is possible to infer the surface of the ROI that underwent a specific LC change from their numerosity.

Mean, Median
In these layers, the mean and the median of the parameters of the linear model describing the LST trend for the considered LC change from class i to class j are stored.The values obtained are very similar for both statistical parameters, pointing out the substantial absence of outliers.Therefore, for the sake of brevity, only the mean value was considered in the analysis.

Permanence/Variation Index
The overall evaluation of the LC changes, which occurred in the investigated ROIs, has been represented through the herein introduced temporal Permanence/Variation Index, defined as follows: where n ij is the number of the ROI pixels involved in the change from the LC class i to the LC class j.Therefore, I v represents an index of permanence for the elements belonging to the main diagonal (here, usually the highest values are observed), whereas it represents an index of variation for the elements in the lower and upper triangles, defining the frequency of occurrence of the changes experienced by each individual LC.

Study Regions of Interest
A preliminary deeper analysis was developed for the Atlanta MA, following the above described procedure and performing also a comparison between the LST and the air temperature measured on the ground at the time of each image acquisition.Then, considering the results obtained on Atlanta MA, the analysis was replicated for all the other five MAs.

Metropolitan Area of Atlanta
It is widely acknowledged that for the past 30 years, the city of Atlanta has undergone heavy alterations in its LC: urban growth, sprawl, loss of forest and cropland have modified drastically its territorial features [48].Such a tendency is confirmed by a visual analysis of the aggregated LC database computed as described in Section 3.2 and shown in Figure 6, where the Urbanized, Forest/Shrubland and Cultivated LC classes are depicted.Obviously, also the I v reflects such urban growth: over only 20 years, about 40% of Cultivated and Forest/Shrubland pixels became Urbanized.Specifically, the developed methodology was applied to an ROI covering an area of about 11,000 km 2 , corresponding to the city of Atlanta and its surroundings.
The elements of the Rate Matrices (Figure 7a) show the presence of the already mentioned global LST increasing trend, which is independent of the specific LC variation.Nonetheless, it can be observed how the highest increasing rate corresponds precisely to the LC change from the class Forest to Urbanized (0.41 • C/year), immediately followed by the Cultivated to Urbanized transformation (0.39 • C/year).Indeed, the substitution of rural and natural areas with dry, impervious surfaces (roofs, asphalted surfaces, etc.) results in less moisture available to keep the soil cool.On the other hand, even considering the LC stable pixels, it is possible to highlight quite high increasing rates, in the range of 0.29-0.37• C/year, which appear unrealistic even in the presence of the well-known global warming phenomenon.Therefore, considering a small portion (about 13.7 km 2 , Figure 8) of the urban area located around a permanent weather station for air temperature measurement [49], an extensive comparison was performed between the LST and the air temperature measured at the time (with less than one hour approximation) of each image acquisition (around 10:00, local time) over the period 1992-2011.It is known [9] that in the morning, the difference between LST and air temperature in urban areas can be positive for up to 10 • C; anyway, the long-term trends of both temperatures should be equal.Actually, the median temperatures difference for 1992 is around 6 • C, whereas the increasing rates were respectively estimated as equal to 0.22 • C/year for LST (probably due to the fact that small green areas are included in this urban area) and to 0.08 • C/year for air temperature, which appears definitely more realistic (Figure 9).A calibration issue in the LST retrieval procedure from Landsat imagery was therefore supposed to explain such a difference.In addition, the LSTs intra-comparison retrieved from Landsat 5 and Landsat 7 in the available common interval (1999-2011) over the same small portion of the urban area evidenced an anomalous increasing rate difference, up to 0.38 • C/year.Consequently, the results related to the LST increasing rates were judged as unreliable; this is the reason why the Rate Matrix was detrended for the general LST increase, in order to highlight the net effects of the LC changes and the specific role of the urbanization.In detail, the LST trend was separately removed from each row of the Rate Matrices, by subtracting the values pertaining to the main diagonal (i.e., the cell representing the permanence in the same LC class during the investigated period) from those related to the other two cells (which represent the changes among the LC classes).In this way, the so-called Detrended Rate Matrix was computed, in which only the LST increase effectively due to the LC variation is considered, and obviously, the values on the main diagonal are zero.It is worth underlining once more that the lower triangle cells (2,1), (3,1) and (3,2) (according to the description given for the Rate Matrix and represented in Figure 5) of the Detrended Rate Matrix represent respectively the net effects of temperature increase due to the LC changes from Cultivated to Urbanized, from Forest/Shrubland to Urbanized and from Forest/Shrubland to Cultivated.On the contrary, the Detrended Rate Matrix upper triangle cells (1,2), (1,3) and (2,3) represent the net effects of the inverse LC changes.
Hence, the values obtained for the Detrended Rate Matrix are shown in Figure 7b, where all the changes represented in the lower triangle of the matrix are positive (0.03 • C/year for Forest/Shrubland to Cultivated LC change) or strongly positive (0.08 • C/year for Cultivated to Urbanized LC change and 0.12 • C/year for Forest/Shrubland to Urbanized LC change).Moreover, as regards the reverse LC changes from Urbanized to the Cultivated or Forest/Shrubland classes, it is possible to notice a decreasing LST trend, as expected (Figure 7b); anyway, given the small number of pixels involved in those two transformations, the results are less significant.
It is worth underlining that the meaning of the Rate Matrix is doubtful for the discussed possible calibration issues of the procedure implemented in CE to retrieve the LST, whereas the Detrended Rate Matrix is clearly able to highlight the connection between LC changes and the increase of urban LST and SUHI.Therefore, only the Detrended Rate Matrix has been considered significant to represent the results of the long-term monitoring of the impact of LC changes on the SUHI.
Finally, since the number of available images is slightly unbalanced between the warmer and the colder periods (Table 1), the analysis was repeated considering them separately.However, the results (Figures 10 and 11) suggest that the seasonal effect does not have a significant influence on the results in the case of the Atlanta MA.

Metropolitan Area of Boston
In the case of Boston MA, the investigated ROI includes both the cities of Boston and Providence (Figure 12).Boston MA is characterized by a much colder climate than Atlanta, and it was selected in order to understand the behavior of the SUHI in a region surrounded by the sea and to analyze if and how such proximity affects the phenomenon.The Detrended Rate Matrix (Figure 13) reveals that the LST again increases for the LC transitions from Cultivated or Forest to Urbanized.On the contrary, when the LC changes from Urbanized to Forest, the LST decreases.Nevertheless, a comparison with Atlanta MA highlights that in the Boston MA, the LST increase due to urbanization is significantly lower, despite the values of the Permanence/Variation Index I v being quite similar (except for the LC change from Forest to Urbanized, which is lower for Boston).This is likely due to the effect of the Atlantic Ocean.

Metropolitan Area of Chicago
In the case of Chicago MA, the developed methodology was applied to an ROI covering an area of about 11,000 km 2 , corresponding to the city of Chicago and its surroundings (Figure 14).Considering the Detrended Rate Matrix (Figure 15), the highest increasing rate (0.04 • C/year) occurs in the case of the LC change from Cultivated to Urbanized, characterized by a high Permanence/Variation Index (I v = 0.34).The highest I v occurs for the Forest/Urbanized change (I v = 0.51), which, however, involves a lower number of pixels; in this case, the LST increasing rate is lower, amounting to 0.03 • C/year (Figure 15).Finally, the changes from urban LC to rural/natural LCs are not relevant since they involve a very limited number of pixels.

Metropolitan Area of Houston
In the last few decades, Houston has been subjected to a significant urban growth [50] that has altered considerably its territorial features, as shown in the aggregated LC maps (Figure 16).For Houston MA, the results are similar to those obtained for Atlanta MA, but with lower increasing rates (Figure 17): the LST increases in the cells aggregating the ROI pixels in which the rural/natural-urban transition occurs and decreases in the opposite transformation, even if in this particular case, the number of pixels involved is low.Indeed, in this case, the LST increasing rate is equal to 0.02 • C/year for the Cultivated/Urbanized changes, and it is much lower than the value obtained in the Forest/Urbanized changes, which amounts to 0.10 • C/year.As regards the Permanence/Variation Index I v , it differs by about 10% in the two cases, even if the final altered surface is slightly wider for the Cultivated/Urbanized transition; indeed, in 20 years, about 40% of the forest has been urbanized against 30% of the cultivated areas, which is evident also in Figure 16.

Metropolitan Area of Phoenix
The developed methodology was applied to the MA of Phoenix [51][52][53][54][55][56][57].Owing to its particular location (it is located in the northeastern reaches of the Sonoran Desert), the Phoenix MA climatic conditions are significantly different from those observed for the other investigated areas.Indeed, the Forest/Shrubland LC class, here mainly related to shrublands, covers wide portions of the ROI (Figure 18).Nonetheless, the developed methodology is once again able to detect effectively the LST trends (Figure 19) connected to LC changes (Figure 19).Indeed, the Cultivated/Urbanized change results are characterized by the highest increasing LST rate (0.11 • C/year).On the contrary, a negative LST increasing rate (−0.06 • C/year) can be observed in the LC change from Shrubland (desert) to Cultivated, while a positive increasing rate (0.09 • C/year) can be noticed in the reverse change.The first situation corresponds to the ROI pixels that become cultivated and thus irrigated, whereas the latter is related to the abandonment of the rural areas, which thus were no longer irrigated.

Metropolitan Area of San Francisco
In the case of San Francisco MA, the investigated ROI includes both the San Francisco Bay area and the city of San Jose (Figure 20).The results show that the number of the ROI pixels involved in the LC change from the Forest or Cultivated classes to the Urbanized class is lower than that observed in the previously investigated MAs (Figure 20).In other words, the urban growth was not so strong during the period considered.The elements of the Detrended Rate Matrix (Figure 21) show the familiar increasing of the LST observed in the transition from rural/natural into urban LU.The highest and significant increasing rate occurs for the Cultivated to Urbanized LC change (0.12 • C/year), anyway involving a small portion of the initially cultivated area (I v = 0.18); I v for the Forest to Urbanized LC transition is even lower, amounting to 0.06, and the increasing rate is also lower (0.05 • C/year).Furthermore, there are not significant reverse LC changes.

Conclusions and Prospects
The aim of this work was to show that it is presently possible to monitor SUHI and its connection with LC changes on a long-term and wide spatio-temporal basis by leveraging the large-scale analysis capabilities of GEE.A general procedure was defined and implemented, thanks to the joint use of GEE and CE, and the indicator Detrended Rate Matrix was introduced to globally represent the net effect of LC changes on SUHI.
The procedure was applied to six diverse U.S. MAs, characterized by different climatic conditions and significant urban expansions over the 1992-2011 period.
Three main simplifications were introduced in this first application of the proposed procedure to test its feasibility: • the spatial analysis was performed at a spatial resolution 10-times lower than the original available information, both regarding LST and LC (a GSD equal to 300 m instead of the original 30-m GSD); • only three main LC classes were considered, suitably grouping those related to Urbanized, Cultivated and Forest/Shrubland areas, which were differently defined in the 1992 and 2011 versions of the USGS NLCD; • a linear model for the LST temporal evolution was adopted.
It is however important to underline that these simplifications are not intrinsic to the procedure and can be released in the future.
A starting analysis to tune the procedure was performed on the Atlanta MA.A general LST increase was detected, also for LC stable pixels, which kept their original LC over the whole period 1992-2011; anyway, the quite high estimated increasing rates, in the range 0.29-0.37• C/year, appeared unrealistic.Therefore, considering a small portion of the urban area located around a permanent station for air temperature measurement, an extensive comparison was performed between the LST and the air temperature measured at the time of each image acquisition over the period 1992-2011, and a significant inconsistency (0.22 • C/year for LST vs. 0.08 • C/year for air temperature) was highlighted.A calibration issue in the LST retrieval procedure from Landsat imagery was supposed to explain such a difference; this was the main reason to introduce the Detrended Rate Matrix, in order to eliminate the effect of the biased LST increases and to highlight the net effect of the LC changes on the SUHI.
The application of the Detrended Rate Matrix evidenced a strong correlation between the highest increasing LST rates and the transformations from rural to urbanized LC found, common to all the studied ROIs.Hence, the results clearly show how urbanization heavily influences the magnitude of the SUHI effects with significant increases in the LST, up to 0.12 • C/year, and demonstrate the effectiveness of the proposed methodology, whose real strength lies in its simplicity, ease of use and rapidity, thanks to the GEE computing capabilities.Indeed, it was possible to analyze a huge amount of data (see Table 1) in a short time (few minutes), not even imaginable until recently.Furthermore, the analysis can be easily extended to other urban areas, provided that a sufficient number of Landsat images is available.
The present feasibility of the proposed procedure and the encouraging obtained results, although preliminary and requiring further investigations, pave the way for a possible global service on SUHI monitoring, able to provide valuable indications to address an increasingly sustainable urban planning of our cities.
In this respect, some possible prospects for the future can be outlined, in order to enhance the proposed procedure: • to understand the better and feasible way(s) to release the above-mentioned simplifications; • to deeply investigate the possible calibration issue of the LST retrieval procedure, widening the comparison with ground air temperature measurements in different areas; • to extend the analysis using also the Sentinel-2 and Sentinel-3 imagery and the CORINE land cover database, already available in the GEE archive.
In this way, one step ahead towards a deeper interpretation of the SUHI phenomena related to specific areas will be possible, but this was not within the objectives of the present investigation.

Figure 1 .
Figure 1.The workflow of the methodology.

Figure 3 .
Figure 3. Maps of the LST annual median for 1994 (left) and 2010 (right), computed through CE over all the Landsat images available for the considered year within the ROI of Atlanta MA.The pixel size is 30 m, as the raw Landsat images from which they derive.

Figure 4 .
Figure 4. Constant (left panel) and rate (right panel) parameters of the LST linear model of Equation (1) for the Atlanta MA (pixel size of 300 m).

•
the Urbanized class includes the classes Low Intensity Residential, High Intensity Residential and Commercial/Industrial/Transportation of the NLCD1992 and the classes Developed Open Space, Developed Low Intensity, Developed Medium Intensity and Developed High Intensity of the NLCD2011; • the Cultivated class includes the classes Pasture/Hay, Row Crops, Small Grains and Fallow of the NLCD1992 and the classes Pasture/Hay and Cultivated Crops of the NLCD2011; • the Forest/Shrubland classes include the classes Deciduous Forest, Evergreen Forest and Mixed Forest for both the considered years and the classes Shrubland of the NLCD1992 and Dwarf Scrub, Shrub/Scrub of the NLCD2011.

Figure 5 .
Figure 5. Representation of the LC changes within the Rate Matrix: the row index (i) represents the initial LC class and the column index (j) the final LC class.As an example, the highlighted cell represents the transition from the Forest class (i = 3) to the Urbanized class (j = 1), with an estimated LST increasing rate of 0.4 • C/year.

Figure 6 .Figure 7 .
Figure 6.Atlanta MA: LC situation in 1992 (a) and in 2011 (b).In red, the Urbanized class, in light green, the Cultivated class, in dark green, the Forest and Shrubland classes.

Figure 8 .
Figure 8.The ROI around the Georgia State University weatherSTEM station, Peoplestown (GA), considered for the LST vs. air temperature comparison (the star denotes the location of the permanent weather station).

Figure 9 .
Figure 9.Comparison of the LST retrieved from the Landsat data with the air temperature measured at the Georgia State University weather station.

Figure 10 .
Figure 10.Atlanta MA warm period: Rate and Detrended Rate Matrices.

Figure 12 .
Figure 12.Boston MA: LC situation in 1992 (a) and in 2011 (b).In red, the Urbanized class, in light green, the Cultivated class, in dark green, the Forest and Shrubland classes.

Figure 14 .
Figure 14.Chicago MA: LC situation in 1992 (a) and in 2011 (b).In red, the Urbanized class, in light green, the Cultivated class, in dark green, the Forest and Shrubland classes.

Figure 16 .
Figure 16.Houston MA: LC situation in 1992 (a) and in 2011 (b).In red, the Urbanized class, in light green, the Cultivated class, in dark green, the Forest class.

Figure 18 .
Figure 18.Phoenix MA: LC situation in 1992 (a) and in 2011 (b).In red, the Urbanized class, in light green, the Cultivated class, in dark green, the Forest and Shrubland classes.

Figure 20 .
Figure 20.San Francisco MA: LC situation in 1992 (a) and in 2011 (b).In red, the Urbanized class, in light green, the Cultivated class, in dark green, the Forest and Shrubland classes.

Table 1 .
Number of Landsat images available in the Climate Engine (CE) platform for each of the investigated U.S. MA, respectively in the warm (w) period (March, April, May, June, July, August)