Improving Mean Minimum and Maximum Month-to-Month Air Temperature Surfaces Using Satellite-Derived Land Surface Temperature

Abstract

Month-to-month air temperature (T_air) surfaces are increasingly demanded to feed quantitative models related to a wide range of fields, such as hydrology, ecology or climate change studies. Geostatistical interpolation techniques provide such continuous and objective surfaces of climate variables, while the use of remote sensing data may improve the estimates, especially when temporal resolution is detailed enough. The main goal of this study is to propose an empirical methodology for improving the month-to-month T_air mapping (minimum and maximum) using satellite land surface temperatures (LST) besides of meteorological data and geographic information. The methodology consists on multiple regression analysis combined with the spatial interpolation of residual errors using the inverse distance weighting. A leave-one-out cross-validation procedure has been included in order to compare predicted with observed values. Different operational daytime and nighttime LST products corresponding to the four months more characteristic of the seasonal dynamics of a Mediterranean climate have been considered for a thirteen-year period. The results can be considered operational given the feasibility of the models employed (linear dependence on predictors that are nowadays easily available), the robustness of the leave-one-out cross-validation procedure and the improvement in accuracy achieved when compared to classical T_air modeling results. Unlike what is considered by most studies, it is shown that nighttime LST provides a good proxy not only for minimum T_air, but also for maximum T_air. The improvement achieved by the inclusion of remote sensing LST products was higher for minimum T_air (up to 0.35 K on December), especially over forests and rugged lands. Results are really encouraging, as there are generally few meteorological stations in zones with these characteristics, clearly showing the usefulness of remote sensing to improve information about areas that are difficult to access or simply with a poor availability of conventional meteorological data.

1. Introduction

The surface air temperature (hereafter T_air) is defined as the temperature measured by a thermometer exposed to the air in a place protected from direct solar radiation, normally located at about 1.5 m above surface level on land [1]. The T_air is a key climatic and meteorological variable that makes it possible to quantify biophysical processes at surface level, such as energy flows, actual and potential evapotranspiration, water stress, and species distribution [2,3,4,5,6]. Further, it is used in many environmental applications, including vector-borne disease bionomics [7], terrestrial hydrology [8], biosphere processes [9] and climate change [10]. In addition to mean T_air, minimum and maximum values are central to understanding these biosphere processes, and are often preferred measurements.

T_air is measured at meteorological stations, which provide accurate (although data filtering and quality control are very important) temporally discrete information that is spatially characterized by the density and distribution of the network of available stations. Hence, these stations have limited ability to describe the spatial heterogeneity of T_air over large areas of the Earth, especially in regions with complex topographies and landscapes. Therefore, there are many studies on climate interpolation methods using meteorological data that may or may not include geographical variables as predictors.

Using information obtained from airborne sensors or artificial satellites (remote sensing) to observe the Earth from space is the only methodology that effectively evaluates the spatial distribution of land surface variables on a regional and global scale. Thanks to the technological developments in the latest generation of spectral band sensors in the thermal infrared region (8–14 μm) (e.g., Sentinel-3, Landsat, Moderate Resolution Imaging Spectroradiometer (MODIS), Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), Advanced Along Track Scanning Radiometer (AATSR), and Advanced Very High Resolution Radiometer (AVHRR)), the capacity of existing Earth observation programs has improved. These sensors incorporate new spectral measurement channels and provide much better spatial and spectral resolution than just a few decades ago. Likewise, the long historical satellite data series that exist today (Landsat, over 40 years; NOAA-AVHRR, over 30 years; MODIS, 16 years) make it possible to combine climate cartography based on data from meteorological stations with satellite information [11,12,13,14,15].

Land surface temperature (LST) retrieved from thermal infrared remote sensing measurements has been frequently used to estimate T_air [13,14,16,17,18,19,20,21,22,23]. The LST is closely related to the surface energy balance and the water status of the Earth’s surface. It mainly depends on the amount of radiative energy that the surface absorbs, the partitioning of heat in sensible and latent heat fluxes, and the characteristics of the atmosphere close to the ground (in particular T_air and turbulence) [24].

The level of precision generally accepted as “accurate” for T_air estimates based on remote sensing is between 1 and 2 °C [25]. A review can be found in [26] of the methods commonly used to estimate T_air based on LST. These authors divided the methods into three distinct groups: (1) statistical approaches based on regression techniques [21,27,28,29]; (2) temperature-vegetation index [4,30]; and (3) energy-balance approaches [12,31]. The methodology applied in our study is based on geostatistical interpolation techniques and we combine long meteorological series of T_air with geographical variables, and LST remote sensing estimates. This methodology has been extensively tested in previous studies [32,33,34] that model Iberian Peninsula T_air and precipitation for monthly and annual periods using only geographical variables. There are also many works that use LST to predict T_air, considering daily or instantaneous T_air for short time scales (generally one year) [4,12,20,21,23,26,27,35,36,37,38,39]. Nevertheless, here we will focus on estimating monthly T_air (similarly to [40]) for a long time series (similarly to recently published works by [41,42,43]).

