## Data and Methodology
### 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. Figure 2: The six investigated U.S. Metropolitan Areas (MAs). 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:  
\[T=T_{0}+r\ (t-t_{0}) \tag{1}\]  
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.