MONITORING URBAN HEAT ISLAND THROUGH GOOGLE EARTH ENGINE : POTENTIALITIES AND DIFFICULTIES IN DIFFERENT CITIES OF THE UNITED STATES

The aim of this work is to exploit the large-scale analysis capabilities of the innovative Google Earth Engine platform in order to investigate the temporal variations of the Urban Heat Island phenomenon as a whole. A intuitive methodology implementing a largescale correlation analysis between the Land Surface Temperature and Land Cover alterations was thus developed. The results obtained for the Phoenix MA are promising and show how the urbanization heavily affects the magnitude of the UHI effects with significant increases in LST. The proposed methodology is therefore able to efficiently monitor the UHI phenomenon.


INTRODUCTION
Google Earth Engine (GEE) is the computing platform recently released by Google "for petabyte-scale scientific analysis and visualization of geospatial datasets".Using a dedicated High Performance Computing (HPC) infrastructure, it enables researchers to easily and quickly access more than thirty years of free and public data archives, including historical imagery and scientific datasets, for global and large scale remote sensing applications.In this way, many of the limitations related to data downloading, storage and processing are effortlessly overcome (Gorelick et al., 2017), (Nascetti et al., 2017).
In particular, the aim of this work is to exploit the large-scale analysis capabilities of such platform in order to develop a new methodology able to investigate the temporal variations of the Urban Heat Island (UHI) phenomenon as a whole.
The term UHI refers to the mesoscale phenomenon associated with higher atmospheric and surface temperatures occurring in urban environments than in the surrounding rural areas due to urbanization (Voogt and Oke, 2003).The effect is most relevant at night when the urban surfaces release the energy stored during the daytime with less efficiency than the nearby rural areas.Therefore, it is easy to understand how the Land Cover (LC) alterations caused by human activities can have an heavy impact on the phenomenon.
For this reason, the work was focused on different USA Metropolitan Areas (MAs) which experienced a significant urban expansion in the last decades.In these cities, soils that were once permeable and wet were transformed in waterproof and dry surfaces, where residential suburbs have replaced forests and/or agricultural hinterlands.Hence, such MAs represent a valuable test site in which to prove the effectiveness of the developed methodology.
The work is thus organized as follows.In Section 2 the methodology is illustrated with great details.In Section 3, the first results obtained for the Phoenix MA (Arizona) are presented and discussed.Indeed, from 1983 to 2010 this city underwent a significant expansion, changing from a mostly agricultural region to a metropolis predominantly characterized by residential suburbs (Fernando et al., 2001), (Doran et al., 2003), (Lee et al., 2003), (Brazel et al., 2007), (Di Sabatino et al., 2009), (Lee et al., 2012), (Wang et al., 2016).Finally, in Section 4 conclusions are drawn and future prospects are outlined.

DATA AND METHODS
The developed methodology implements a large-scale correlation analysis between the Land Surface Temperature (LST) and LC alterations through the joint use of GEE and Climate Engine (CE), a free web application powered by GEE, enabling users to process, visualize, and download various global and regional climate and remote sensing datasets and products in real-time (Huntington et al., 2017).Specifically: • the annual median of the LST was computed through CE (Figure 1) from the Landsat Top of Atmosphere Reflectance Data for every year of the temporal period comprised between the 1992 and the 2011 over the Region Of Interest (ROI)1 corresponding to the considered MA; • the LC data were directly retrieved through GEE from the USGS National Land Cover Database (NLCD, Figure 5) on the same ROI for the 1992 and 2011 years, respectively .
In this way, 20 thermal images (each one relative to a single year) were obtained in which the value stored in every pixel is the median of the LST computed over all the Landsat images available in CE for the considered year2 .Figure 2 reports two examples of such images, respectively for the years 1992 and 2011, where it is possible to notice a general increase of the LST.Then, for every pixel of the ROI, the parameters of a simple linear model (Eq. 1) able to describe the temporal LST variation were robustly estimated.
Thus, the two maps illustrated in Figure 3 were computed: they show the same increasing LST trend (m > 0) observed in Figure 2, characterized by variable rates within the ROI.Nevertheless, there is also a very small number of pixels, statistically not significant, presenting a decreasing or constant LST trend (Figure 4, m ≤ 0).Hence, it is important to investigate the possible relation of the observed LST trends, and their relative growth rates, with the corresponding LC alterations.
As previously mentioned, the LC alterations were retrieved from the NLCD which adopts a LC classification scheme based on several functional classes.
In particular, the NLCD (Fry et al., 2009), (Homer et al., 2015) is a Landsat-based LC database covering four epochs (1992, 2001, 2006 and 2011).It adopts a LC classification system based on several functional classes but, while NLCD2001, NLCD2006, and NLCD2011 products share the same classification scheme (Figure 6(b)), the one adopted by NLCD1992 is slightly different (Figure 6(a)).This can cause a not complete correspondence between the LC classes, as it happens for the classes Commercial Industrial/Transportation of the NLCD1992 and Medium Intensity Residential of the NLCD2011 (Table 1), for which it does not exist an equivalent class in the other classification scheme.For this reason, only the classes present in both the 1992 and 2011 products and more representative of the Phoenix MA territorial features were considered: • Low Intensity Residential: includes areas characterized by a mixture of constructed materials (30%-80%) and vegetation (70%-20%).These areas most commonly include singlefamily housing units.Population density is lower than in high intensity residential areas.
• High Intensity Residential: includes highly developed areas where people reside in high numbers such as apartment complexes and row houses.Vegetation accounts for less than 20% of the cover; the remaining part corresponds to constructed materials.• Shrubland: areas dominated by shrubs; shrub canopy accounts for 25%-100% of the cover.Considering that the city of Phoenix is located in the north-eastern reaches of the Sonoran Desert, this class corresponds to the (scarce) desert vegetation, i.e. to the desert itself.
• Row Crops: rural areas used for the production of crops, such as corn, soybeans, vegetables, tobacco, and cotton.
Therefore, the LST and LC data thus obtained were spatially analysed to cluster the ROI pixels with similar LC alterations, in order to investigate the correlation with their LST growth rates.In this way it is indeed possible to highlight the UHI effect occurred in the most significant urban expansion areas.Specifically, following the adopted linear temperature growth model (Eq.1), two matrices were computed, one for the constants (c) and one for the slopes (m).Every cell ij of the two matrices aggregates all the ROI pixels presenting the specific LC i in the initial year (1992) and the specific LC j in the final year (2011).Moreover, such ij cell contains the mean of the two model parameters computed considering the values of c (constant matrix) or slope (slope matrix) of all the ROI pixels subjected to the LC transformation from i to j.A schematic representation of these matrices is shown in Figure 7.