The present work has the general aim of progressing in the understanding and prediction of future environmental changes. The main objective is to improve the monthly T_air spatial estimates (minimum and maximum) currently provided by Geographic Information Systems (GIS)-based models that combine statistical (multiple regression) and spatial interpolation of the regression residuals. T_air data have been obtained from meteorological stations but we have added LST from satellite data as a predictor. In this study, we use a long and robust time series (from 2003 to 2015) that includes data from both conventional meteorological observations and remote sensing as well as geographic factors such as altitude, latitude and continentality, etc., and that also takes into account LST products from different satellites. Our objective is to demonstrate that long-term LST averages could be useful to improve long-term averaged climate maps. Some previous works [4,12,20,21,23,26,27,35,36,37,38,39] incorporate instantaneous, hourly or daily LST values to build detailed time-resolution climate surfaces, and LST information is very useful in these scenarios because the geographic predictors are less important in a daily or sub-daily basis where specific atmospheric states could be more significant to quantify temperature.

We have focussed particularly on detecting when and in which cases LST could improve the GIS-based maps, especially for minimum and maximum T_air for which it is not enough to use just DEM-related variables to obtain reliable digital temperature surfaces. Changes in the T_air estimates at the spatial and temporal levels were evaluated considering the analysis of data from the different times of the year as well as the differentiation between land cover and levels of orographic complexity of the territory.

2. Study Area

The study area corresponds to Catalonia, northeastern Iberian Peninsula on the Mediterranean coast, centered at 1°40′ of east longitude and 41°44′ of north latitude and with an approximate area of 32,106 km² (see Figure 1). This region is characterized by a heterogeneous orography with elevations ranging from sea level to about 3000 m. It has variable water availability both intra-annually and inter-annually due to repeated drought periods. The climate is typically Mediterranean. The Pyrenees and the neighboring areas have a high-altitude climate, with abundant snow and minimum temperatures below 0 °C in winter, and an annual rainfall above 1000 mm. Along the coast, the climate is mild and temperate with temperatures increasing from north to south, while the rain follows the opposite gradient. The inland areas have a continental Mediterranean climate, with colder winters and hotter summers than coastal areas [44].

Catalonia is composed by four main vegetation units. The first is Mediterranean vegetation, which covers most of the surface area of Catalonia and is dominated by sclerophyllous evergreen forests ranging from sea level to 1000 m. The second is Eurosiberian vegetation dominated by deciduous broad-leaved forests ranging from 1000 to 1600 m. The third is subalpine vegetation dominated by needle-leaf evergreen conifers from 1600 to 2200 m, and the fourth is the Boreoalpine vegetation dominated by grasslands in areas higher than 2200 m [45]. Almost 38% of Catalonia’s surface is covered by forests (including evergreen forests, deciduous broad-leaved forests and mixed forests). Another important category is crops, which represent 34.4% of the surface. The third most abundant cover is shrublands (16.4%) followed by grassland (4%), urban (3%), inland water (1.6%), and other minor categories such as roads, bare soils, etc. (2.6%) [46].

3. Materials

The common temporal range used in the study, in accordance with the availability of satellite data, is from 2003 to 2011, while it extended to 2015 in the case of MODIS data. We selected four months to analyze the seasonal dynamics that correspond to a Mediterranean climate (detailed in the previous section), corresponding to the months of the equinoxes and solstices: March, June, September and December. The different datasets used are explained below and summarized in

Table 1. Data used in the approach, together with its source or the name of the product, the date and the resolution/density.

Data		Source or Product Name	Dates	Density/Resolution
Meteorological station data		AEMET (www.aemet.es) SMC (www.meteo.cat)	2003–2015	179–329 stations for each month
Other geographic data	Digital Elevation Model	Cartographic Institute of Catalonia (www.icgc.cat)	2004–2007	15 m
Satellite LST	MODIS	MOD11B3 and MYD11B3	2003–2015	5568 m
	AATSR and ATSR-2	ATSR LST Climate Data Record Level-3	2003–2011	0.05°

.

3.1. Meteorological Station Data

Month-to-month temperature records were obtained from the meteorological stations of the Spanish National Meteorological Agency (AEMET) [47] and the Catalan Meteorological Service (SMC) [48]. For long-term climatology, it is important to submit the data series to a rigorous quality control before mapping; this control has to be highly accurate to avoid unnecessary losses of spatial coverage while removing outliers that could negatively affect the resulting maps. In our case, spatio-temporal outlier detection was based on removing records outside three standard deviations from the mean values, combining the temporal perspective (one record was compared with all the records for the same station and for a specific month) with the spatial perspective (one record was compared with surrounding stations for a specific month and year). Finally, records that did not fulfill these two criteria were removed from the database.

In addition to removing outliers, we used expert knowledge (quality of the station series, position accuracy, etc.) and statistical indicators (robustness of the series based on the standard deviation of the number of monthly records) for the complete period of the series and for the studied period (2003–2015). Finally, to avoid artifacts in the interpolation step due to the proximity of some stations, a criterion of minimum distance (1 km) was applied. The criteria used to choose one of the stations closer than 1 km was based on the series that was longest and most robust. Once stations and/or records had been filtered, we selected the meteorological stations with temperature series longer than 15 years (considering the complete meteorological series for each station) or at least longer than 10 years in the studied period. This criterion was relaxed for recent stations with less than 10 years of records.

After applying the filter criterion, an average of 224 stations for each month were considered, varying between a minimum of 179 and a maximum of 329. Figure 1 shows the spatial distribution of these meteorological stations over the study area for the month with the maximum number of stations, as well as the 80 meteorological stations common for all cases. The minimum and maximum T_air values were observed around 5:00–6:00 a.m. and 3:00–4:00 p.m. local time respectively, although there was some variability depending on the location and time of the year. The accuracy of the T_air measurements was 0.1 K.

