Spatial Distribution of Surface Temperature and Land Cover: A Study Concerning Sardinia, Italy

Land surface temperature (LST) is a key climate variable that has been studied mainly at the urban scale and in the context of urban heat islands. By analyzing the connection between LST and land cover, this study shows the potential of LST to analyze the relation between urbanization and heating phenomena at the regional level. Land cover data, drawn from Copernicus, and LST, retrieved from Landsat 8 satellite images, are analyzed through a methodology that couples GIS and regression analysis. By looking at the Italian island of Sardinia as a case study, this research shows that urbanization and the spatial dynamics of heating phenomena are closely connected, and that intensively farmed areas behave quite similarly to urban areas, whereas forests are the most effective land covers in mitigating LST, followed by areas covered with Mediterranean shrubs. This leads to key policy recommendations that decision-makers could implement to mitigate LST at the regional scale and that can, in principle, be exported to regions with similar climate and land covers. The significance of this study can be summed up in its novel approach to analyzing the relationship between LST and land covers that uses freely available spatial data and, therefore, can easily be replicated in other regional contexts to derive appropriate policy recommendations.


Introduction
Fast changes aimed at promoting social and economic development have characterized international land cover dynamics in the last few decades [1]. Land cover features are related not only to rapid urbanization, but also to processes of higher anthropization, such as the transformation of forests into agricultural areas or transitions from woodlands and shrubs to new predominantly homogeneous agrarian land [2]. According to the National Institute for Environmental Protection [3], artificial land cover in Italy increased by 0.21% in 2018, that is by 23,033 km 2 .
Land cover changes, and especially transitions from natural and semi-natural areas to artificial land covers, affect regional and local temperatures [4]. In fact, although urban areas and their surroundings receive the same amount of solar radiation, the local temperatures differ because different surface materials have different heat capacity [5]. From this perspective, land surface temperature (LST) is a significant parameter for investigating the effects of land covers on local temperatures.
Hofierka et al. [6] define LST as "the radiative skin temperature of the ground", affected by solar reflectance, thermal emissivity, and heat capacity. In other words, LST combines interactions between land surface and atmosphere with ground-atmosphere energy fluxes. Therefore, understanding how land cover changes affect climate represents a key element in international debates [7] that highlight that land cover can significantly affect quality of life, that is human health and safety, through its influence on LST [8]. The impacts on climate conditions generated by land-cover change processes can be effectively analyzed by assessing the relations between the spatial distribution of LST and that of land covers. The influence of land use/land cover changes on LST variation has been studied by various authors [4,[9][10][11][12]. Feizizadeh et al. [4] analyze the relations between LST and land use/land cover in Maraqeh County (Iran), using a method based on the application of the Surface Energy Balance Algorithm for Land (SEBAL) to Landsat Enhanced Thematic Mapper (ETM+) imagery. Chaudhuri and Mishra [9] compare land use/land cover and LST differences between India and Bangladesh by analyzing satellite images corresponding to different time periods and by using intensity analysis. Zullo et al. [10] study the relationship between LST variations and the increase in urbanized areas from 2001 to 2011 in the Po Valley. Guha et al. [11] investigate the relationship between LST and two indexes, that is the normalized difference vegetation index (NDVI) and the normalized difference built-up index (NDBI) in Florence and Naples by using Landsat 8 data. Stroppiana et al. [12] focus on the relationship between spatio-temporal LST variation and land cover, topography, and potential solar radiation in four southern Italian regions.
Although various authors have studied the relations between land use/land cover and LST, further analysis may provide more accurate parameters in order to understand this complex relationship [13] at a wider scale. Moreover, a mere understanding of this relation is not sufficient to deal with climate change, while new policies and strategies should be included within spatial planning at the regional and local levels.
In this study, the connections between anthropization and natural conditions and the spatial taxonomy of LST are analyzed in order to assess whether, and to what extent, land covers and their transitions affect the spatial layout of heating phenomena at the regional level, in the case study of Sardinia, Italy. The methodological approach applied in this study, which is based on a cross-sectional spatial inferential analysis, provides results that can be used in two ways. On the one hand, the differential impacts of different land cover types on LST are characterized with reference to the LEAC classification, and the most critical land covers related to heating phenomena are identified. On the other hand, a number of planning policies concerning land cover changes are defined, with the ultimate aim of decreasing LST.
The study is structured into five sections. The first section identifies the wider context and the debates to which the study is contributing. The second section describes the study area, data, and methodology used in this study. The third section presents the results that are discussed in the fourth section, where the main findings are interpreted and compared with findings from similar studies. The fifth section provides recommendations in terms of strategies and policies, and proposes concluding remarks and future directions of the research.

Study Area
With a size of 24,000 km 2 and approximately a 1,850-km long coastline, Sardinia is one of the major Mediterranean islands. The island, which from an administrative point of view is an Italian autonomous region populated by around 1.6 million inhabitants, is taken as the spatial context for this study because its climate homogeneity and self-containment allow for a pretty straightforward identification of the regional boundaries.

Data
Landsat 8 OLI (Operational Land Imager) and TIRS (Thermal Infrared Sensor) images are freely available from the US Geological Survey website [14]. For spatial criteria, a bounding box having minimum latitude 38.8, maximum latitude 41.4, minimum longitude 8, maximum longitude 10 was selected; for temporal criteria, the spring interval spanning from 15 April 2019 to 31 May 2019 was chosen. The five selected multi-band images are listed in Table 1. The cell size of the thermal bands (10 and 11), for each image, was 30 meters. Finally, a Digital Terrain Model (DTM) is available for Sardinia from the regional geoportal [17] and its cell size equals 10 meters.

Methodology
The methodological approach of this study develops as follows. First, LST was extracted and mapped. Secondly, land covers and elevation data were processed so as to build a spatial dataset comprising the three layers; such dataset was next used to feed a multiple regression analysis to estimate the relationships between different land covers and LST. A graphical summary of the methodology adopted in this study is provided in Fig. 1. LST Extraction and Mapping at the Regional Scale A QGIS plugin implemented by Ndossi and Avdan [18] was used in this study. For each image listed in Table 1, a five-step process was therefore performed through the plugin. In the first step, the top-of-atmosphere spectral radiance was calculated for each pixel [19]; subsequently, the top-of-atmosphere spectral radiance was converted into the top-of-atmosphere brightness temperature [19]. Next, the NDVI was calculated using Landsat 8's bands 4 and 5 images as inputs [20]; NDVI makes it possible to calculate, in the fourth step, Land Surface Emissivity (LSE) through various algorithms. Among those implemented in the plugin, Zhang, Wang, and Li's algorithm [21], which builds upon van de Griend and Owe's findings concerning the correlation between LSE and NDVI [22], was chosen here because it was reported to yield the best results [18] for LST retrieval. Finally, in the fifth step, LST was calculated using the so-called "Planck function" [23], which has been reported to be "easier to use in comparison to the other algorithms as it does not require atmospheric variables" [18] (p. 28). Through these steps, five LST raster maps (one for each Landsat image listed in Table 1 and each having resolution 30 meters) were obtained. The five LST images were then merged, so as to obtain a single LST image for the whole region (in the case of overlapping pixels, the maximum LST value, consistently corresponding to pixels belonging to scene 193 produced on May 23 was retained), which was next resampled to a 300-m cell size so as to lower the resolution, hence reducing computational efforts in the subsequent steps.

Land Covers and Elevation
From the European 2018 CLC vector dataset, polygons concerning Sardinia were first extracted and next reclassified based on the LEAC classification defined by the EEA [24] (p. 98) and used for land cover accounts. Table 2 shows the LEAC-related macroclasses as groups of the 44 third-level classes of the CLC dataset. The asterisk (*) marks any sub-classes of a given class, or any sub-sub-classes of a given sub-class.
Next, the resulting layer was converted into a raster map with a cell size of 300 meters (where the grid was the same as that of the resampled LST). Each cell was assigned a value corresponding to the prevailing land cover group.
As for elevation, the regional DTM was resampled so as to obtain a new raster file with the same cell size and grid as the LST and land cover group raster images. As for the resampling technique, linear interpolation, appropriate for continuous data [25] as elevation, was used to derive the cell's altitude.

Input Table for the Regression
A vector layer (shapefile) was created. Polygons in this shapefile spatially coincide with pixels belonging to the three raster maps derived as per the previous subsubsections; hence, every polygon is a 300 m by 300 m square in the projected reference system EPSG 32632 [26]. Each polygon was then assigned the following attributes: the LST value of the corresponding cell, its prevailing land cover group, and its altitude. Furthermore, the latitude of the polygon's centroid (in the reference system EPSG 32632) was added, since the temperature is expected to be dependent on the latitude. Finally, a field taking the value 1 (pixels pertaining to scene 193) or 0 (pixels only pertaining to scene 192) was added. The pixels featuring in both scenes were regarded as belonging to scene 193.

Land Covers and LST
The polygons described in the previous subsubsection were identified as the spatial units for estimating a multiple linear regression which takes the following form: where: • explanatory variables from URB through to OPSP, which represent the LEAC groups, are dichotomous variables; each variable can take either of the two values, 1 or 0, according to whether the area size of a LEAC group in a cell takes the largest value with respect to the area sizes of the other groups; therefore, if in a cell the URB group shows the largest area size, the variable URB equals 1, otherwise it equals 0, and so on; each coefficient estimated by regression (1), β i , i = 1, …, 6, identifies the change in LST related to a cell in case it shows the largest area size identified by the variable associated to the coefficient β i (i.e., URB, APC, etc.) with respect to the basic condition that the largest area size of the cell is identified by the variable WAT (Wetlands and water bodies); the coefficients estimated by regression (1), β i , i = 1, …, 6, define a taxonomy of the zone types based on the quantitative contribution to LST expressed by the values of β i , i = 1, …, 6; • ALTIT is the altitude in meters related to the polygon; • LATD is the latitude in meters of the polygon's centroid; • W (standing for "West") is a dichotomous variable which can take either of the two values, 1 or 0.
The results from the multiple linear regression were used to identify the quantitative effects of different land covers on LST and their ordinal ranking; for each LEAC group, the rank depended on the value of the coefficients βi, i = 1, …, 6, of regression (1).
That being so, the surface, whose equation was unknown, representing an ndimensional phenomenon, e.g., a spatial phenomenon, as a functional relation between n variables, was approximated, point by point, by the point neighborhood identified by the tangential hyperplane. The neighborhood shared by the surface, whose equation was unknown, and by the tangential hyperplane, whose equation was known, was identified by the linear relation between the covariates, as a locally identified approximation of the unknown general relation between the factors. A multiple linear regression such as (1) provided the estimates of the coefficients of the linear equation, which demonstrated the trace of a tangential hyperplane over a nine-dimensional unknown surface which represented the functional relation between LST and its covariates [27].
The variables ALTIT and LATD were used to control the impacts of the altitude and latitude of a cell on the LST; as a consequence, if the estimates of the coefficient β 7 and β 8 were significant with reference to their p-values, this would entail that the altitude and latitude generate effects on the LST, whose expected signs are negative, since it is expected that the higher the altitude, or the greater the latitude, the lower the LST, everything else being equal. Finally, the variable W was used to control the impact of the date of the satellite images, which should show a systematic difference between May 23, 2019 and May 16, 2019, since May 23 was a clear sunny day whereas May 16 was partly cloudy. That being so, the cells whose LSTs were extracted on May 23, belonging to scene 193 and hence located on the western side of the island, should show higher LSTs than the others, everything being equal, and the expected sign of W should be positive.
A standard significance test based on p-values was implemented as regards the coefficients of the regressions, to detect whether they are significantly different from zero. If p-values were lower than 5%, then the explanatory variables' effects on LST could confidently be considered important.

Results
The outcomes of the methodological approach discussed in the second section are presented as follows. In the first subsection, the spatial layout of the LST taxonomy is described; in the next subsection, the estimates of regression model (1) are reported and the effects of different LEAC macroclasses on LST are presented, which define an ordered scale whose lowest level is represented by water bodies and wetlands and upper level by urbanized land.

LST Spatial Taxonomy Results
Some outputs of the LST extraction process are shown in Fig. 2, which clearly shows that some parts of the island were covered by clouds when the Landsat images where produced-this locally affected both NDVI and LST values. Hence, such pixels (25,099 out of 266,818) were subsequently removed; Fig. 3 shows the spatial distribution of LST minus the cloudy cells, of LEAC groups, and of elevation that were used to develop the matrix that fed the regression.

Regression Results
The estimates of the regression identify the impacts of the land cover macroclasses of the LEAC on the LST. The estimates of the coefficients of the Boolean covariates define  the differential effects of each variable on LST change with respect to "Water bodies and wetlands", whose effect on LST is the lowest among the LEAC macroclasses. Table 3 shows the results of the regression model.
The estimates of the coefficients of the altitude and latitude variables are significant in terms of p-value hypothesis testing, and they show the expected sign. As a consequence, it can be stated that the higher the altitude, or the larger the latitude, the lower the LST. On average, an increase of 100 m in altitude will imply a decrease of 0.64 K in temperature, whereas an increase of 10 km in latitude will entail a decrease of 0.16 K in temperature. Moreover, the coefficient of the dichotomous variable W is significant and it shows the expected sign as well, which implies that, on average, a cell whose LST was related to the satellite images taken on May 23 is 4.3 K higher than the LST of a cell whose LST was related to the satellite images taken on May 16. Therefore, the estimates of the three control variables' coefficients are significant and show the expected sign, and, that being so, it is pretty straightforward to identify the implications concerning the differential effects on LST derived from the six LEAC macroclasses on the basis of the estimated coefficients of the six dichotomous variables from URB through OPSP.
The estimated coefficients of the six dichotomous variables are significant with respect to the standard test based on p-values described at the end of Sect. 2. The estimates of the coefficients entail the following implications, which all assume the condition "everything else being equal". First, the highest effect on LST is shown by URB. The estimated coefficient implies that, on average, urbanized (artificialized) cells present an LST higher by: 8.9 K than WAT cells (water bodies and wetlands); 0.4 K than APC cells (arable and permanent crops); 1.5 K than MCP cells (mosaic crops and pastures) and OPSP (open spaces with sparse or absent vegetation); 2.6 K than HNGS cells (heathland, natural grasslands, and sclerophyllous vegetation); 4.1 K than FSW (forests, shrubs, and woodlands). Secondly, intensive farming (APC) generates impacts on LST similar to urbanization, even though it is not characterized by the soil-sealing phenomenon. Thirdly, extensive farming and pastures (MCP), and bare land or poorly vegetated soils (OPSP) are correlated to lower LST, and their impacts on LST are similar to each other. Fourthly, the two LEAC non-agricultural vegetated macroclasses are the most effective as regards mitigation of the land surface heating phenomenon. Forests, shrubs, and woodlands (FSW) show the largest positive impact on LST mitigation.

Discussion
The ordered scale of the effects of the LEAC macroclasses on LST derived from the estimates of the regression model shows a number of aspects consistent with the findings reported in the current literature. The most relevant impact on increases in LST is represented by urbanized land, whose differential effect on LST with respect to water bodies and wetlands is about 9 K. Several reasons can explain this outcome. Sealed soils either limit or prevent air circulation and the impact of downwind cooling [28]. In heavily urbanized areas, the thermal comfort of vegetated areas is almost completely missing [29]. Artificial surfaces make evapotranspiration almost impossible, and, as a result, LST and air temperature are comparatively higher than in non-artificial surfaces [30]. In addition, the surface materials used in urban areas have higher radiant temperatures [13]. The increasingly frequent heat-wave and heat-island phenomena would be mitigated by the canopies of woodland trees and shrubs, which are rare or absent in urbanized areas [31].
In relation to intensive farming, the outcomes of the regression model show that the impact on LST increase generated by intensive farming areas is only slightly lower than that of urbanized land, that is, about 8.5 K higher than that of the water bodies and wetlands LEAC macroclass. Indeed, the areas belonging to the APC macroclass reveal conditions that are very similar to those of the URB macroclass. Even though the soil is not sealed, arable land and permanent crops generate negative effects on downwind cooling and air circulation which make LST comparatively higher than those of the other LEAC macroclasses except artificial land. Even the thermal comfort and evapotranspiration provided by areas characterized by intensive agriculture are very low with respect to that of the other LEAC macroclasses, with the exception of the URB macroclass, since in artificial areas trees are almost totally absent and soils are generally covered by dense low-growing vegetation [32].
HNGS shows lower LST compared to other macroclasses due to the presence of vegetation that reduces the amount of stored heat in the soil through transpiration [33]. In fact, vegetation affects the microclimate through shading that prevents incoming solar radiation reaching the land surface through evapotranspiration and by affecting air movement and heat exchange [29].
The results of this study reveal that forests, shrubs, and woodlands show the larger positive impact on LST mitigation. According to a study carried out by Walawender et al. [34] in relation to Krakow, forests and waters show low values of LST. However, their behavior changes during the year. In March, before the period of growth between germination and flowering, waters show considerably lower values of LST compared to forests. In May, when the vegetative period is just beginning in Poland, the LST values of forests are slightly lower than the values of waters. In June and August, at the end of the vegetative period, forests and waters take the same values. This different behavior is caused by evapotranspiration that is superior at the peak of the vegetative period when air temperature goes up. The same study shows that the LST of arable land, pastures, and permanent crops is affected by seasons. At the beginning of the vegetative period, the LST of arable land is slightly higher than the mean values. Before the vegetative period, the LST values of permanent crops and grassland are higher than the mean values. This temperature variation can be explained through the state in which arable land is kept after the planting period, that is not covered by any vegetation. In April-May, during the early vegetative season, arable land shows the same thermal properties as bare soil, and its LST values are higher than areas covered with vegetation.
On the other hand, pastures are perennial; therefore, the early vegetative period is characterized by a closed canopy different from the bare soil that characterizes arable land in the same vegetative period. Moreover, the capacity of regenerating living biomass is strongly connected with precipitations that limit a constant forage production during the vegetation period [35].
It must be noted that lowering LST in rural areas has relevant implications in terms of mitigation of water shortages and of the entailed negative economic impacts and social problems [36]. As per Sruthi and Aslam [37], LST and NDVI are strictly negatively correlated in agricultural areas, which implies that the higher the LST, the lower the vegetation density. This outcome is particularly important with regards to periods of drought, characterized by productivity decline determined by poor rains and soil humidity [37]. In turn, the decline in productivity can possibly generate significant consequences in terms of economic and social unrest [38]. This is an outstanding issue with reference to Sardinia, where agricultural area comprises about one half of the regional land (11,500 km 2 ) and employment in agriculture is about 7.5% of total employment (41,000 people) [39].

Conclusions
A number of policy recommendations are implied by the outcomes of this study.
First, high LSTs in the consolidated fabric of cities, towns, and small villages, which may possibly develop into heat islands and waves, can be mitigated through the implementation of widespread urban microscale operations, such as the planting of trees, the increase in the endowment of urban green zones through the plantation of new areas and/or the enlargement of existing ones [29]. Other microscale measures consist of the realization of green facades and walls, and green and blue urban grids. There are several examples of these planning policies implemented in different urban contexts. A pioneer and outstanding primer is the London Green Grid [40].
The microscale planning measures aimed at decreasing LST in urbanized contexts are based on the provision of ecosystem services regulating LST. The implementation of these measures implies the integration of different kinds of planning policies, in order to induce virtuous behaviors on behalf of the local communities, organized groups of residents, building enterprises, and public administrations. A central issue is that the value of urban land is strictly related to its permitted building volume, whether it is residential or devoted to services or infrastructure. That being so, the plantation of new green areas and/or the enlargement of existing ones, which imply that these areas would decrease their value substantially by losing their building potential, should be implemented through the integration of appropriate planning measures. On the one hand, it is necessary to enact a regulatory system whose rules should state that new settlements, and, perhaps, existing ones, have to be endowed with a certain amount of vegetated areas, either by devoting part of the facades and roofs to greenery, or by binding more-or-less large areas to be part of blue or green grids, such as in the case of East London [41]. On the other hand, a system of financial incentives should be put in place so as to make buildings endowed with green facades and roofs, as well as the realization of green and blue grids in the existing urban neighborhoods and in the new developments, attractive in terms of private investment [42]. Finally, public commitment to the increase in the urban endowment of blue and green grids and infrastructure should be identified clearly through the implementation of public planning policies, for example through compulsory purchase orders concerning municipal areas where urban greening plans, such as the East London Green Grid, are to be implemented [43].
Secondly, as regards non-artificial land, the most effective LEAC macroclasses in lowering LST are forests, shrubs, and woodlands (FSW) and heathland, grasslands, and sclerophyllous vegetation (HNGS),whereas the impacts of intensive and extensive agriculture-related land covers (APC and MCP), and of open spaces (OPSP) are very close to the negative impact of urbanized land. A pretty straightforward implication is that the LEAC macroclasses characterized by the lowest degree of anthropization, such as FSW and HNGS, should be targeted by public policies in order to implement the mitigation of the land surface heating phenomenon. Policy measures should aim at supporting stepwise transitions from APC, MCP, and OPSP to HNGS and FSW.
Afforestation is the main road to identify the implications of the regression model results as regards lowering LST in rural areas [44]. Moreover, for afforestation, a distinction should be made between high-rent intensive arable land (ARA) and low-rent mosaic farmland and pastures (MCP) [45]. Transitions from agriculture to forests are not likely to take place in cases of ARA, whereas they are more likely in cases of MCP, since the comparatively high rents of forest farming are a relevant incentive towards land cover change [45]. Policy measures aiming at implementing a decrease in LST should be based on targeting both ARA and MCP, as follows. A system of afforestation incentives should be implemented, targeting rural areas on the basis of agricultural rent, with the aim of encouraging low-rent farmers to turn into forest farmers. It is highly likely that these incentives are effective in cases of rural areas classified as either MCP or OPSP, whereas high-rent farming, which mainly takes place in ARA, should be less interested in implementing land cover changes [46]. From this perspective, a fundamental role has to be played by national, regional, and local administrations in order to identify the optimal size of afforestation-related land cover transitions and of public investment in terms of financial feasibility [47]. Public commitment to decrease LST in rural areas is very important as regards the effectiveness of planning policies, and it should be visible to the local communities through, for instance, the direct purchase of private rural land, such as low-rent croplands or abandoned, poorly vegetated open spaces, in order to develop afforestation processes [48].
The methodological approach defined and implemented in this study as regards the Sardinian region provides estimates of the quantitative size of the impact of public policies concerning urban and rural contexts. Moreover, the methodology defined and implemented in this study can be easily exported to other national and European regional contexts, since worldwide and free satellite images the enable the identification of the spatial distribution of LST are available, and the CORINE land-cover-based LEAC taxonomy is available for all European countries.
As regards the outcomes of this study, there are two main directions for future research in the field of spatial planning policies to decrease LST and mitigate the connected problems, which consist of heat waves and heat islands, and quality of life in rural and urban areas. On the one hand, a network of direct, on-site observations would be of fundamental importance as regards the validation of the spatial layout of LST, which would make feasible the experimental implementation of planning policies aimed at generating impacts on LST. No planning experiment could be developed in the absence of a validated baseline dataset related to LST. On the other hand, tests and pilot projects involving municipal administrations would be very useful to test the effectiveness of policies such as the afforestation or urban blue and green grids. Researchers in the field of sustainability-oriented planning and local public administrations should therefore lobby local, national, and international funding agencies so as to secure enough funds to implement experimental programs concerning decreasing LST, including pilot actions and local tests.