Theoretical and Applied Climatology https://doi.org/10.1007/s00704-019-02973-y ORIGINAL PAPER Air temperature in Barcelona metropolitan region from MODIS satellite and GIS data C. Serra 1 & X. Lana 1 & M. D. Martínez 1 & J. Roca 2 & B. Arellano 2 & R. Biere 2 & M. Moix 2 & A. Burgueño 3 Received: 28 September 2018 / Accepted: 7 August 2019 # Springer-Verlag GmbH Austria, part of Springer Nature 2019 Abstract The metropolitan region of Barcelona (BMR) is one of the most densely populated areas in the Mediterranean countries. The estimation of air temperature at a short scale from satellite measurements would contribute to a better understanding of the varied and complex spatial distribution of temperatures in BMR. This estimation would be a first step to study several patterns of the thermometric regime affecting population life quality and health. Taking advantage of MODIS data, air temperature measurements at 48 thermometric stations along the year 2015, together with their geographic and topographic data, multiple regression analyses have permitted to obtain fine spatial distributions (pixels of 1 km2) of minimum, mean and maximum daily air temperatures. Previous to the multiple regression, Pearson coefficients and principal component analysis offer a first overview of the relevance of the variables on the empiric temperatures. The most relevant variables on the multiple regression process at annual and seasonal scale are land surface temperatures, latitude, longitude and calendar day. At a monthly scale, altitude (maximum temperature) and continentality (cold months for minimum and hot months for maximum temperatures) are also relevant. The best fits between empiric temperatures and those derived from the multiple regression processes have square regression coefficients within the range (0.92–0.96) for the annual case, (0.70–0.92) at seasonal scale and (0.52–0.87) at monthly scale. The root mean square error varies from 1.5 to 2.0 °C (annual case), from 1.3 to 2.0 °C (seasonal scale) and from 1.2 to 2.1 °C (monthly scale). In agreement with these regression coefficients and mean square errors, the obtained spatial distribution of temperatures is of notable quality. As an outstanding application, the detection of several urban heat islands on different conurbations within BMR along the Mediterranean coast becomes possible. 1 Introduction The estimation of the land surface temperature (LST) from the thermal infrared radiation (TIR) emitted by the Earth became possible a few years after the arrival of remote sensing from the space. Rao (1972) was the first to apply TIR to estimate temperature patterns for the cities along the USA mid-Atlantic coast from data collected by the Improved TIROS Operational Satellite (Gallo et al. 1995). TIR measurement from the space * X. Lana francisco.javier.lana@upc.edu 1 Department of Physics, Technical University of Catalonia, Barcelona, Spain 2 Department of Architectural Technology, Technical University of Catalonia, Barcelona, Spain 3 Department of Applied Physics–Meteorology, University of Barcelona, Barcelona, Spain permits to know the LST assigned to a pixel (Dash et al. 2001, 2002), without considering the multiple parts contained in this pixel, as vegetation, sunlit and shadowed soils, irregular urban surfaces or soil moisture, with different albedo and emissivity values (Tomlinson et al. 2011; Benali et al. 2012). In fact, TIR is derived from the top of atmosphere radiances, from which LST is obtained after applying corrections due to atmospheric attenuation, angular effects and emissivity values at the heterogeneous surface. Water vapour and aerosols are the main agents causing variable attenuation in the TIR signal. This reduces the LST availability to only under cloud-free conditions, to avoid a systematic bias toward colder-than-true values (Williamson et al. 2013). In this sense, the integration of synergistic information from satellite optical-IR and passive microwave remote sensing has been proved recently to permit consistent and reasonable temperature estimations with cloudy skies (Jang et al. 2014). The right estimation of the temperature of the air at ≈ 2-m height above ground (Ta) from LST is possible but complex. The vertical lapse rate to be applied is function of the surface C. Serra et al. energy balance, which varies in function of the nature of the surface and of the instant of the day, as also of advection, adiabatic processes, turbulence and latent heat fluxes, all of them affected by cloud cover, water vapour content and vegetation (Benali et al. 2012). During the night, the estimation of Ta becomes simpler because the earth surface behaves almost as homogeneous surface (Didari et al. 2017). This Ta estimated from satellite measurements would solve the weather stations scarcity in wider regions, where the geospatial interpolation methods, as kriging or splines, cannot provide accurate estimations, as happens, for instance, in mountainous terrain (Lin et al. 2016) or undeveloped countries. In this way, Ta estimation becomes of crucial importance to solve spatial gaps for a wide range of applications, in such a way that it is accepted that TIR produces better Ta estimations than those obtained by interpolating ground-station temperatures (Mendelsohn et al. 2007). In the first years of remote sensing, this strategy based on TIR permitted to obtain extended and automatic LST for large regions, but with a limited spatial resolution, as in the case of Meteosat satellite (e.g. Cresswell et al. 1999). Since the year 2000, the MODerate resolution Imaging Spectroradiometer (MODIS) sensor in Terra and Aqua polar satellites (http:// modis.gsfc.nasa.gov/) have reduced the spatial resolution for LST to 1 km per pixel. In this way, two images per satellite per day are generated (Terra satellite passes daily over the equator close to 10:30 UTC and 22:30 UTC; Aqua at 13:30 UTC and 01:30 UTC), using both the 10.78–11.28 μm and 11.77–12. 27 μm spectral bands, together with split-window algorithms (Wan et al. 2002). These MODIS products have been submitted to consistent validation (Coll et al. 2005; Wang et al. 2008). In addition to climatological applications, other biological and physical processes on the land and the ocean may be derived using 36 electromagnetic spectral bands from visible to TIR available from MODIS (Zhang et al. 2003; Wan et al. 2004; Wang et al. 2009) or also mapping the global distribution of urban land (Schneider et al. 2009). Given the practical impossibility of Ta direct determination from MODIS Terra LST, different estimation methods have been applied (Zaksek and Schroedter-Homscheidt 2009). The simplest one consists in supposing a linear relation between Ta and LST by distinguishing different land cover types, as Shen and Leptoukh (2011) have applied to Central and Eastern Eurasia, or without this distinction for smaller regions (Fu et al. 2011; Sohrabinia et al. 2015). The consideration of the called temperature-vegetation index (TVX), proposed by Nemani and Running (1989), by applying the normalised difference vegetation index (NDVI) and ignoring its seasonal, ecosystem type and soil moisture variability, has permitted the inclusion of the vegetation cover as a relevant factor (Prihodko and Goward 1997; Vancutsem et al. 2010; Cristóbal et al. 2008; Nieto et al. 2011; Wenbin et al. 2013; Shah et al. 2013; Bustos and Meza 2015). Modifications have been also added, as the differential TVX method (Sun et al. 2014). Certainly, the vegetation cover is determinant by its transpiration cooling and latent heat fluxes, as also through their low albedo and roughness which aides efficient sensible heat dissipation (Benali et al. 2012). Alternatively, the multiple linear regression applies different variables in addition to LST and NVDI to estimate Ta, as latitude, distance from coast, altitude and solar radiation (Cristóbal et al. 2008) or albedo and solar radiation (Xu et al. 2014). Nevertheless, after having considered different predictors, Lin et al. (2012) have proved that just the altitude and LST permit to obtain Ta for East Africa. Zhang et al. (2011) also apply the solar declination variable along the year and LST to derive Ta in China. Kloog et al. (2017) derive daily Ta estimations from LST, NVDI, elevation and the grid cell percentage of urbanicity for France. But, in general, for extended regions, spatialtemporal variables as Julian day of the year, latitude, longitude, height above sea level, slope, curvature and distance to the coast use to be considered (Recondo et al. 2013; Peón et al. 2014; Good 2015; Thanh et al. 2016; Yang et al. 2017). A relation with the variables applied by different authors can be found in Janatian et al. (2017). Also, the use of spatiotemporal regression-kriging and incorporation of time-series of remote sensing images have been proved to permit significantly more accurate maps of temperature than if plain spatial techniques were used (Hengl et al. 2012). At the planetary scale, the new dataset of spatially interpolated monthly climate data for global land areas at a very high spatial resolution (approximately 1 km2) has considered LST observations to cover areas with a low station density (Fick and Hijmans 2017). This effort improves the first LST map at the planetary scale from MODIS measures (Kilibarda et al. 2014). Besides different valuable applications of the estimation of Ta from LST, this study should be the base of future detailed analyses of the urban heat island (UHI) of Barcelona considering its entire Metropolitan Region (BMR). Unless otherwise indicated, UHI intensity is derived as the difference in spatially averaged surface temperatures between urban and non-urbanised surroundings, as a measure of the excess of warmth of the urban atmosphere (Voogt and Oke 2003). UHI phenomenon has been usually analysed from air temperature measurements of a short number of gauges across the city, sometimes with emplacements submitted to criticism, as gardens or roofs, and others outside the city influence (Stewart 2011). This humble departing state (Landsberg 1981), due to the scarcity of points with measurements, has led to a first spatial and temporal characterisation of the phenomenon (Arnfield 2003). Later on, thermometers installed in automobile have permitted to extend the analysis for selected transects across the city (Caselles et al. 1991; Moreno-Garcia 1994) to derive a thinner description of the anomalous urban thermal behaviour in relation to the rural proximity of the city. Urban Air temperature in Barcelona metropolitan region from MODIS satellite and GIS data networks of stations together with different rural temperature observatories around the city have been also undertaken permitting detailed descriptions (Giannaros and Melas 2012; Yang et al. 2013). The multiplicity of urban internal configurations, due to the diversity of geometry, morphology and size of the cities, together with local or regional air dynamics, as the mesoscale sea breeze in the case of littoral cities, makes the urban climate a difficult objective to be rightly modelled. This shortcoming is notably solved by considering Ta series derived from satellite remote sensing. Voogt and Oke (2003) have reviewed the research done with thermal remote sensing before MODIS application. They conclude that the complexity of the urban surface should be analysed through couple canopy radiative transfer models with both sensor view models and surface energy balance models to simulate air temperature in and above the urban canopy layer. With this purpose, Miao et al. (2009) have applied MODIS observations with dynamical models to simulate urban weather features for comparison with observations in Beijing. Nevertheless, most of the studies analyse the spatial UHI taking advantage of the high spatial resolution provided by MODIS but without considering air dynamics for single cities (Cheval et al. 2009; Cheval and Dumitrescu 2009; Fabrizi et al. 2011; Tomlinson et al. 2012; Ma et al. 2016), as also for sets of cities (Jin et al. 2005; Hung et al. 2006; Yasuoka 2006; Pongrácz et al. 2006, 2010; Imhoff et al. 2010), for selected episodes, months or a few years. 2 Database 2.1 Study area The metropolitan region of Barcelona, BMR, with an extension of 3242.2 km2 and a population density of 1566.2 inhabitants/km 2 , according to IdesCat-2017 (Institut Català d’Estadística), is a crowded area close to the Mediterranean Sea. Particularly, Barcelona city, with a population of 1.6 million inhabitants, covers an area close to 100 km2 with a population density close to 16,000 inhabitants/km2. The orography of the analysed region is characterised by the Littoral and Pre-Littoral chains, with moderate altitudes up to 1700 m a.s.l., both parallel to the Mediterranean coast. Between both chains are placed the Vallès valley and Penedès Basin. The most extended urban area (Barcelona city) is constrained among the Mediterranean shoreline, the Littoral chain and Llobregat and Besós rivers. The main orographic features of the region and the distribution of altitude in meters are shown in Fig. 1a,b. Details of the spatial distribution of CORINE land cover classes (http://land.copernicus.eu/pan-european/corineland-cover/clc-2012), at level 3 for the year 2012, are shown in Fig. 1c and Table 1. 2.2 Meteorological station data Observed daily minimum, Tmin; mean, Tmean and maximum, Tmax temperatures are the dependent variables of this study. These come from 48 meteorological stations, 37 of them belonging to the Servei Meteoròlogic de Catalunya (www.meteocat.cat) and 11 to the Agencia Estatal de Meteorología, (www.aemet.es) for the year 2015. Data are obtained in both cases from automatic weather stations, and their quality is guaranteed by periodic instrumental controls of the two governmental institutions. Additionally, the 48 thermometric records are free of perturbations, such as sharp changes or artificial trends, in agreement with the results of the Buishand (1982) and Pettitt (1979) tests, as proposed by Wijngaard et al. (2003). Given that these possible perturbations cannot be detected analysing only 1 year, the two mentioned tests have been applied to longer records including the year 2015, and the results have been also compared with previous analysis of the thermometric regime in a wider area of Catalonia (Martínez et al. 2010). Figure 1 d shows the spatial distribution of the stations, where five of them are outside but very close to the studied region. The stations are well spread over BMR except in the north, where they are scarce. Table 2 gives the main geographical and topographic variables of the thermometric station emplacements. Figure 2 a and b show the histograms of altitudes for the 1 km2 pixels covering all the BMA and for the set of meteorological stations respectively. Altitudes of the available thermometric stations are mostly emplaced below 600 m a.s.l. 2.3 Satellite data The daily MOD11A1 LST measured by MODIS Terra, including daytime 10:30 UTC surface temperature, LSTd, and night time 22:30 UTC, LSTn, with 1 km2 spatial resolution, has been used in this study. The normalised difference vegetation index, NDVI, is obtained from the 16-day MOD13Q1 product (with a resolution of 250 m): NDVI ¼ NIR−RED NIR þ RED ð1Þ where NIR is the near-infrared reflectance band-2 (841– 876 nm) and RED, the reflectance of the red band-1 (620– 670 nm). The normalised difference build-up index, NDBI, is calculated as: NDBI ¼ MIR−NIR MIR þ NIR ð2Þ where MIR is the surface reflectance band-6 (1628– 1652 nm) from MOD09A1 product of 8-day average (with C. Serra et al. b a SIN BA T TO R AL AL AI N 41.8 EY LL VA CH Barcelona AI N AN NE A R ER IT ED M 41.6 Altitude -9999 (m) -0 A SE R VE RI LI OR ÈS LL VA T GA RE OB LL S DÈ NE PE E TT IVER ÒS R BES PR I -L CH 0 - 100 100 - 200 200 - 300 300 - 400 400 - 500 500 - 600 600 - 700 700 - 800 800 - 900 900 - 1000 1000 - 1300 1300 - 1700 41.4 41.2 1.6 c 1.8 2.0 2.2 2.4 2.6 d 45 11 27 26 10 25 30 12 2 34 35 18 21 15 14 28 39 375 38 6 46 3 47 31 23 8 48 42 17 29 7 32 41 9 44 43 33 24 40 13 19 22 20 1 4 16 36 Fig. 1 Spatial distribution of main orographic features (a), altitude above sea level (b), CORINE land cover classes (c) and the thermometric stations (d) a resolution of 500 m). Figure 2 c and d show the histograms of NDVI and NDBI for emplacements of temperature stations. A half of these coefficients are within the 0.3–0.5 interval, corresponding to emplacements with low vegetation cover. For the whole set of 1 km2 pixels, the modal value of NDVI is shifted toward 0.6. In consequence, rural areas are slightly predominant in comparison with urban domains. With respect to NDBI, it is worth mentioning that most of the station emplacements have coefficients ranging from − 0.2 to 0.1. For the whole domain, the mode of NDBI is − 0.1, suggesting a slight predominance of nonurban areas. 2.4 GIS data and calendar day Besides Satellite variables LSTd, LSTn, NDVI and NDBI, other six geographical and topographic variables are considered. These are latitude (lat), longitude (lon), distance to coast or continentality (con), altitude (alt), orientation (ori) and slope (slp) of the terrain for every meteorological station and pixel. The first three are derived from ArcGIS software (Geographic Information Systems, GIS). Altitude, orientation and slope are obtained from the Ground Digital Model (MDT—Institut Cartogràfic I Geològic de Catalunya, ICGC) with a 15 × 15 m2 resolution. Table 3 summarises the minimum, mean and maximum of LSTd and LSTn, NDVI and NDBI, together with geographic and topographic variables. Furthermore, the calendar day, cd, has been transformed into a new calendar day, cd*, according to: cd* ¼ cos 2πðcd−cdmax Þ 365 ð3Þ to obtain the linearity respect to the air temperature (Janatian et al. 2017). cdmax is the calendar day for which the mean temperature along the year is the highest. Figure 3 a and b show the relationship between Tmean and cd or cd* respectively. cdmax for the year 2015 is equal to 200 (July 19th). Figure 3 b shows more signs of linearity between the air temperature and the transformed calendar day, cd*. Figure 4 a shows the dependence of Tmax on the orientation. This dependence is unclear and a linear relationship should be discarded. Trying to solve this lack of linearity, the orientation is given as sine and cosine compounds. Figure 4 b shows the case for Tmax against sine compound, where a small linear increasing tendency is observed. Conversely, the cosine compound does not show signs of linear tendency. 3 Methodology The estimation of surface air temperatures is based on the relationships between variables obtained from satellite (LSTd, LSTn, NDVI and NDBI), geographic and topographic data (latitude, longitude, altitude above sea level, orientation, Air temperature in Barcelona metropolitan region from MODIS satellite and GIS data Table 1 Types of land cover (percentage and areas) on BMR Land cover 2 Area (km ) (%) Coniferous forest Broad-leaved forest 866.91 493.41 26.7 15.2 Discontinuous urban fabric 333.53 10.3 Vineyards Sclerophyllous vegetation 293.50 291.79 9.0 9.0 Non-irrigated arable land 173.43 5.3 Industrial or commercial units Continuous urban fabric 166.17 132.75 5.1 4.1 Permanently irrigated land 95.16 2.9 Occupied by agriculture Transitional woodland-shrub 68.76 68.26 2.1 2.1 Pastures 43.71 1.3 Complex cultivation patterns Green urban areas 36.72 26.34 1.1 0.8 Mixed forest 22.44 0.7 Fruit trees and berry plantations 19.66 0.6 Sport and leisure facilities 19.48 0.6 Road and rail networks Mineral extraction sites 18.54 14.05 0.6 0.4 9.34 9.33 8.35 7.44 5.66 5.12 4.22 3.45 3.37 3244.72 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 100 Construction sites Airports Moors and heathland Port areas Olive groves Beaches, dunes, sands Salt marshes Sparsely vegetated areas Watercourses Total slope and distance to coast), the modified calendar day (cd*) and empiric data (Tmin, Tmean and Tmax) recorded at the thermometric stations. The first step consists of computing the Pearson correlation coefficient for all possible pairs of data, including empiric data. In this way, possible relationships between assumed independent variables can be detected. Additionally, the dependence of empiric data on the set of independent variables can be determined. The rotated principal component analysis (RPCA) (Jolliffe 1986; Richman 1986; Preisendorfer 1988) is the second step. In this way, more detailed characteristics of relationships between independent variables can be established; particularly, the ratio of data variance explained by every rotated principal component, RPC, and the contribution (factor loading) of every independent variable in the RPCs. These strategies, Pearson correlation and PCA, have been also applied by Thanh et al. (2016). The third step consists of a multiple regression process with software from Statistic package for Social Sciences, IBMSPSS, with an assumed linear relationship between empiric data and independent variables. The multiple regression goodness of fit is quantified by the square regression coefficient, R2, including all the significant independent variables according to P values and α = 0.05 (Harrell 2001) and by residuals between empirical, emp, and estimated, est, temperatures computed from the root mean square error, RMSE, and the mean average error, MAE. These errors are computed as: n o1=2 RMSE ¼ N −1 ∑Ni¼1 ðempi −est i Þ2  MAE ¼ N −1 ∑Ni¼1 jempi −est i j ð4aÞ ð4bÞ The multiple regression process is repeated adding step by step a new independent variable. In this way, the relevance of every variable on the multiple regression process is contrasted by observing the changes in R2 and RMSE. In addition to the coefficients of the multiple regression equations, the beta weights (standardised coefficients) are also computed to determine the relevance of every independent variable on the multiple regression equation. As a summary, in comparison with other similar researches based on satellite and thermometric stations data, it should be mentioned that in this paper, first, cross-correlation and principal component analysis (PCA) permits the detection of possible redundant variables in the multiple regression process. Second, new geographic variables (orientation and slope) are tested. Third, both thermometric satellite data (LSTd and LSTn) are used in the multilinear regression, whatever minimum, Tmin; mean, Tmean and maximum, Tmax, daily temperatures are deduced from the multiple linear regression. It should be also remembered that only variables with absolute values of beta weight exceeding 0.01 have been finally considered for the multiple regression equations. 4 Results and discussion 4.1 Correlation coefficients and PCA The Pearson correlation and the RPCA, based on the principal component analysis, PCA, have permitted to detect the degree of dependence between empiric data (Tmin, Tmean and Tmax) and the rest of parameters (geographic and topographic variables and data from a satellite). Table 4 shows the Pearson correlation coefficients among all variables. The high correlations, ranging from 0.86 to 0.97, between Tmin, Tmean and Tmax with daytime and nighttime LST and also with cd* are outstanding. The correlation of the empiric temperatures with the other parameters is always inside ± 0.23. C. Serra et al. Table 2 GIS characteristics of the stations: longitude (Lon), latitude (Lat), altitude above sea level (Alt), distance to shoreline (Con), orientation of the slope (Ori) and topographic slope (Slp). Ori equal to − 1.0 indicates flat terrains (slope equal to 0.0) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 Station Lon (°) Lat (°) Alt (m) Con (m) Ori (°) Slp (%) Barcelona—Airport Pontons Vilafranca del Penedès Sitges—Vallcarca Barcelona—CMT Barcelona—Drassanes Sabadell—Aeródromo Vilassar de Dalt Arenys de Mar Santa Susana Fontmartina els Hostalets de Pierola Vacarisses Vallirana Barcelona Observatori Fabra el Vendrell Font-rubí Sant Martí Sarroca PN del Garraf—el Rascler Viladecans el Montmell Sant Pere de Ribes—Garraf Cabrils Dosrius—Montnegre Corredor Rellinars Sant Llorenç Savall Tagamanent—PN del Montseny la Granada Vilanova del Vallès Montserrat—Sant Dimes la Bisbal del Penedès Canaletes Malgrat de Mar Badalona Museu Sant Sadurní d’Anoia Cunit Barcelona Zoo Barcelona—el Raval Barcelona Z. Universitària Caldes de Montbui la Llacuna Castellbisbal Sabadell—Parc Agrari Parets del Vallès Puig Sesolles el Prat de Llobregat Canyelles Sant Cugat CAR 2.070000 1.519167 1.676944 1.852500 2.200000 2.173889 2.103056 2.362500 2.540000 2.696944 2.431111 1.808131 1.914997 1.935642 2.123885 1.521214 1.623863 1.630325 1.907752 2.037870 1.487694 1.804796 2.377015 2.445317 1.917178 2.026470 2.302911 1.728574 2.300134 1.837509 1.467169 1.693370 2.756574 2.247574 1.794292 1.633462 2.188469 2.167751 2.105397 2.168358 1.535278 1.975463 2.069524 2.226185 2.437738 2.080219 1.721949 2.079558 41.292778 41.416944 41.330278 41.243889 41.390556 41.375000 41.523611 41.505000 41.587500 41.650833 41.760000 41.531094 41.592518 41.381968 41.418432 41.215534 41.432921 41.374910 41.288317 41.299278 41.341706 41.278610 41.517731 41.619917 41.632863 41.681290 41.747610 41.366193 41.544401 41.595390 41.271511 41.485181 41.647064 41.452149 41.433861 41.201869 41.389433 41.383899 41.379197 41.612653 41.479475 41.478924 41.565680 41.567346 41.773622 41.340456 41.288007 41.483110 4 632 177 58 6 5 146 56 74 40 936 316 343 252 411 59 415 257 573 3 545 161 81 460 421 528 1030 240 126 916 185 325 2 42 164 17 7 33 79 176 589 147 258 123 1668 8 148 158 1802.41 25,990.09 13,508.15 925.42 185.5 357.9 14,686.96 1586.01 1406.52 2814.06 22,575.42 31,647.54 31,910.67 13,153.66 6527.97 3946.43 25,729.52 19,374.48 3446.25 3275.54 18,210.02 4823.22 2228.08 8362.11 34,368.7 31,792.68 26,303.64 16,069.11 7461.61 37,477.48 11,123.03 29,463.36 198.76 620.04 21,985.71 990.45 914.09 1324.86 5269.22 18,811 32,688.25 20,531.47 19,910.4 12,099.08 23,695.86 5401.52 7951.78 14,117.55 − 1.00 241.39 243.43 206.57 225.00 − 1.00 − 1.00 150.95 71.57 158.20 182.86 150.52 90.00 160.35 230.31 26.57 161.57 90.00 237.53 − 1.00 153.43 199.80 − 1.00 209.74 315.00 155.56 221.82 78.69 270.00 0.00 180.00 101.31 45.00 135.00 65.56 203.20 135.00 0.00 191.31 45.00 165.96 71.57 135.00 180.00 233.13 − 1.00 21.80 258.69 0.00 20.88 3.73 7.45 1.18 0.00 0.00 8.58 5.27 8.98 33.37 22.02 11.67 24.78 50.89 3.73 7.91 3.33 10.87 0.00 3.73 22.14 0.00 13.44 9.43 10.07 21.25 4.25 11.67 75.00 8.33 4.25 1.18 2.36 20.14 6.35 1.18 1.67 4.25 3.54 6.87 15.81 3.54 6.67 8.33 0.00 8.98 12.75 NDVI and NDBI are strongly correlated (− 0.79) as expected. There is also some correlation between these indices and the geographical variables latitude, longitude, altitude and slope, with values for NDVI between 0.26 and 0.52 and for NDBI between − 0.48 and − 0.36. It is important to mention that the correlations between the NDVI and NDBI indices with temperatures are low, with absolute values around 0.20. Latitude shows its highest correlation with longitude (0.62) and also with altitude (0.55), continentality (0.50), NDVI (0.48) and NDBI (− 0.48). Altitude has the highest correlation with continentality (0.60) and with latitude and NDVI index (0.52). The orientation has low correlations with all the other variables, although the highest ones correspond to those of the NDVI and NDBI indices. Finally, the slope, with low correlations in general, presents the highest ones with NDVI, NDBI, continentality and altitude. Figure 5 shows some examples of the possible linear relationship between Tmin, Tmean or Tmax and some of the independent variables. The variables with the clearest linear relationship are, as expected, LSTd and LSTn according to the Air temperature in Barcelona metropolitan region from MODIS satellite and GIS data a b 12 600 11 10 9 Number of stations Number of areas 1kmx1km 500 400 300 200 8 7 6 5 4 3 100 2 1 0 0 0 200 400 600 800 1000 Altitude (m) 1200 1400 1600 1800 0 400 600 800 1000 Altitude (m) 1200 1400 1600 1800 d 0.4 0.4 0.3 0.3 Relative frequency Relative frequency c 200 0.2 0.1 0.2 0.1 0 0 0 0.1 0.2 0.3 0.4 0.5 NDVI 0.6 0.7 0.8 0.9 1 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 NDBI 0.1 0.2 0.3 0.4 0.5 Fig. 2 Histograms of altitude for all the pixels (a) and meteorological stations (b). Histograms of NDVI (c) and NDBI (d) only for meteorological stations high correlation value obtained between these variables and the temperatures. The first four RPC selected with the eigenvalue exceeding 1.0 criterion have a similar percentage of explained variance (from 23.4 to 14.6%). Whereas they explain 77.2% of data variance, the remaining eight RPC are associated with the 22.8% of data variance. Consequently, the revision of the results offered by the PCA is centred on these four first RPCs. Tables 5b and 6b show the RPC factor loadings and the explained variance by the components, for the set of thermometric stations. The first component, RPC1, is strongly related to LSTd and LSTn and cd*, explaining 23.4% of data variance. The second component, RPC2, is mainly correlated with continentality and altitude, and also related to slope, latitude and NDVI. This component explains 22.6% of data variance. RPC3 is highly correlated with longitude and latitude, and slightly negative with NDBI. The third component explains 16.6% of data variance. Finally, RPC4 is notably correlated with the sine of the orientation, moderately with the cosine of the orientation and with NDVI and NDBI. This last component explains the variance of 14.6%. It is worth mentioning that NDVI and NDBI have middle weights on the second, third and fourth components, this fact suggesting that these two variables will probably have a lesser role in the multiple regression equations than expected. Given that PCA is applied to 12 variables for the 48 thermometric stations, this relatively low number of samples in comparison with the high number of 1 km2 pixels could mask the relevance of some variables on the spatial distribution of Tmin, Tmean and Tmax. In order to detect these possible differences, RPCA has been also applied C. Serra et al. Table 3 Mean, minimum and maximum of recorded diurnal (LSTd) and nocturnal (LSTn) temperatures, normalised difference vegetation index (NDVI) and normalised difference built-up index (NDBI), latitude Pixels Mean LSTd (°C) LSTn (°C) NDBI (Lat), longitude (Lon), distance to shoreline (Con), altitude (Alt), orientation (Ori) and slope (Slp) for all 1 × 1 km pixels and for the 48 stations NDVI Lat (°) Lon (°) Con (m) Alt (m) Ori (°) Slope (%) 23.00 11.09 − 0.09 0.52 41.52 2.09 15.565.01 294.53 165.86 23.40 − 13.93 48.13 − 8.39 30.65 − 0.57 1.00 − 0.20 0.98 41.19 41.81 1.47 2.78 1.92 41.103.86 0.00 1646.00 0.00 359.26 0.00 163.91 Mean Min 23.54 − 3.35 11.47 − 5.67 − 0.04 − 0.75 0.41 0.10 41.46 41.20 2.01 1.47 13.313.41 185.50 279.17 2.00 148.30 0.00 10.66 0.00 Max 45.31 27.19 0.44 0.85 41.77 2.76 37.477.50 1668.00 360.00 75.00 Min Max Station to pixels with 1 km2 resolution. Tables 5b and 6b summarise the RPC factor loadings and the explained variance by the components, for the 4042 pixels of 1 km2. By comparing with results of Tables 5a and 6a, the first RPC is quite similar for both cases. The second RPC for Table 5a is equivalent to the third RPC for Table 5b, almost disappearing the contribution of the slope, NDVI and latitude. The third component (Table 5a), linked to latitude, longitude and NDBI, corresponds to the fourth RPC (Table 5b). The fourth RPC (Table 5a), with the contribution of NDVI, NDBI and orientation, is substituted by the fifth RPC (Table 5b), basically sinus of the orientation. Finally, the 2nd RPC (Table 5b) is constituted by NDBI, NDVI and slope. Conversely to Table 5a, where NDVI and NDBI are linked with moderate weights to more than one component, for Table 5b, these parameters are clearly related only to the second RPC. As a summary, the degrees of independence of the multiple regression variables are quite similar considering data from the set of thermometric stations or from a denser network of 1 km2 pixels. Only some discrepancies are detected comparing factor loadings corresponding to NDVI and NDBI for both spatial resolutions. Fig. 3 Evolution of Tmean with the calendar day, cd (a) and the transformed calendar day, cd* (b) 4.2 Multiple regression 4.2.1 Annual case The first multiple regression is applied to the annual case, including all days along the year 2015 accomplishing two conditions: first, only not cloudy days can be selected, given that LSTn, LSTd, NDVI and NDBI cannot be accurately computed for cloudy days; second, days with missing Tmin, Tmean or Tmax are not chosen for the multiple regression procedure. The set of days accomplishing both conditions are designed as complete data days and the same dataset with the same constraints are used at seasonal and monthly scale. Table 7 gives the different models obtained in every stepwise regression, being added one more variable until all the significant variables are used. This table also gives the R2 coefficients and RMSE. The goodness of fit improves step by step, but with minor differences. For Tmin, Tmean and Tmax in the first step with a single variable, values of R2 from 0.860 to 0.935 and RMSE from 1.8 to 2.7 °C are obtained. When all significant variables are involved, R2 coefficients of 0.920, 0.955 and 0.918 and RMSE of 1.9 °C, 1.5 °C and 2.0 °C are reached. Air temperature in Barcelona metropolitan region from MODIS satellite and GIS data Fig. 4 Dependence of Tmax on orientation, Ori (a) and sinus compound, ORIsin (b) Table 8 gives the standardised (beta weight) and nonstandardised multiple regression coefficients for the last models of every dependent variable. In the case of Tmin, the variable with the highest beta weight is LSTn, while LSTd does not appear in the multiple regression equation because it is not significant. Latitude, longitude and calendar day also have prominent standardised coefficients. For Tmean, the most important variables are LSTn and LSTd, with the rest of the variables having small or non-significant coefficients. Finally, in the case of Tmax, the variables LSTd and LSTn have the highest beta weights and are also quite similar. Latitude, longitude and altitude also play a significant role in Tmax. It should be noted that NDVI and NDBI indices have small or non-significant coefficients. Orientation and slope are the variables that present the lowest beta weight coefficients for the three temperatures, possibly because a much denser network of stations would be necessary to represent all the variety of slopes and orientations. Figure 6 shows the bar histogram of the beta coefficients of the last model, including all the significant variables. It is outstanding the positive values for LSTn and, not so relevant, the positive values of LSTd. Among the negative values, the beta weight corresponding to latitude for Tmin is the most relevant. The relationships between the temperatures obtained from the multiple regression equations and the observed Tmin, Tmean and Tmax are plotted in Fig. 7. While R2 for Tmean is 0.96, for Tmin and Tmax is 0.92. The RMSE (MAE) ranges from 1.5 Table 4 Pearson correlation coefficients for empiric temperature data, satellite data, geographic and topographic variables and transformed calendar day. The not significant coefficients (α = 0.05) are codified as ns Tmax Tmean Tmin LSTd LSTn NDBI NDVI Lat Lon Con Alt orisin oricos Slp cd* Tmax Tmean Tmin LSTd LSTn 1000 0.961 0.881 0.928 0.913 0.961 1000 0.968 0.930 0.967 0.881 0.968 1000 0.864 0.953 0.928 0.930 0.864 1000 0.907 0.913 0.967 0.953 0.907 1000 0.214 0.194 0.163 0.247 0.139 − 0.210 − 0.221 − 0.205 − 0.262 − 0.181 − 0.127 − 0.139 − 0.146 − 0.143 − 0.102 − 0.064 ns 0.070 − 0.043 0.087 − 0.071 − 0.162 − 0.227 − 0.113 − 0.203 − 0.229 − 0.200 − 0.157 − 0.204 − 0.166 0.100 0.042 − 0.018 0.066 ns 0.036 0.036 0.036 0.060 0.039 − 0.112 − 0.079 − 0.035 − 0.145 − 0.060 0.868 0.897 0.874 0.855 0.902 NDBI NDVI lat lon con alt orisin oricos slp cd* 0.214 − 0.210 − 0.127 − 0.064 − 0.071 − 0.229 0.100 0.036 − 0.112 0.868 0.194 − 0.221 − 0.139 ns − 0.162 − 0.200 0.042 0.036 − 0.079 0.897 0.163 − 0.205 − 0.146 0.070 − 0.227 − 0.157 − 0.018 0.036 − 0.035 0.874 0.247 − 0.262 − 0.143 − 0.043 − 0.113 − 0.204 0.066 0.060 − 0.145 0.855 0.139 − 0.181 − 0.102 0.087 − 0.203 − 0.166 ns ns − 0.060 0.902 1000 − 0.792 − 0.482 − 0.463 − 0.066 − 0.379 0.269 0.257 − 0.357 0.136 − 0.792 1000 0.481 0.255 0.281 0.524 − 0.300 − 0.291 0.445 − 0.103 − 0.482 0.481 1000 0.624 0.502 0.551 0.037 − 0.028 0.242 ns − 0.463 0.255 0.624 1000 − 0.331 0.034 − 0.110 ns − 0.011 ns − 0.066 0.281 0.502 − 0.331 1000 0.599 0.182 0.037 0.369 ns − 0.379 0.524 0.551 0.034 0.599 1000 − 0.281 − 0.177 0.435 ns 0.269 − 0.300 0.037 − 0.110 0.182 − 0.281 1000 0.210 − 0.187 ns 0.257 − 0.291 − 0.028 ns 0.037 − 0.177 0.210 1000 ns ns − 0.357 0.445 0.242 − 0.011 0.369 0.435 − 0.187 ns 1000 ns 0.136 − 0.103 ns ns ns ns ns ns ns 1000 C. Serra et al. Fig. 5 Dependence of Tmax and Tmin on some of the variables of the multiple regression process Air temperature in Barcelona metropolitan region from MODIS satellite and GIS data Table 5 Rotated principal components (a) RPC (stations) 1 (b) RPC (pixels) 2 3 4 1 2 3 4 5 LSTdia 0.962 − 0.103 − 0.070 0.088 0.957 − 0.194 − 0.036 − 0.036 0.001 LSTnit NDBI NDVI 0.964 0.137 − 0.136 − 0.131 − 0.343 0.552 0.031 − 0.558 0.388 − 0.022 0.546 − 0.544 0.969 0.156 − 0.153 − 0.025 − 0.851 0.878 − 0.102 − 0.033 0.160 0.014 − 0.306 0.245 0.000 0.055 − 0.046 Lat − 0.036 0.531 0.784 0.127 0.018 0.232 − 0.300 0.899 − 0.022 Lon Con 0.032 − 0.067 − 0.228 0.897 0.933 − 0.083 − 0.098 0.283 − 0.020 − 0.047 0.248 0.045 0.439 0.946 0.819 − 0.016 0.002 0.025 Alt Slp − 0.067 − 0.030 − 0.012 0.009 0.817 0.664 − 0.099 0.012 0.152 0.000 0.073 − 0.065 − 0.224 − 0.255 0.777 0.570 − 0.035 0.015 0.010 0.022 0.457 0.789 0.073 0.224 0.769 0.170 − 0.069 − 0.167 0.061 − 0.058 − 0.096 − 0.291 − 0.016 0.035 0.932 − 0.378 0.978 0.046 − 0.019 0.009 0.977 − 0.028 0.043 − 0.002 − 0.004 orisin oricos cd* (a) 48 thermometric stations (b) 4042 pixels of 1 km × 1 km (1.2) °C for Tmean to 2.0 (1.6) °C for Tmax. The histogram for Tmean residuals is also shown in this figure. Fifty-one per cent of differences between estimated and observed temperatures are lower or equal to 1.0 °C. 4.2.2 Seasonal and monthly cases Table 9 shows the R2, RMSE and the number of samples, N, with complete data for seasonal multiple regressions and for the different temperatures. The highest correlations correspond to Tmean for spring and autumn, possibly due to the moderate range of the temperatures in these seasons. The lowest R2 corresponds to the winter Tmin. The RMSE values do not exceed 2.0 °C, which corresponds to Tmax in summer. The standardised beta coefficients for each variable are given in Table 10, with only significant coefficients. The italic entries correspond to coefficients greater than 0.10, seeing at a glance the most important variables in the multiple regression equations. Tmean is the variable that depends on the minimum number of variables in any season of the year, especially in spring and autumn. For example, in autumn, only the LSTn and LSTd temperatures have notable beta coefficients, and in spring, the cd* is also important. The NDVI and NDBI indices Table 6 Total variance and percentage of variance and cumulated variance for every RPC do not have high coefficients in any of the seasons. Only NDVI has values slightly higher than 0.10 in winter and summer for Tmax, while NDBI has a negative coefficient in spring for Tmax. The geographical variables that more contribute to the multiple regression are latitude, longitude, continentality and altitude, especially for Tmax and Tmin. Slope and orientation have small or non-significant beta weights. Table 11 shows R2, RMSE and N for every monthly multiple regression and for each temperature. For monthly cases, R2 is lower than for seasonal or annual cases, ranging from 0.503 for July (T max ) to 0.867 for November (T mean ). However, RMSE have lower values than seasonal or annual cases, especially for Tmean, which ranges from 1.2 to 1.5 °C. The highest values correspond to April (Tmin) and July (Tmax), both with 2.1 °C. Figure 8 shows the RMSE for every month and Fig. 9 the estimated versus observed Tmean for November. Table 12 summarises the significant variables on the multiple regression process for every month and the different temperatures. A significant difference in comparison with annual and seasonal scale is the relevance of LST. Whereas LST is the most relevant at annual and seasonal scales, some differences are detected at monthly scale depending on the specific month and Tmin, Tmean and Tmax. The first multiple regression RPC (a) Total Variance (%) Cumulated (%) RPC (b) Total Variance (%) Cumulated (%) 1 2 2811 2712 23,427 22,600 23,427 46,027 1 2 2862 2539 23,849 21,158 23,849 45,006 3 4 1992 1753 16,600 14,606 62,627 77,232 3 4 5 1872 1736 1019 15,599 14,470 8492 60,606 75,076 83,568 (a) 48 thermometric stations. (b) 4042 pixels of 1 km × 1 km C. Serra et al. Table 7 Models obtained with the stepwise regression analysis for Tmin, Tmean and Tmax at annual scale Model Tmin 1 2 3 4 5 6 7 8 9 10 11 Model Tmean 1 2 3 4 5 6 7 8 Model Tmax 1 2 3 4 5 6 7 8 9 10 11 12 R2 RMSE Variables 0.908 0.910 0.912 0.914 0.916 0.917 0.919 0.920 0.920 0.920 0.920 2075 2047 2027 2004 1981 1968 1947 1937 1932 1930 1928 LSTn LSTn, lat LSTn, lat, cd LSTn, lat, cd, lon LSTn, lat, cd, lon, slp LSTn, lat, cd, lon, slp, NDVI LSTn, lat, cd, lon, slp, NDVI, alt LSTn, lat, cd, lon, slp, NDVI, alt, orisin LSTn, lat, cd, lon, slp, NDVI, alt, orisin, con LSTn, lat, cd, lon, slp, NDVI, alt, orisin, con, NDBI LSTn, lat, cd, lon, slp, NDVI, alt, orisin, con, NDBI, oricos .935 .950 .952 .953 .954 .955 .955 .955 1780 1551 1525 1509 1497 1484 1480 1474 LSTn LSTn, LSTd LSTn, LSTd, lon LSTn, LSTd, lon, orisin LSTn, LSTd, lon, orisin, cd LSTn, LSTd, lon, orisin, cd, alt LSTn, LSTd, lon, orisin, cd, alt, slp LSTn, LSTd, lon, orisin, cd, alt, slp, oricos .861 .890 .896 .900 .903 .913 .916 .917 .918 .918 .918 .918 2659 2365 2299 2248 2222 2096 2069 2053 2046 2044 2042 2041 LSTd LSTd, LSTn LSTd, LSTn, lon LSTd, LSTn, lon, orisin LSTd, LSTn, lon, orisin, lat LSTd, LSTn, lon, orisin, lat, alt LSTd, LSTn, lon, lat, alt, NDVI LSTd, LSTn, lon, lat, alt, NDVI, cd LSTd, LSTn, lon, lat, alt, NDVI, cd, oricos LSTd, LSTn, lon, lat, alt, NDVI, cd, oricos, orisin LSTd, LSTn, lon, lat, alt, NDVI, cd, oricos, orisin,con LSTd, LSTn, lon, lat, alt, NDVI, cd, oricos, orisin,con, slp Table 8 Standardised and not standardised multiple linear regression coefficients for the annual case. Discarded variables for the multiple regression are codified by ns Constant LSTd LSTn NDBI NDVI lat lon con alt orisin oricos slp cd Standardised coefficients Coefficients Tmin Tmax Tmin Tmean Tmax 0.414 0.473 ns 0.057 0.123 − 0.149 0.062 − 0.189 0.015 − 0.028 − 0.014 0.091 532.394 ns 0.786 1.818 − 1.481 − 13.089 5.015 5.723E−5 0.002 0.368 − 0.175 0.023 1.398 5.023 0.172 0.660 ns ns ns − 0.703 ns − 0.001 0.326 − 0.221 0.011 0.914 − 228.921 0.298 0.484 ns 2.812 5.886 − 3.341 3.839E−5 − 0.004 0.178 − 0.315 − 0.007 0.921 ns 0.801 0.027 − 0.031 − 0.286 0.234 0.096 0.083 0.033 − 0.016 0.047 0.144 Tmean 0.245 0.660 ns ns ns − 0.032 ns − 0.043 0.029 − 0.020 0.022 0.092 variable for Tmin is the latitude for 8 months, the longitude (1 month) and LSTn only for 3 months. Conversely, for Tmean, LSTn is the most relevant for 11 months. Only for August is detected a slightly higher relevance of LSTd in comparison with LSTn. Finally, for Tmax, the first multiple regression variable is the altitude for 5 months, LSTd (5 months) and LSTn (2 months). 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 Tmin Tmean Tmax Fig. 6 Histogram of the beta coefficients of the last model at an annual scale for Tmin, Tmean and Tmax Air temperature in Barcelona metropolitan region from MODIS satellite and GIS data Fig. 7 Estimated versus observed Tmin, Tmean and Tmax, and the histogram for residual Tmean at an annual scale 4.3 Spatial distribution of Tmin, Tmean and Tmax Some examples of the spatial distribution of temperatures on BMR obtained by multiple regression are shown in Figs. 10, 11, and 12. Figure 10 corresponds to T min , Table 9 Square regression coefficient, R2, root mean square error, RMSE and number of samples, N, at seasonal scale Season Variable R2 RMSE N Winter Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax 0.696 0.853 0.799 0.857 0.921 0.852 0.765 0.829 0.738 0.808 0.865 0.821 1.95 1.32 1.83 1.96 1.35 1.94 1.73 1.42 1.97 1.83 1.41 1.67 1241 1241 1241 1368 1368 1368 1146 1146 1146 1183 1183 1183 Spring Summer Autumn Tmean and Tmax obtained for November 29. The UHI phenomenon is quite evident for Tmin, being associated with Barcelona city and a neighbouring area at the south of the city with the highest temperatures. Yellow, green and blue areas represent zones of lower temperatures, which correspond to the Littoral and Pre-Littoral chains (yellow and blue areas respectively) and the Vallès valley (green area). The combined effect of the vicinity to the Mediterranean coast and the UHI phenomenon is detected on the Tmean map, with the highest temperatures along a narrow littoral fringe. It is also worth mentioning the detection of two nuclei of high temperatures, spatially coincident with those observed for Tmin. For Tmax, this effect of the vicinity to the littoral disappears and high temperatures cover a good part of the metropolitan region. Only at the northern extreme of BMR (Pre-Littoral chain) small green and blue areas are detected with lower temperatures. At seasonal scale, some examples of the spatial distribution of average Tmin for the winter season are shown in Fig. 11. Figure 12 depicts the monthly average of Tmin for winter months (January, February and March). In spite of the spatial distributions, the other three seasons are obtained with a notable degree of accuracy; the winter case C. Serra et al. Table 10 Standardised and not standardised multiple linear regression coefficients for the seasonal scale. Discarded variables for the multiple regression are codified by ns Winter Spring Summer Autumn Tmin Tmean Tmax LSTd 0.277 0.457 0.643 − 0.055 0.190 0.349 0.061 0.331 0.534 0.188 0.282 0.476 LSTn 0.641 0.625 0.409 0.712 0.611 0.458 0.497 0.301 0.068 0.709 0.663 0.452 NDBI NDVI Lat ns − 0.085 − 0.598 ns ns ns 0.062 0.121 0.334 ns − 0.052 − 0.387 − 0.053 ns ns − 0.114 ns 0.104 ns ns − 0.478 ns ns ns − 0.074 0.120 ns ns − 0.077 − 0.623 ns − 0.057 ns ns 0.052 0.284 Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax Lon 0.457 ns − 0.264 0.320 ns − 0.221 0.455 ns ns 0.465 ns − 0.350 Con Alt 0.305 0.128 0.115 − 0.119 ns − 0.412 0.129 0.101 ns − 0.088 0.167 − 0.206 ns 0.060 ns − 0.193 0.201 − 0.323 0.256 0.237 ns ns ns − 0.390 Slp orisin oricos cd* 0.077 0.064 − 0.052 0.049 ns − 0.054 ns − 0.038 ns 0.063 0.039 ns ns 0.044 − 0.021 ns 0.065 − 0.048 0.117 0.043 ns 0.061 0.056 − 0.030 ns 0.067 − 0.064 0.085 0.059 ns 0.075 0.036 ns − 0.036 ns ns − 0.131 − 0.138 − 0.181 0.234 0.213 0.176 0.209 0.310 0.297 ns ns − 0.040 Month Variable R2 RMSE N January Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax Tmin Tmean Tmax 0.671 0.815 0.726 0.664 0.824 0.774 0.707 0.824 0.788 0.648 0.769 0.692 0.517 0.651 0.604 0.692 0.748 0.605 0.606 0.606 0.503 0.665 0.727 0.666 0.735 0.779 0.746 0.706 0.765 0.797 0.825 0.867 0.822 0.669 0.728 0.721 1.86 1.24 1.72 1.84 1.24 1.72 1.96 1.34 1.92 2.06 1.32 1.68 1.77 1.20 1.76 1.67 1.21 1.90 1.66 1.29 2.05 1.67 1.25 1.60 1.51 1.22 1.40 1.74 1.45 1.40 1.81 1.37 1.70 1.76 1.28 1.22 609 609 609 356 356 356 276 276 276 393 393 393 455 455 455 520 520 520 504 504 504 423 423 423 219 219 219 301 301 301 615 615 615 267 267 267 February March April May June July August September October November December for Tmin is introduced here given that the characteristics of these four maps clearly manifest the UHI phenomenon in the BMR. The map of average winter Tmin reproduces the two nuclei of UHI on the downtown of the Barcelona city and at the south along the coast. If the spatial analysis is revised at a monthly scale, very similar spatial patterns to those observed for November 29 (Tmin) are now found for the cold months of January and February. In the case of March (a more temperate month), the UHI is not so clear, but the two nuclei appear again. 3.0 Tmin Tmean Tmax 2.5 2.0 RMSE (oC) Table 11 Square regression coefficient, R2, root mean square error, RMSE and number of samples, N, at monthly scale 1.5 1.0 0.5 0.0 1 2 3 4 Fig. 8 RMSE for every month 5 6 7 Month 8 9 10 11 12 Air temperature in Barcelona metropolitan region from MODIS satellite and GIS data 20 4.4 Discussion of the results November Estimated Tmean (ºC) 16 12 8 4 0 0 4 8 12 Observed Tmean (ºC) Fig. 9 Estimated versus observed Tmean for November Table 12 Significant variables in decreasing order exceeding Beta equal to 0.20 for the multiple regression process at monthly scale of Tmin, Tmean and Tmax Month January February March April May June July August September October November December Month January February March April May June July August September October November December Month January February March April May June July August September October November December 16 20 With respect to the correlations among dependent and independent variables used in the multiple regression, it is noticeable, as expected, the high correlation (0.97) between the three daily temperatures (minimum, mean and maximum) and the two LST temperatures and also the calendar day. The correlations are notably small for the rest of independent variables, sometimes achieving values lower to 0.23. Another relevant question is that results obtained by PC analysis, taking into account the thermometric stations or the dense network of 1 km2 pixels, are very similar. Consequently, the relatively sparse distribution of the 48 thermometric gauges would not be a shortcoming to obtain a reliable spatial distribution of temperatures, being then defined a relatively good image of the thermometric variability on the BMR. With respect to the multiple linear processes, the square regression coefficients obtained at annual scale are notably good, in spite of the RMSE varies within the (1.5–2.0 °C), results quite similar to those obtained by Cristóbal et al. (2008) Tmin significant standardised beta coefficients > 0.20 lat(− 0.926), lon (0.754), con(0.524), LSTn(0.486), LSTd (0.358), alt(0.254) lat(− 0.723), LSTn(0.690), lon(0.456), con(0.425) LSTn(0.661), LSTd(0.205) LSTn(0.701), lat(− 0.683), lon(0.531), con(0.626) lat(− 0.709), lon(0.686), LSTn(0.361), NDBI(0.231) lat(− 0.564), lon(0.533), LSTn(0.499), NDVI (− 0.242), alt (0.236), cd*(0.231), slp(0.217) lon(0.628), lat(− 0.582), LSTn(0.467), NDVI(− 0.207) lat(− 0.600), lon(0.476), LSTn(0.474) lat(− 0.893), lon(0.835), LSTn(0.439) LSTn(0.690), lat(− 0.569), lon(0.462) lat(− 0.770), LSTn(0.600), lon(0.584), con(0.311), alt(0.270) lat(− 0.598), LSTn(0.531), alt(0.309), NDVI(− 0.238) Tmean significant standardised beta coefficients > 0.20 LSTn(0.588), LSTd (0.457), con(0.213) LSTn(0.482), alt(− 0.323),LSTd (0.286) LSTn(0.650), LSTd(0.393) LSTn(0.685) LSTn(0.376),alt(− 0.366), LSTd (0.250), cd*(− 0.202) LSTn(0.451), LSTd (0.409), NDBI(− 0.219) LSTn(0.435), alt(− 0.333), cd*(− 0.318), LSTd (0.208) LSTd(0.358), LSTn (0.282), alt(− 0.273) LSTn(0.369), alt(− 0.377), LSTd (0.265), NDVI(0.212) LSTn(0.629), LSTd (0.219) LSTn(0.645), LSTd (0.223) LSTn(0.577), LSTd (0.271), slp(0.227) Tmax significant standardised beta coefficients > 0.20 LSTd(0.530), alt(− 0.532), LSTn (0.398), lat(0.484), lon(− 0.309) LSTd(0.557), alt(− 0.547), lat (0.357), lon(− 0.322) LSTd(0.556), LSTn(0.489), con(0.329), alt(− 0.317) alt(− 0.459), lon(− 0.490), LSTn(0.415), lat(0.360), LSTd(0.260) alt(− 0.409), LSTd(0.395), con(0.387), cd*(− 0.313) LSTd(0.510), alt(− 0.422), con (0.330), NDBI(− 0.329) alt(− 0.403), LSTd(0.399), con(0.390), cd*(− 0.347) LSTd(0.535), alt(− 0.468), con (0.243) alt(− 0.604), LSTd(0.456), NDVI(0.364) LSTn(0.503), LSTd(0.363), lon(− 0.345), alt(− 0.329), lat(0.275) LSTn(0.434), LSTd(0.430), alt(− 0.410), alt(− 0.329), lon(− 0.359), lat(0.284) alt(− 0.814), lon(− 0.508), lat(0.481), LStd(0.355) C. Serra et al. Fig. 10 Example of the spatial distribution of Tmin, Tmean and Tmax derived by multiple regression for November 29, 2015 Temperature (ºC) Tmin 29/11/2015 -6 to -5 -5 to -4 -4 to -3 -3 to -2 -2 to -1 -1 to 0 0 to 1 1 to 2 2 to 3 3 to 4 4 to 5 5 to 6 6 to 7 7 to 8 8 to 10 Temperature (ºC) Tmean 29/11/2015 2 to 3 3 to 4 4 to 5 5 to 6 6 to 7 7 to 8 8 to 9 9 to 10 10 to 11 11 to 12 Temperature (ºC) Tmax 29/11/2015 for the whole Catalonia. For the results at seasonal and monthly scales, even though the square regression coefficients are lesser than those obtained at an annual scale, the RMSE values are very similar, being not exceeded 2.1 °C. Consequently, the images of the spatial distribution of temperatures should be of similar quality at annual, seasonal and monthly scales. Nevertheless, the highest RMSE values are detected at a monthly scale (Tmin for April and Tmax for July). 7 to 8 8 to 9 9 to 10 10 to 11 11 to 12 12 to 13 13 to 14 14 to 15 15 to 16 16 to 17 17 to 18 Another noticeable characteristic is the low weight of NDVI and NDBI on the multiple linear regression, which has been also detected in other similar analysis around the world. In spite of both coefficients could be relevant, as they represent the type of vegetation and building respectively, the LST obtained from satellite data could itself include a great percentage of the information concerning NDVI and NDBI. Air temperature in Barcelona metropolitan region from MODIS satellite and GIS data Tmin Winter 2015 Fig. 11 Winter season spatial distribution of average Tmin derived by multiple regression Temperature (ºC) -6 to -5 -5 to -4 -4 to -3 -3 to -2 -2 to -1 -1 to 0 0 to 1 1 to 2 2 to 3 3 to 4 4 to 5 5 to 6 6 to 7 7 to 8 8 to 9 9 to 10 From an applied point of view, given that the obtained thermometric maps are submitted to a maximum RMSE of 2 °C, a dense network of minimum, mean and maximum temperature data could be possible to analyse thermometric phenomena (UHI and hot and cold outbreaks) affecting life quality and health of BMR population. Additionally, data obtained with smaller pixels, for instance with LANDSAT satellite, would permit a notable increase in the spatial resolution of temperatures. Unfortunately, the available data from LANDSAT is nowadays minor than that found from MODIS satellite, and the accuracy of the results would be then questionable. 5 Conclusions Previous to the multiple regression process, the Pearson correlation coefficient and the PCA have permitted to detect links between empiric temperatures recorded at 48 meteorological stations and satellite data, geographic and topographic data and transformed calendar days. The PCA has also permitted to validate if the set of the thermometric stations are appropriate for a good multiple regression process by comparing RPCs and factor loadings corresponding to 48 stations dataset and 1 km2 pixel network. In spite of the very different spatial data density for thermometric stations and pixel coverage, a few discrepancies are found with respect to the factor loadings of NDVI and NDBI. In this way, a denser network of thermometric data would improve the role of NDVI and NDBI on the multiple regression. It is also worth mentioning the substitution of the topographic parameter of orientation (Ori) by cos(Ori) and sin(Ori), being detected a slight improvement on the relevance of sin(Ori) when it is used instead of Ori in the multiple regression process. A revision of the multiple regression analyses results manifests the strong relevance of LSTn for Tmin and Tmean and LSTd for Tmax at the annual scale, as obtained by Thanh et al. (2016). A similar pattern is observed at a seasonal scale. The relevance at a monthly scale of LSTd on Tmin is not significant for February, April, May, June and September. Additionally, the relevance of LSTn on Tmax is not significant for September. With respect to specific results at the annual scale, first of all, it is noticeable that the best fit between empiric temperatures and those generated by multiple regression is usually found for Tmean, being obtained the worst for Tmax in terms of R2 and RMSE. It is also noticeable that cd* only plays a relatively important role for Tmin and the set of relevant variables are not the same for Tmin, Tmean and Tmax. At monthly scale, whereas the best fits are effectively obtained since January to December for Tmean, the worst fit is obtained for Tmax (June, July) and Tmin (April). At a seasonal scale, the results are quite different. Whereas the minimum residual for Tmean is detected in winter and spring, for Tmin it is found in summer and for Tmax in autumn. In short, the reasonably good results of the multiple regression process would permit: – – Describing with detail (pixels of 1 km2) the spatial distribution of temperatures, notably improving the spatial data density on BMR derived from the thermometric network and without applying interpolations. Obtaining detailed maps of UHI phenomenon on urban areas. It has to be remembered that these details of the UHI intensity could not be obtained from a few thermometric stations. In particular, two clear focus of high UHI intensity for T min in winter have been detected. Additionally, the smooth temperatures along the Mediterranean coast are verified by observing Tmax and Tmean maps. A systematic and detailed spatial description of temperatures with the methodology used in this paper could be a significant improvement in the analysis of cold and hot outbreaks. It should be underlined that these analyses could be C. Serra et al. Fig. 12 Monthly spatial distribution (winter season) of average Tmin derived by multiple regression Tmin January 2015 Temperature (ºC) -6 to -5 -5 to -4 -4 to -3 -3 to -2 -2 to -1 -1 to 0 0 to 1 1 to 2 2 to 3 3 to 4 4 to 5 5 to 6 6 to 7 7 to 8 8 to 10 Tmin February 2015 Temperature (ºC) -6 to -5 -5 to -4 -4 to -3 -3 to -2 -2 to -1 -1 to 0 0 to 1 1 to 2 2 to 3 3 to 4 4 to 5 5 to 6 6 to 7 7 to 8 8 to 10 Tmin March 2015 Temperature (ºC) -6 to -5 -5 to -4 -4 to -3 -3 to -2 -2 to -1 -1 to 0 0 to 1 1 to 2 2 to 3 3 to 4 4 to 5 5 to 6 6 to 7 7 to 8 8 to 9 9 to 11 very useful to study the effects on the life quality and health of the Barcelona city and metropolitan area inhabitants. Conversely to a relatively scarce distribution of thermometric stations used to analyse these outbreaks, the multiple linear regression method provides a more detailed (1 km × 1 km) spatial distribution of air temperatures. Finally, comparing R2 and RMSE for BMR with those obtained by Cristóbal et al. (2008) for the whole Catalonia, it is observed that better results have been obtained for BMR. Whereas for Tmin and Tmax, the best results are achieved for BMR (R2 equal to 0.92 in front of 0.54–0.57 and RMSE equal to 1.5–1.6 °C in front of 2.3–1.8 °C), for Tmean, R2 is again better for BMR (0.96 in front of 0.66). The RMSE for Tmean obtained for the whole Catalonia (1.3 °C) is slightly better than that obtained for BMR (1.5 °C). One reason for these differences could be that BMR area is ten times smaller than Catalonia, implying a minor variability of temperatures and geographical characteristics and permitting a better Air temperature in Barcelona metropolitan region from MODIS satellite and GIS data description of the spatial distribution of temperatures. A higher spatial density of thermometric stations could be another factor favouring the BMR results. Acknowledgements Temperature data were gently provided by the Agencia Estatal de Meteorología, AEMET, and Servei Meteorològic de Catalunya, SMC. Funding information This research has been supported by the Spanish Government through the project BIA2015-68623-R (Ministerio de Economía y Competitividad, MINECO and European Regional Development Fund (ERDF). References Arnfield AJ (2003) Two decades of urban climate research: a review of turbulence, exchanges of energy and water, and the urban heat island. Int J Climatol 23:1–26 Benali A, Carvalho AC, Nunes JP, Carvalhais N, Santos A (2012) Estimating air surface temperature in Portugal using MODIS LST data. Remote Sens Environ 124:108–121 Buishand TA (1982) Some methods for testing the homogeneity of rainfall records. J Hydrol 58:11–27 Bustos E, Meza FJ (2015) A method to estimate maximum and minimum air temperature using MODIS surface temperature and vegetation data: application to the Maipo Basin, Chile. Theor Appl Climatol 120:211–226. https://doi.org/10.1007/s00704-014-1167-2 Caselles V, Lopez Garcia MJ, Melia J, Perez Cueva AJ (1991) Analysis of the heat-island effect of the city of Valencia, Spain, through air temperature transects and NOAA satellite data. Theor Appl Climatol 43:195–203 Cheval S, Dumitrescu A (2009) The July urban heat island of Bucharest as derived from Modis images. Theor Appl Climatol 96:145–153. https://doi.org/10.1007/s00704-008-0019-3 Cheval S, Dumitrescu A, Bell A (2009) The urban heat island of Bucharest during the extreme high temperatures of July 2007. Theor Appl Climatol 97:391–401. https://doi.org/10.1007/s00704008-0088-3 Coll C, Caselles V, Galve JM, Valor E, Raquel NR, Sánchez JM, Rivas R (2005) Ground measurements for the validation of land surface temperatures derived from AATSR and MODIS data. Remote Sens Environ 97:288–300 Cresswell MP, Morse AP, Thomson MC, Connor SJ (1999) Estimating surface air temperatures, from Meteosat land surface temperatures, using an empirical solar zenith angle model. Int J Remote Sens 20(6):1125–1132. https://doi.org/10.1080/014311699212885 Cristóbal J, Ninyerola M, Pons X (2008) Modeling air temperature through a combination of remote sensing and GIS data. J Geophys Res 113:D13106. https://doi.org/10.1029/2007JD009318 Dash P, Gottsche FM, Olesen FS, Fischer H (2001) Retrieval of land surface temperature and emissivity from satellite data: physics, theoretical limitations and current methods. J Indian Soc Remote Sens 29:23–30 Dash P, Göttsche FM, Olesen FS, Fischer H (2002) Land surface temperature and emissivity estimation from passive sensor data: theory and practice-current trends. Int J Remote Sens 23(13):2563–2594. https://doi.org/10.1080/01431160110115041 Didari S, Norouzi H, Zand-Parsa S, Khanbilvardi R (2017) Estimation of daily minimum land surface air temperature using MODIS data in southern Iran. Theor Appl Climatol 130:1149–1161. https://doi.org/ 10.1007/s00704-016-1945-0 Fabrizi R, De Santis A, Gómez A (2011) Satellite und ground based sensors for the urban heat island analysis in the city of Madrid. Joint Urban Remote Sensing Event, Munich, pp 349–352 Fick SE, Hijmans RJ (2017) WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int J Climatol 37:4302–4315 Fu G, Shen Z, Zhang X, Shi P, Zhang Y, Wu J (2011) Estimating air temperature of an alpine meadow on the Northern Tibetan Plateau using MODIS land surface temperature. Acta Ecol Sin 31:8–13 Gallo KP, Tarpley JD, McNab AL, Karl TR (1995) Assessment of urban heat islands: a satellite perspective. Atmos Res 37:37–43 Giannaros TM, Melas D (2012) Study of the urban heat island in a coastal Mediterranean City: the case study of Thessaloniki, Greece. Atmos Res 118:103–120 Good E (2015) Daily minimum and maximum surface air temperatures from geostationary satellite data. J Geophys Res Atmos 120:2306– 2324. https://doi.org/10.1002/2014JD022438 Harrell FE (2001) Regression modelling strategies, With applications to linear models, logistic regression models, and survival analysis. Springer, New York, 568 pp Hengl T, Heuvelink GBM, Percec TM, Pebesma EJ (2012) Spatiotemporal prediction of daily temperatures using time-series of MODIS LST images. Theor Appl Climatol 107:265–277. https:// doi.org/10.1007/s00704-011-0464-2 Hung T, Uchihama D, Ochi S, Yasuoka Y (2006) Assessment with satellite data of the urban heat island effects in Asian mega cities. Int J Appl Earth Obs Geoinf 8:34–48 Imhoff ML, Zhang P, Wolfe RE, Bounoua L (2010) Remote sensing of the urban heat island effect across biomes in the continental USA. Remote Sens Environ 114:504–513 Janatian N, Sadeghi M, Sanaeinejad SH, Bakhshian E, Farid A, Hasheminiaa SM, Ghazanfari S (2017) A statistical framework for estimating air temperature using MODIS land surface temperature data. Int J Climatol 37:1181–1194 Jang K, Kang S, Kimball JS, Hong SY (2014) Retrievals of all-weather daily air temperature using MODIS and AMSR-E data. Remote Sens 6:8387–8404. https://doi.org/10.3390/rs6098387 Jin M, Dickinson RE, Zhang DL (2005) The footprint of urban areas on global climate as characterized by MODIS. J Clim 18:1551–1565 Jolliffe IT (1986) Principal component analysis. In: Springer series in statistics. Springer, New York 271 pp Kilibarda M, Hengl T, Heuvelink GBM, Gräler B, Pebesma E, PercecTadic M, Bajat N (2014) Spatio-temporal interpolation of daily temperatures for global land areas at 1 km resolution. J Geophys Res Atmos 119:2294–2313. https://doi.org/10.1002/ 2013JD020803 Kloog I, Nordio F, Lepeule J, Andrea PA, Lee M, Auffray A, Schwartz J (2017) Modelling spatio-temporally resolved air temperature across the complex geo-climate area of France using satellite-derived land surface temperature data. Int J Climatol 37:296–304 Landsberg HE (1981) The urban climate. International Geophysics Series, 28. Academic Press, New York 275 pp Lin S, Moore NJ, Messina JP, DeVisser MH, Wu J (2012) Evaluation of estimating daily maximum and minimum air temperature with MODIS data in east Africa. Int J Appl Earth Obs Geoinf 18:128–140 Lin X, Zhang W, Huang Y, Sun W, Han P, Yu L, Sun F (2016) Empirical estimation of near-surface air temperature in China from MODIS LST data by considering physiographic features. Remote Sens 8: 629. https://doi.org/10.3390/rs8080629 Ma W, Zhou L, Zhang H, Zhang Y, Dai X (2016) Air temperature field distribution estimations over a Chinese mega-city using MODIS land surface temperature data: the case of Shanghai. Front Earth Sci 10:38–48. https://doi.org/10.1007/s11707-015-0510-y Martínez MD, Serra C, Burgueño A, Lana X (2010) Time trends of daily maximum and minimum temperatures in Catalonia (ne Spain) for the period 1975–2004. Int J Climatol 30:267–290 C. Serra et al. Mendelsohn R, Kurukulasuriya P, Basist A, Kogan F, Williams C (2007) Climate analysis with satellite versus weather station data. Clim Chang 81:71–83. https://doi.org/10.1007/s10584-006-9139-x Miao S, Chen F, Lemone MA, Tewari M, Li Q, Wang Y (2009) An observational and modeling study of characteristics of urban heat island and boundary layer structures in Beijing. J Appl Meteorol Climatol 48:484–501 Moreno-Garcia MC (1994) Intensity and form of the urban heat island in Barcelona. Int J Climatol 14:705–710 Nemani RR, Running SW (1989) Estimation of regional surface resistance to evapotranspiration from NDVI and thermal IR AVHRR data. J Appl Meteorol 28:276–284 Nieto H, Sandholt I, Aguado I, Chuvieco E, Stisen S (2011) Air temperature estimation with MSG-SEVIRI data: calibration and validation of the TVX algorithm for the Iberian Peninsula. Remote Sens Environ 115:107–116 Peón J, Recondo C, Calleja JF (2014) Improvements in the estimation of daily minimum air temperature in peninsular Spain using MODIS land surface temperature. Int J Remote Sens 35:5148–5166 Pettitt AN (1979) A non-parametric approach to the change-point problem. J Appl Stat 28(2):126–135 Pongrácz R, Bartholy J, Dezsö Z (2006) Remotely sensed thermal information applied to urban climate analysis. Adv Space Res 37:2191– 2196 Pongrácz R, Bartholy J, Dezsö Z (2010) Application of remotely sensed thermal information to urban climatology of Central European cities. Phys Chem Earth 35:95–99 Preisendorfer RW (1988) Principal component analysis in meteorology and oceanograph. Elsevier, New York Prihodko L, Goward SN (1997) Estimation of air temperature sensed surface observations from remotely. Remote Sens Environ 60: 335–346 Rao PK (1972) Remote sensing of urban “heat islands” from an environmental satellite. Bull Am Meteorol Soc 53:647–648 Recondo C, Peón JJ, Zapico E, Pendás E (2013) Empirical models for estimating daily surface water vapour pressure, air temperature, and humidity using MODIS and spatiotemporal variables. Applications to peninsular Spain. Int J Remote Sens 34:8051–8080. https://doi. org/10.1080/01431161.2013.828185 Richman RB (1986) Rotation of principal components. Int J Climatol 6: 293–335. https://doi.org/10.1002/joc.3370060305 Schneider A, Friedl MA, Potere D (2009) A new map of global urban extent from MODIS satellite data. Environ Res Lett 4:044003. https://doi.org/10.1088/1748-9326/4/4/044003 Shah DB, Pandya MR, Trivedi HJ, Jani AR (2013) Estimating minimum and maximum air temperature using MODIS data over IndoGangetic Plain. J Earth Syst Sci 122:1593–1605 Shen S, Leptoukh GG (2011) Estimation of surface air temperature over central and eastern Eurasia from MODIS land surface temperature. Environ Res Lett 6:045206 Sohrabinia M, Zawar-Reza P, Rack W (2015) Spatio-temporal analysis of the relationship between LST from MODIS and air temperature in New Zealand. Theor Appl Climatol 119:567–583 Stewart ID (2011) A systematic review and scientific critique of methodology in modern urban heat island literature. Int J Climatol 31:200– 217 Sun H, Chen Y, Gong A, Zhao X, Zhan W, Wang M (2014) Estimating mean air temperature using MODIS day and night land surface temperatures. Theor Appl Climatol 118:81–92 Thanh NP, Kappas M, Degener J (2016) Estimating daily maximum and minimum land air surface temperature using MODIS land surface temperature data and ground truth data in northern Vietnam. Remote Sens 8:1002. https://doi.org/10.3390/rs8121002 Tomlinson CJ, Chapman L, Thornes JE, Baker C (2011) Remote sensing land surface temperature for meteorology and climatology: a review. Meteorol Appl 18:296–306 Tomlinson CJ, Chapman L, Thornes JE, Baker CJ, Prieto-Lopez T (2012) Comparing night-time satellite land surface temperature from MODIS and ground measured air temperature across a conurbation. Remote Sens Lett 3:657–666. https://doi.org/10.1080/01431161. 2012.659354 Vancutsem C, Ceccato P, Dinku T, Connor SJ (2010) Evaluation of MODIS land surface temperature data to estimate air temperature in different ecosystems over Africa. Remote Sens Environ 114:449– 465 Voogt JA, Oke TR (2003) Thermal remote sensing of urban climates. Remote Sens Environ 86:370–384 Wan ZM, Zhang YL, Zhang QC, Li ZL (2002) Validation of the landsurface temperature products retrieved from Terra Moderate Resolution Imaging Spectroradiometer data. Remote Sens Environ 83:163–180 Wan Z, Wang P, Li X (2004) Using MODIS land surface temperature and normalized difference vegetation index products for monitoring drought in the southern Great Plains, USA. Int J Remote Sens 25(1):61–72 Wang W, Liang S, Meyers T (2008) Validating MODIS land surface temperature products using long-term nighttime ground measurements. Remote Sens Environ 112:623–635 Wang L, Koike T, Yang K, Yeh PJF (2009) Assessment of a distributed biosphere hydrological model against streamflow and MODIS land surface temperature in the upper Tone River Basin. J Hydrol 377: 21–34 Wenbin Z, Aifeng L, Shaofeng J (2013) Estimation of daily maximum and minimum air temperature using MODIS land surface temperature products. Remote Sens Environ 130:62–73 Wijngaard JB, Klein Tank MG, Können GP (2003) Homogeneity of 20th century European daily temperature and precipitation series. Int J Climatol 23:679–692 Williamson SN, Hik DS, Gamon JA, Kavanaugh JL, Koh S (2013) Evaluating cloud contamination in clear-sky MODIS Terra daytime land surface temperatures using ground-based meteorology station observations. J Clim 26:1551–1560 Xu Y, Knudby A, Ho HC (2014) Estimating daily maximum air temperature from MODIS in British Columbia, Canada. Int J Remote Sens 35:8108–8121. https://doi.org/10.1080/01431161.2014.978957 Yang P, Ren G, Liu W (2013) Spatial and temporal characteristics of Beijing urban heat island intensity. J Appl Meteorol Climatol 52: 1803–1816 Yang YZ, Cai WH, Yang J (2017) Evaluation of MODIS land surface temperature data to estimate near-surface air temperature in Northeast China. Remote Sens 9:410. https://doi.org/10.3390/ rs9050410 Yasuoka Y (2006) Assessment with satellite data of the urban heat island effects in Asian mega cities. Int J Appl Earth Obs Geoinf 8:34–48 Zaksek K, Schroedter-Homscheidt M (2009) Parameterization of air temperature in high temporal and spatial resolution from a combination of the SEVIRI and MODIS instruments. ISPRS J Photogramm Remote Sens 64(4):414–421 Zhang X, Friedl MA, Schaaf CB, Strahler AH, Hodges JCF, Gao F, Reed BC, Huete A (2003) Monitoring vegetation phenology using MODIS. Remote Sens Environ 84(3):471–475 Zhang W, Huang Y, Yu Y, Sun W (2011) Empirical models for estimating daily maximum, minimum and mean air temperatures with MODIS land surface temperatures. Int J Remote Sens 32(24):9415–9440. https://doi.org/10.1080/01431161.2011.560622 Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.