3.2. Other Geographic Data

The geographic variables were chosen taking into account their influence on the climate. We used variables whose influence is well documented in the bibliography (altitude and latitude). In addition, continentality, potential solar radiation and a topographic wetness index were also included to model the complexity of the climate in mountain areas. All variables were obtained from a Digital Elevation Model (DEM) and from different GIS analysis procedures.

Altitude (Alt.) is the nominal altitude above sea level of the stations, and it gives information on the variability caused by the relief. It was extracted from a DEM of the Cartographic Institute of Catalonia with 15 m of spatial resolution (with the ASTER Global DEM at 30 m resolution densified at 15 m to fill the contour zone). Finally, the whole set was densified to 90 m, the final resolution of the climatic maps to be obtained. Latitude (Lat.) is the Euclidean distance taking the 0° parallel (Equator) as a reference. Continentality was modeled by using the Euclidean distance to the Mediterranean coastline (Cont.), as well as the cost distance to the Mediterranean (CostMed) or the Cantabric Sea (CostCant) to account for the influence of mountain barriers. Monthly potential solar radiation (Rad) was extracted from a model proposed by [49] and extensively tested in [50], using the DEM and Sun position coordinates. The Topographic Wetness Index (TWI) [51] is the result of a function between the slope and the top-up contribution area. Although it was designed for hydrological applications, this index can be used as an indicator of the accumulation of cold air masses [52].

3.3. Satellite Land Surface Temperature

Operational Land Surface Temperature (LST) satellite products from MODIS, ATSR-2 and AATSR were considered in this study. They were derived by means of the split-window (SW) method that uses two channels within the 10.5–12.5 μm atmospheric window [53,54]. The main advantage of this method is that it does not require an accurate description of the atmosphere (in the form of vertical profiles of temperature and water vapor, ideally obtained from local, concurrent, radiosonde data) nor radiative transfer calculations based on these profiles. Initially created to estimate sea surface temperatures, the split-window technique was extended to land surfaces by accounting for the effects of surface emissivity (e.g., [55,56,57]). LST algorithms explicitly include a dependence on the surface emissivity in the channels considered, or alternatively different coefficient sets are provided for each land cover type. In any case, the characteristics of the surface must be well known (via the emissivity, or via the class and amount of vegetation cover) to obtain the LST, which is the main drawback of LST retrieval algorithms.

3.3.1. MODIS Data

The MODIS LST is derived from two thermal infrared band channels, 31 (10.78–11.28 μm) and 32 (11.77–12.27 μm). The atmospheric effects are corrected using the split-window algorithm [56] considering that the signal difference in the two thermal infrared bands is caused by differential absorption of radiation in the atmosphere [58]. The algorithm also corrects for emissivity effects, assuming they are known in each 1 km² pixel, using prior knowledge of the land cover type classification [59,60] based on the MODIS Land Cover product (MOD12C1). The MODIS sensor is onboard the Terra and Aqua spacecrafts, which acquire data both during the day and night (see the detailed time of acquisition of the images of this study in

Table 2. Local time interval (and average local time) of acquisition of satellite images after considering all the images and pixels of the study area.

Product	Daytime	Nighttime
MOD11B3	10:12–12:00 (11:11)	21:20–23:00 (21:59)
MYD11B3	12:24–13:36 (13:08)	01:24–02:36 (2:04)
ATCDR	10:11–10:59 (10:26)	21:08–22:01 (21:24)

).

The latest MODIS LST version, C6, is a higher quality product than previous versions [61]. This product aims to retrieve LST with an error lower than 1 °C (±0.7 °C standard deviation) in the range of −10 to 50 °C, assuming that surface emissivity is known [59,62]. Validation of the LST algorithm showed errors lower than 1 °C over homogeneous surfaces such as water, crop and grassland surfaces [62,63]. Detailed validation and major refinements implemented in the C6 MODIS LST product are given in [64].

In the present work, MODIS LST products onboard the Terra (MOD11B3) and Aqua (MYD11B3) spacecraft were considered [65]. Their nominal spatial resolution is approximately 5568 m and the products include the associated quality indicator layers. Each LST pixel value is a simple average of all the corresponding valid values in the MOD11B1 or MYD11B1 daily product collected during the month period. To include them in the regression models, images were adapted to 90 m spatial resolution by means of bilinear interpolation (or the nearest neighbor, in the case of quality layers). From now on, we will refer to products “MODd” or “MODn”, and “MYDd” or “MYDn”, where d and n refer to the daytime and nighttime products, and MOD and MYD to the Terra and Aqua products, respectively.

3.3.2. ATSR Data

The AATSR and ATSR-2 LST products are based on the split-window method [66] using the nadir view brightness temperatures at 11 μm and 12 μm. The algorithm expresses the LST as a linear combination of the brightness temperatures depending on the land cover type, the fractional vegetation cover and the precipitable water. It should be noted that land surface emissivity is implicitly handled within the retrieval coefficients dependent on the fractional vegetation. The regression coefficients can be implemented as a look-up table [66] and be updated in a routine manner; the method considers 13 land cover classes or biomes, covering from desert conditions to pine forests.