RESULTS
For the city of Phoenix, the developed methodology was applied to a ROI of about 13340 km 2 , corresponding to the MA itself plus its rural and desert surroundings.In conclusion, the results achieved in this work are pretty consistent, in terms of LST rates, with those obtained in (Wang et al., 2016), where a similar spatio-temporal approach was developed for modelling the UHI in the Phoenix MA.

CONCLUSIONS
In this work, an intuitive methodology was developed to investigate the temporal variations of the UHI effects as a whole, based on the large-scale analysis capabilities of GEE.
Specifically, the promising results regarding the Phoenix MA were presented: they clearly show how the urbanization heavily influences the UHI magnitude with significant increases in LST.The proposed methodology is therefore able to efficiently monitor the UHI phenomenon and in an increasingly precise and fast way.
Nevertheless, it is however necessary to validate such methodology on other MAs, characterized by different weather conditions and not located in desert regions.Moreover, it could be also worth to analyse the LC data in an aggregate way, in order to avoid losing those LC classes not completely equivalent in the two considered classification schemes.

Figure 1 :
Figure 1: Selected ROI for the Phoenix MA.

Figure 2 :
Figure 2: LST 1992 (a) and 2011 (b) obtained through CE for the considered ROI.

Figure 3 :
Figure 3: Map of the constant (a) and slope (b) parameters of the LST linear model Eq. 1 on the investigated ROI.

Figure 4 :
Figure 4: Trend of the least squares for different pixels: positive trend, constant trend, negative trend.

Figure 7 :
Figure 7: Matrix representation of the LC changes: the LC transition from the class Row Crops to the class High Intensity Residential (from i = 51 to j = 24) is indicated with the red arrow.

Figure 9 and
Figure 9 and Figure 10 illustrate separately the spatial distribution of the four considered LC classes, in 1992 and 2011 respectively: it is evident that the (desert) Shrubland class is predominant, because of the desert nature of Phoenix.Conversely, the High Intensity Residential class included a limited number of pixels in 1992 which however increased in 2011, as also those belonging to the Low Intensity Residential class.Such urban expansion was prevalently obtained at expenses of the rural areas which decreased significantly in the investigated temporal period.The obtained results are shown in Figure 8, where the Slope Matrix and the Constant Matrix are reported.On their main diagonals, there are the ROI pixels whose original LC remained unaltered: they form, as expected, the more numerous group.Then there are the pixels which underwent the LC transitions from the classes Row Crops or Shrublands to the classes Low Intensity Residential or High Intensity Residential, changing their nature from purely rural to urban.The abandonment of the rural areas is instead reflected in the LC transformation from the class Row Crops to the class (desert) Shrublands, which involved a consistent number of pixels.Lastly, there are the "urban-rural" transitions, i.e. from Low Intensity Residential or High Intensity Residential to Shrublands or Row Crops, but they are not statistically significant since they involve a very limited number of pixels.Specifically, analysing the Constant matrix (Figure 8(a)), it is evident how in 1992 the higher LST values are located in the first three rows, which correspond to the areas already urbanized or belonging to the desert.The lowest initial LST values is instead in correspondence of the LC class Row Crops, i.e. in the cultivated areas, systematically irrigated and thus colder.For what regards the values contained in the Slope matrix, they highlight the presence of a general LST increasing trend, independent of the considered LC variation (Figure 8(b)).At the (a) LC code 21: Low Intensity Residential (b) LC code 22: High Intensity Residential (c) LC code 51: Shrublands (d) LC code 82: Row Crops

Figure 9 :
Figure 9: 1992 spatial distribution of LC classes within the ROI.

Figure 10 :
Figure 10: 2011 spatial distribution of LC classes within the ROI.