In this study, we used the ATSR LST Climate Data Record (CDR) Level-3 data, made available through the GlobTemperature Data Portal [67] and generated through the ESA DUE GlobTemperature Project with the support of the UK National Centre for Earth Observation (NCEO) [68]. These data are from the ATSR-2 onboard ERS-2 (from June 1995 to July 2002 inclusive) or from AATSR onboard Envisat (from August 2002 to March 2012), with a spatial resolution of about 0.05°. This prioritizes ATSR-2 data during the commissioning period of AATSR. Moreover, only full months of data were processed.

This product, hereafter called “ATCDRd” or “ATCDRn” for referring to the day or nighttime product respectively (see the detailed time of acquisition of the images of this study in Table 2), is constructed from the data in the Level-2 products (or fitted Level-2 data from the overlap analysis). It is made by mapping the Level-2 pixels onto equal-angle grids, and averaging the LST values—separately for day and night—of overlapping pixels in each grid using the percentage of overlapping area as weight. The input Level-2 products have been validated both within the framework of GlobTemperature and externally. For example, the AATSR LST datasets have been validated against several permanent stations of the Atmospheric Radiation Measurement (ARM) network. These show a bias lower than ±1.0 K with RMSE between 1.6 K and 2.2 K for daytime retrievals, and a bias lower than ±1.0 K with RMSE between 1.1 K and 1.6 K for nighttime retrievals depending on the surface [69]. The product includes the associated LST uncertainty for each grid cell, derived by adding in quadrature the components of the LST uncertainty budget (i.e., random, locally-correlated, systematic).

4. Methodology

4.1. Including LST in the Air Temperature Interpolation Scheme

Although there are different spatial interpolation techniques, we opted to apply a methodology based on multiple linear regression analysis combined with the spatial interpolation of residual errors using inverse distance weighting [12,32,70]. This interpolator offers the possibility of incorporating new geographic variables. In addition, this method has provided better results than other methodologies, at least in the case of T_air modeling [32]. The indicated methodology has been applied to produce the Digital Climatic Atlas of the Iberian Peninsula, which includes mean, minimum and maximum monthly and annual temperature layers; obtaining a RMSE of less than 1 °C in all months [33].

Multiple linear regression provides a predictive model of the climatic variable (T_air) using the variables that influence the climate of the zone (in our case, the LST has been added to the classic geographic predictors based on topography and distance surfaces). The result is a potential map obtained from the equation for fitting the regression that explains the general behavior of the climate. Once this potential mapping is available, it is possible to interpolate the residuals of the regression itself to obtain the potential surface closer to the observed data, and therefore, in general, for all the mapped territory. In other words, the use of unexplained variance from the regression model allows, once it is combined with the potential map, to obtain maps that are as realistic as possible. The residuals (difference between the observed value in the meteorological stations and the value predicted by the regression model) show local aspects of the climate not explained by the variables used in the regression or, in other words, they show how distinct the climate is at each point with respect to the general model.

Several reliability indexes were used to describe the thematic quality of each map, using a leave-one-out cross-validation procedure to compare the predicted and observed values from the meteorological stations.

In this study, to analyze the importance of using LST in T_air modeling, we used a combined model that includes both geographic and a remote sensing predictor (the LST) in the multiple regression analysis. Moreover, to quantify the improvement resulting from using LST, a multiple regression analysis was also carried out using only geographic variables. In addition, to analyze the implications of model complexity for the T_air estimate, two different sets of independent variables were employed to run the model: a first set that only included Alt., Rad., Lat., Cont. and LST as predictors (hereafter called the “simple” model); and a second set that also included TWI, CostMed and CostCant (hereafter called the “complex” model).

4.2. Segmentation Based on Land Cover and Orographic Complexity

To analyze the performance of the model in relation to certain land features that could be relevant, the study area was segmented based on land cover and orographic complexity (see Figure 1).

The study area was segmented into three main land cover categories: Forests, Agricultural and livestock, and Urban areas. This was carried out with version 4 of the Land Cover Map of Catalonia [71], which has 2 m spatial resolution and 25 categories (Level 2). Stations were then classified into categories by considering the mode within a radius of 500 m around each meteorological station. Depending on the month and year, a set of 33 to 38 stations were classified as Forests, 99 to 110 as Agricultural and livestock, and 44 to 59 as Urban areas.

Three levels of orographic complexity of the territory were considered based on the standard deviation of the altitude within a radius of 10 km around each meteorological station. Stations were then classified as Flat lands (deviation <150 m), Middle grounds (between 150 and 250 m) or Rugged lands (>250 m), which included sets of 120–131, 33–46, and 26–34 meteorological stations, respectively.

4.3. Analysis of LST Spatial Scale

The remote sensing LST products used in this study have a spatial resolution of around 6 km, while the regression models were applied to a spatial resolution of 90 m. The spatial scale of LST data was changed accordingly. This change in pixel size is required because the study presented here aims to analyze the improvements gained by including an additional variable in relation to the previous result. Therefore, although some precision was clearly lost by changing the scale of the LST data, it allowed us to achieve our objective. Notwithstanding, two different analyses were carried out to study the possible impact of the original coarse spatial resolution of the LST data on the T_air estimates, by segmenting the stations depending on: (1) the heterogeneity in solar radiation; and (2) the heterogeneity in land cover.

As the potential solar radiation did not follow a normal distribution within a ~6 km cell, its variability within a cell was quantified by considering the mean absolute deviation around the median (explained later in Section 4.5). A similar approach was carried out for land cover, but considering the percentage of the mode present within a ~6 km cell. It was then possible to create one class, named “homogeneous in solar radiation or land cover” (about a third of the stations), that had the lowest deviation values, and one class, named “heterogeneous in solar radiation or land cover” (another third), that had the highest values. Analysis of the performance metrics of each class and their comparison with the performance metrics of models that did not include LST provided some conclusions.

4.4. Weighted Linear Regression Based on LST Quality Bands

To use remote sensing products to derive information about the Earth’s surface it is necessary to know well, understand and, if possible, correct or minimize their limitations. In the present work, to eliminate low quality data, LST quality bands provided by the remote sensing products were used in the calibration step to weight the stations included in the model.

Quality information corresponding to each station was assigned by considering the mode (in the case of MODIS) or the mean (in the case of ATCDR, as the uncertainty is given as a continuous variable) within a radius of 10 km (to include at about 3 × 3 pixels) around each meteorological station from quality layers at the resampled spatial resolution (90 m). Then, the stations were weighted following the indications given in

Table 3. Weighting of the meteorological stations according to the quality masks provided by the MOD11B3 and MYD11B3 products. Land surface temperature (LST) quality control (QC) dataset, QC-LST, called as “QC_Mandatory_QA_flag”, where QA refers to quality assurance; Emissivity (Emis) quality control dataset, QC-Emis, called as “QC_Emis_error_flag”.

QC-LST	QC-Emis	Weighting (%)
1	1,2,3,4	100
3	1	70
3	2	50
3	3	30
3	4	10

and

Table 4. Weighting of the meteorological stations according to the uncertainty (δ) of the land surface temperature (LST) product from ATCDR.

ATCDR
δLST (K)	Weighting (%)
≤1	100
>1 to ≤1.2	90
>1.2 to ≤1.4	80
>1.4 to ≤1.6	70
>1.6 to ≤1.8	60
>1.8 to ≤2.0	50
>2.0 to ≤2.3	40
>2.3 to ≤2.6	30
>2.6 to ≤3.0	20
>3.0	10

for MODIS and ATCDR data, respectively.

Note that, as observed by [23], MODIS LST uncertainties increase with the MODIS viewing angle, with RMSE being greater than 2 or 3 K for viewing angles ≥40°. Thus, they only used MODIS LST measured with angles <40° in their study. Nevertheless, this filtering was not necessary in our study because only three images had viewing angles ≥40° and the corresponding stations were already weighted with a small value.

Benali et al. [16] also considered correcting low quality pixels due to non-detected cloud contamination and an error in the optical leak correction that occurred at the end of scan lines and led to anomalous high LST values. However, an accurate examination of our images did not reveal cloud contamination data or anomalous high LST values, probably because the monthly products considered here only averaged valid values.

4.5. Performance Metrics

Absolute root mean square error (RMSE_A) quantifies the scatter between measured (Meas) and estimated (Estim) values; thus, it can be used to make a quantitative assessment of the accuracy of our estimates:

$$RMS{E}_A=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(Esti{m}_i-Mea{s}_i\right)}^2}$$

(1)

where i refers to each of the N data points (meteorological stations or pixels) analyzed. On the other hand, the coefficient of determination (R²) provides a measure of how well measures are inferred by the model based on the proportion of total variation of measures explained by the model:

$$R^2=1p\frac{{\sum }_{i=1}^N{\left(Mea{s}_i-Esti{m}_i\right)}^2}{{\sum }_{i=1}^N{\left(Mea{s}_i-\overline{Mea{s}_i}\right)}^2}$$

(2)

When the distribution of variable x is not Gaussian, the median ($\widehat{x}$) instead of the mean is used as a measure of central tendency, and the mean absolute deviation around the median (σ′) is used as a measure of statistical dispersion. σ′ was calculated using:

$${\sigma }^{\prime }=\frac{{\sum }_{i=1}^N\left|x_i-\widehat{x}\right|}{N}$$

(3)

For the analysis based on the segmentation of meteorological stations according to their land cover and orographic complexity, or the heterogeneity in potential solar radiation and land cover, only the meteorological stations classified in each class were considered to calculate the metrics. However, all available meteorological stations were considered to analyze the overall performance, as well as to run the model.

5. Results

5.1. Effect of Considering the LST Quality

Performance metrics obtained from weighting linear regression based on LST quality bands were compared to those obtained when no weighting was applied. Detailed results from the analysis (not shown here for brevity), revealed that the impact on the R² and RMSE_A was minimal: R² remained the same and RMSE_A varied in general by <0.01 K, except for some months with variations up to 0.04 K as maximums when MODIS LST estimates were considered.

This analysis showed it was not necessary to add complexity to the model by weighting linear regressions, with the added bias downplaying relevant stations for the model due to the corresponding quality information of the LST estimates, which is merely indicative. Therefore, the following results did not include the LST quality data.

5.2. Effect of LST Spatial Scale

For this analysis, we used the case with the maximum number of meteorological stations (maximum T_air for March 2014). The statistical significance of the selected predictors and performance metrics are summarized in

Table 5. Selected model characteristics and performance metrics by considering the meteorological stations classified into two sets (homogeneous and heterogeneous) and regarding their surrounding potential solar radiation (Rad) or land cover attributes. The case analyzed corresponds to the maximum T_air observed during March 2014, including or not the LST provided by the daytime MYD11B3 product as an independent variable in the model. LST and Rad denote the standardized regression coefficients, × meaning non-significant (p > 0.05). N: number of meteorological stations considered.

	N	LST	Rad		RMSEA (K)		R2
			noLST	withLST	noLST	withLST	noLST	withLST
Homo. in Rad	109	0.28	×	×	0.80	0.77	0.55	0.61
Hetero. in Rad	109	0.08	0.069	0.074	0.78	0.78	0.97	0.97
Homo. in land cover	109	0.17	0.037	0.035	1.00	0.97	0.94	0.95
Hetero. in land cover	109	×	0.106	0.106	1.08	1.08	0.91	0.91
All meteo. stations	327	0.10	0.066	0.063	0.80	0.79	0.93	0.94

. A linear correlation analysis between the LST and the analyzed variables (Rad, σ′_Rad and the percentage of the mode in land cover) showed that there was no linear dependence between them. Despite this, these characteristics were chosen as they allowed us to analyze the effect of the LST spatial scale, as shown below.

5.2.1. Analysis in Terms of the Heterogeneity in Solar Radiation

In the analysis shown in Table 5, we observe that when meteorological stations classified as heterogeneous in potential solar radiation (Rad) within a cell of 6 km were considered, the Rad was included in the model as an independent variable. However, it was not the case when meteorological stations that were homogeneous in Rad were considered, which makes sense. Consequently, when models that included LST were considered, the weight of the LST in the regression was smaller in the case of heterogeneous than in homogeneous meteorological stations. A possible reason would be that, in the first case, part of the information was already provided by Rad, while in the second case, the LST provided much more new information. Another reason could be the loss of precision in the LST estimation for heterogeneous meteorological stations due to the change in the LST spatial scale. This effect was also observed in the improvement in the RMSE_A provided by including LST, which was greater for homogeneous (~0.03 K) than for heterogeneous (<0.01 K) stations. Nevertheless, what is important is that, regardless of the accuracy of the LST estimates, the inclusion of this variable was significant (as it was considered as an independent variable by the regression model) and the RMSE_A associated with the T_air estimates did not get worse.

The R² was much higher when meteorological stations that were heterogeneous in Rad were considered. This is understood to be due to the increase in the dispersion of values, thus providing a better adjustment.

5.2.2. Analysis in Terms of the Heterogeneity in Land Cover

Similarly, we observed that the heterogeneity in land cover downplayed the information provided by the LST to the point that it ceased to be an independent variable in the model (see Table 5). As was the case for Rad, we think that the heterogeneity in land cover unified the LST (it lost precision), and hence, this variable did not provide new information to the regression model. However, the opposite happened in the case of homogeneous meteorological stations. Therefore, it is obvious that an improvement in the spatial resolution of LST would improve the T_air estimates. Nevertheless, once again, we observe that, regardless of the accuracy of the LST estimates, the inclusion of this variable was significant (as it was considered as an independent variable by the regression model when all meteorological stations were considered) and the RMSE_A associated with the T_air estimates did not get worse (in fact, it decreased by 0.02 K).

5.3. Statistical Comparisons of Model Performance

5.3.1. Models with a Set of Common Meteorological Stations

To compare accuracy metrics more precisely, independently of the months and the number of data points considered, we selected the years with LST data from all sensors (2003–2011) and the 80 meteorological stations common for all cases (see their location in Figure 1a). This allowed us to study the impact of including LST in the model (Figure 2a) and to decide whether to use it in the research presented here. In addition, it permitted us to examine whether the daytime or nighttime LST product was the most appropriate for estimating maximum or minimum T_air (Figure 2b), as well as deciding which LST product usually provides the best T_air estimates (Figure 3).

From the analysis in Figure 2a, we observe that including LST in the models always implied an improvement (in terms of RMSE_A differences) in the T_air accuracy, especially for the minimum T_air (being up to 0.35 K in December and over 0.15 K in the other months). In the case of maximum T_air, the improvement was not as large (<0.12 K), except in June (up to 0.17 K).

In addition, Figure 2b revealed that, except for maximum T_air in June, both maximum and minimum T_air were clearly better estimated by the nighttime LST products. More specifically, as indicated by Figure 3, the nighttime ATCDR product provided the best proxy for minimum T_air for all seasons and for maximum T_air in March, while it was MYDn for minimum T_air in September and December. Quantitatively, the differences in the RMSE_A given by the other nighttime products were up to 0.10 K, with a mean value of 0.03 K. As an exception, in June the maximum T_air was better estimated by the daytime product from Aqua (MYDd), showing RMSE_A differences with the other products of up to 0.11 K.

5.3.2. Analysis Segmented by Land Cover and Orographic Complexity

According to Figure 4, maximum T_air estimates were more accurate than minimum T_air independently of the land cover or orographic complexity. However, there was no clear pattern between months, while the highest errors for maximum and minimum T_air occurred for June and December respectively. Regarding the land cover, Figure 4 shows that meteorological stations classified as Agricultural and livestock had the best accuracy in all seasons, with mean RMSE_A values of 0.82 K and 0.98 K for maximum and minimum T_air respectively, and maximum RMSE_A values <1 K. Forests and Urban areas showed similar accuracies for minimum T_air, while for maximum T_air Forests were less accurate than Urban areas. In relation to the orographic complexity, there was no clear dependence for the T_air accuracy, although we could say that Rugged lands showed the lowest and highest accuracies for maximum and minimum T_air, respectively.

Figure 5 depicts how the improvement in T_air due to including LST was clearly higher for minimum than maximum T_air independently of the land cover or orographic complexity. Meteorological stations classified as Forests showed the most improvement, with maximum differences between 0.18 and 0.29 K for the maximum T_air (in comparison to differences <0.11 K for the other land covers), and between 0.29 K and 0.43 K for the minimum T_air (in comparison to differences <0.13 K for Agricultural and livestock and <0.28 K for Urban areas). Meteorological stations classified as Rugged lands had the biggest and clearest improvement for maximum T_air, but the smallest improvement (except for December) for minimum T_air, for which Flat lands showed the best improvement.

5.3.3. Implication of Model Complexity

As shown by the statistics in Figure 6a,b, the accuracy of T_air improved (ΔRMSE_A > 0 K) when all considered independent variables were included, while it was more significant for maximum T_air (with differences up to 0.17 K, and mean values from 0.06 K to 0.09 K, except in December). However, the improvement for minimum T_air was slight (with differences <0.05 K) or even inexistent in some cases. The R² decreased with the complexity of the model by up to 0.11 K for maximum T_air (Figure 6c), although there was less of a decrease (0.03 K as maximum) for minimum T_air (Figure 6d).

5.3.4. Overall Performance

To give an idea of the overall performance achieved with our method, RMSE_A and R² value metrics corresponding to all years and months are summarized in Figure 7. To simplify, for these plots the model considered the selected LST products indicated in Figure 3. As already seen in Section 5.3.2, maximum T_air had better accuracies than minimum T_air, even though in neither case did they exceed 1.6 K. Mean RMSE_A values ranged from 0.8 to 1.1 K for maximum T_air, and from 1.1 to 1.2 K for minimum T_air. Similarly, the goodness of fit (explained by the R²), was much better for maximum than minimum T_air, especially for some of the warmest months (March, June and September). Mean values varied from 0.90 to 0.95 and from 0.83 to 0.88 for maximum and minimum T_air, respectively.

In addition, to test the consistency of the model, we analyzed the linear relationship between the calculated and ground measured T_air. As an example, Figure 8 shows the dispersion of values observed for maximum and minimum T_air, together with the linear fit and the statistics, which were slightly better for maximum T_air. We also carried out a more in-depth analysis that considered the months separately, and segmented by land cover or by orographic complexity. The analysis results (not shown here for brevity) revealed that the performance metrics were practically equivalent to the mean values shown in Figure 4 and Figure 7.

Figure 9 gives evidence of the improvement in accuracy for T_air estimates gained by including LST. The improvement was greater for minimum T_air, except in June. Differences in mean RMSE_A for maximum T_air ranged from 0.015 to 0.06 K, and from 0.05 to 0.12 K for minimum T_air. Nevertheless, differences were up to 0.1 K and 0.2 K for maximum and minimum T_air, respectively. From the comparison with Figure 6, we observe that including LST as a predictor of minimum T_air led to a much larger improvement than that provided by additional geographic data (i.e., TWI, CostMed and CostCant), while its contribution to maximum T_air was less important. The regression coefficients of the model (Figure 10) also give evidence of this effect, as for minimum T_air the LST variable is the most important one, and for maximum T_air its weight is moderate, being the altitude the most relevant parameter. The analysis of results observed by running the model considering only the LST as independent variable (Figure 8c,d) gives again evidence of the improvement provided by such variable. Specifically, LST by itself can provide a high accuracy in minimum T_air (1.41 °C) in comparison with results from a model considering all independent variables (1.16 °C), while the achieved accuracy is not so high for minimum T_air (2.16 °C vs. 1.01 °C).

While all these results only concerned values corresponding to the meteorological stations, Figure 11 depicts the spatial and temporal variation over all the study area during a selected year. Although these maps do not show which model is the best (as there is no comparison with ground measurements), they provide a spatial and temporal quantitative view of the differences between the results of the model that considers LST and the model that does not. These differences were smaller in March and much more pronounced in the other months, especially over high altitude areas in the northwest of Catalonia, which corresponded to Rugged lands (see Figure 1b), with absolute differences exceeding 1.2 K. More details are observed in Figure 12, which shows the spatial difference over areas with a dominant land cover. The highest T_air differences were observed over Urban areas (black and blue pixels in Figure 12c), followed by some Forest areas (green pixels in Figure 12c and blue pixels in Figure 12a) and Agricultural and livestock regions (cyan pixels in Figure 12c and the different shades of green in Figure 12b).

6. Discussion

6.1. Improvements from Including LST

This study shows the strong physical relationship between LST and T_air and evidence the potential of using MODIS and ATCDR LST monthly data to estimate month-to-month T_air dynamics at a detailed spatial resolution (90 m), despite their rough original spatial resolution (~6 km).

According to [16], daytime and nighttime LST are the most relevant predictors. From comparisons of Figure 6 and Figure 9, we observe that our study is clearly in agreement with these authors for minimum T_air, while it is not so evident for maximum T_air. Nevertheless, whatever the range of improvement achieved, it is clear that including remotely sensed LST estimates as a predictor in classical T_air modeling is positive (Figure 2 and Figure 5).

6.2. Differences between Daytime and Nighttime LST

As observed by different authors [16,17,18,19,20], besides data availability, cloud contamination and angular effects, the time of overpass of the remote sensor has an impact on the accuracy of T_air estimates. However, [21] stated that the difference between satellite overpass and the timing of maximum or minimum T_air was not critical.

Regarding our considered LST products, MODIS Aqua time of overpass is closer than MODIS Terra and ATCDR to the time of minimum and maximum T_air (Table 2). Similarly to [18,19] and contrary to [22], we found that the RMSE_A of T_air estimation was higher using MODIS Terra than Aqua, although it should be noted that ATCDR provided the best estimates (due to differences in the accuracy of the products).

Furthermore, it is important to note that, generally, most of the studies [16,18,23] consider nighttime LST data as a covariate to predict minimum T_air, and daytime LST data to predict maximum T_air. Our research, in agreement with findings by [20,37] over complex terrain and by [19] over the Corn Belt in the U.S., showed that nighttime LST provided a good proxy of minimum T_air, while models with daytime LST had the lowest accuracy. However, in agreement with findings by [17] and [19] and following recommendations by [20], our study also showed that nighttime LST provided a good proxy of maximum T_air. Interestingly, by combining both daytime and nighttime LST data, it is possible to improve the accuracy of maximum T_air, while no significant improvement in accuracy is found for minimum T_air estimates [19], indicating that daytime LST is not relevant for estimating minimum T_air. In addition, as observed by [16], combining MODIS Terra and Aqua data also has great potential for T_air estimation. Therefore, future research could focus on combining data sources, as well as daytime and nighttime LST estimates.

6.3. Temporal and Spatial Patterns

The derived datasets were representative of T_air dynamics for the whole study region, but also for regions characterized by a predominant land cover and by orographic complexity. Our analysis showed that model performance did not follow a clear seasonal pattern for either maximum or minimum T_air, while there were significant differences between accuracies for maximum T_air but not for minimum T_air. Benali et al. [16] and Zeng et al. [19] observed higher accuracy for daily T_air estimates in summer, which they attributed to the prevalence of clear skies. On the contrary, we observed the lowest accuracy for maximum month-to-month T_air estimates in summer (June), followed by those in winter (December), and the highest accuracy in spring (March).

From an analysis at station level, Benali et al. [16] and Zeng et al. [19] found that, in general, the meteorological stations located in forests had a better performance than those in crop and urban areas, with the lowest performance for stations in urban areas. From the analysis of Figure 11 and Figure 12, we observed that our performances followed the same pattern. Nevertheless, it is worth mentioning that the statistics in Figure 5 are probably not strictly comparable between them, except for the seasonal variation within each class, due to the different number of stations classified in each class.

Similarly, the analysis results in Figure 11 and Figure 12 showed that the month-to-month T_air accuracy improved with the decrease in terrain roughness, in accordance with the findings by [20,19]. Further, in agreement with [18], we also observed high errors in minimum T_air over water surfaces due to the inertia of water, which may lead to an underestimation or an overestimation of minimum T_air according to the temperature of the previous months, respectively cooler or warmer than the present [18].

7. Conclusions

The good results of this study can be considered operational given the feasibility of the models employed (linear dependence on predictors that are already easily available), the robustness of the leave-one-out cross-validation procedure and the improvement in accuracy achieved compared to classical T_air modeling results.

The results were within or above the accuracies reported in the literature considering the methodologies for T_air estimation based on remote sensing LST. The model had a general tendency to present lower performance but higher improvements for minimum T_air in comparison with maximum T_air. This is especially relevant as month-to-month minimum T_air prediction using only GIS-based geographical variables often provides unexpected results. Quantitatively, on the one hand, the models estimated maximum and minimum month-to-month T_air with mean RMSE_A within 0.8–1.1 K and 1.1–1.2 K, respectively, and mean R² within 0.90–0.95 and 0.83–0.88, respectively. On the other hand, the inclusion of LST into the models implied an accuracy improvement up to 0.35 K in relation to classical T_air modeling on December and beyond 0.15 K on the other months for minimum T_air, while for maximum T_air it was not as important (<0.12 K), except on June (up to 0.17 K). The improvement was higher over Forests and Rugged lands. Results are really encouraging as there are generally few meteorological stations in zones with these characteristics, clearly showing the usefulness of remote sensing to improve information for areas that are difficult to access or simply with a poor availability of conventional meteorological data. Such improvement was thanks to the LST, which provides additional information in relation to phenomena as the thermal inversion, not reflected so far, or not enough, by geographic variables or in situ measurements, probably due to the few meteorological stations located on rugged lands.

Although it might seem that the improvement was not very high, it is important to note that the improvement provided by the new model was clear and robust and that can be applied routinely. Currently, within the study of trends in the context of climate change, such improvement is relevant, as spatial and temporal finer resolutions are often required and, of course, higher accuracies are needed. Moreover, this is in line with the increase of operational LST products developed by some Spatial Agencies. A wide scope of applications can benefit from the improved remote sensing based T_air estimations presented in this work. Further, as an added value, we would like to mention that this work will imply the update and improvement of the “Digital Climate Atlas of Catalonia” developed by the Autonomous University of Barcelona and available online [72].