Impact of Urban Land-Cover Changes on the Spatial-Temporal Land Surface Temperature in a Tropical City of Mexico

: Climate change has severe consequences on ecosystem processes, as well as on people’s quality of life. It has been suggested that the loss of vegetation cover increases the land surface temperature (LST) due to modiﬁcations in biogeochemical patterns, generating a phenomenon known as “urban heat island” (UHI). The aim of this work was to analyze the effects of urban land-cover changes on the spatiotemporal variation of surface temperature in the tropical city of M é rida, Mexico. To ﬁnd these effects we used both detected land-cover changes as well as variations of the Normalized Difference Vegetation Index (NDVI). M é rida is ranked worldwide as one of the best cities to live due to its quality of life. Data from satellite images of Landsat were analyzed to calculate land use change (LUC), LST, and NDVI. LST increased ca. 4 ◦ C in the dry season and 3 ◦ C in the wet season because of the LUC. In addition, a positive relationship between the LST and the NDVI was observed mainly in the dry season. The results conﬁrm an increase in the LST as a consequence of the loss of vegetation cover, which favors the urban heat island phenomenon.


Introduction
Urbanization is a global phenomenon widely associated with processes of land-cover change with negative consequences on the ecosystem functioning. Particularly, under climate change context, the consequences of urbanization on local weather variables (e.g., rainfall intensity and surface temperature) are increasingly important [1]. During the last three decades, the world's urbanization has seen a severe transformation in urban landscapes [2]. According to a United Nations report [3], urban populations increased from 751 million people in 1950 to 4.2 billion in 2018. In this context, Latin America is the second region, after Africa, with the highest rate of urban growth, approximately 133 million people live in urban areas. In Mexico, approximately 90 million people (77% of the population) live in urban areas [4].
One of the main negative consequences of the increase in urban infrastructure is the creation of impermeable and dry surfaces that reduce the availability of moisture for natural evapotranspiration from the soil [5], causing a phenomenon called urban heat island UHI [5][6][7]. The UHI is defined as the difference in temperature between the hottest sector of the city and the non-urban, or rural, space surrounding the city at a given time and is influenced by climatic factors (solar radiation, wind and clouds, for example), and the variability in structures and materials typical of an urban area [6][7][8] define seven which has increased the demand for housing and the inherent growth of the urban area. However, this population growth, as well as poor urban landscape planning, can put the quality of life, as well as the ecosystem functioning of the region, at risk.
In summary, the LST can be a proxy to know the health status of ecosystems. The LST provides information to explore physical and chemical processes at the surface, as well as to understand the variation in energy flows due to land use change (LUC), deforestation, and environmental degradation. Furthermore, under a context of accelerated urban growth, as is the case of the city of Mérida, Yucatán, Mexico, there is an increase in the intensity of human activities, as well as in land use change. Therefore, the main objective of this research was to analyze the effects of land-cover changes on land surface temperatures during the rainy and dry seasons, using remotely sensed tools. To accomplish this goal, we have the four following objectives: (1) To determine the urban growth rate of the city of Mérida over the last 24 years. (2) To evaluate if there is a variation of land surface temperature after a land-cover change happened. (3) To assess if this variation in temperature is attributable to a loss of land-cover vegetation. (4) To analyze the relationship between NDVI and land surface temperature among cover and seasonality classes.

Study Area
The study was carried out in a mosaic of different uses and land covers 1,051,011.4 ha of human impacted landscape, located in the municipalities of Mérida, Conkal, and Kanasin, Mexico (coordinates: 21 • 11'18.42" N; 89 • 43'17.76" W and 20 • 41'50.02" N; 89 • 39'8.02" W) ( Figure 1). The region is characterized by a calcareous platform, with karstic soils of biogenic origin, and absence of major topographic features, as well as surface rivers [39]. The combination of geomorphological, climatic, hydrological, soil, and ecological elements of the region encourages a particular development of biodiversity. The climate of the city of Mérida is the driest of subhumid climates with summer rains Ax'(w)w, a high percentage of winter rains and intraannual drought [40]. It is also a region prone to tropical cyclones during the rainy season, and cold fronts during the winter months due to its proximity to the sea and the effects of an ocean current from the Yucatán Channel [41,42]. Two climatic seasons were identified to this study, dry season (March-May) and wet season (July to February); the study area is the driest in the Yucatán Peninsula with an average annual rainfall of 500 mm [43,44]. According to the classification of Miranda and Hernandez-X [45], the dominant vegetation type is tropical dry forest. However, anthropogenic activities such as cattle raising, agriculture, open-pit mining, and the growth of the urban patch modified the vegetation matrix of the tropical dry forest towards secondary succession with different ages of abandonment [46].

Data Acquisition and Image Preprocessing
The effect of the loss of vegetation cover on the land surface temperature (LST), as a consequence of the increase in urban infrastructure, was evaluated based on a historical analysis of land-cover changes in the municipalities of Conkal, Kanasin, and Mérida. Images were acquired from Landsat-5(MSS/TM), Landsat-7 (ETM+) and Landsat-8 (OLI/TIRS) sensors (USGS; http://www.usgs.gov/ accessed on 2 June 2020) for a period of 24 years (1995, 2000-2001, 2014-2016, and 2018-2019). The analyses were carried out for dry season (February-May) and wet season (November-January) ( Table 1). The spectral and spatial resolution of the images varied between the sensors; Landsat TM and ETM+ have seven bands with the same spatial resolution (30 m) and a thermal band of 100 and 60 m, respectively; while Landsat 8 has nine bands at 30 m resolution (OLI), one thermal band at 100 m spatial resolution (TIRS), and one panchromatic band at 15 m (Table 2). Particularly, the following bands were used for the LST calculation: for the Landsat 5 sensor, band 6 was used; while for Landsat 7, band 6_VCID_1 was utilized, here we named as band 6. This last band is recommended for tropical studies because it does not present radiometric saturation in the hottest areas [47]. For Landsat 8, we selected band 10 instead of band 11 to calculate land surface temperature. This is because band 10 has the best adjustment in terms of temperature calculation, while band 11 presents an overestimation [48]. Based on the suggestions of Chuvieco [49], the acquired images were radiometrically and atmospherically corrected to reduce the effect of atmospheric scattering using Semi-Automatic Classification Plugin (SCP) Version 7.0.11-Metera [50], in QGIS Version 3.12 "București" software [51]. This software is an open-source geographic information system for all operative systems. Both QGIS and SCP are widely used for the study of the Earth, supported by remote sensing [52][53][54][55][56]. The general schematic workflow for this study is shown in Figure 2.

Calculation of Land-Cover Changes
To detect land-cover changes, a supervised classification was performed using the dzetsaka classification tool plugin, with Gaussian Mixture Model [57], with corrected Landsat images for the years 1995, 2000-2001, 2014-2016, and 2018-2019. The years were selected based on two criteria: (1) the end of the 1990s represented an important point in the "modernization" of the state of Yucatán, and Mérida began to have relevance at the national level; (2) the availability of satellite images. The regions of interest (ROI) for each of the classifications were made based on our experience in the area of study and vector information from the National Institute of Statistics and Geography of Mexico [58], as well as with the support of the Google Earth tool. We determined 100 ROI per type of coverage, giving a total of 300 ROI per year. The types of coverage defined were as follows: urban (e.g., houses, buildings, roads, and avenues), vegetation, and bare soil with dispersed vegetation (thereafter, scattered vegetation). Subsequently, the thematic accuracy of each of the classifications was assessed using overall precision, and the accuracy of the classifications was calculated using the Kappa index [59] based on an independent subset of 30 ROI per cover class. A 3-pixel by 3-pixel filter was applied to remove the salt-and-pepper effect from the monitored ratings. Changes of class covers were calculated through a confusion matrix.

Computation of Normalized Differences Vegetation Index (NDVI) and Land Surface Temperature (LST)
NDVI was calculated following the general normalized difference between band TM4 (near infrared-NIR) and band TM3 (visible red-RED) from the Landsat preprocessed images Equation (1). NDVI is a key element to calculating, afterward, the proportion of the vegetation (P V ) and emissivity (ε) necessary to calculated LST.
LST calculation is a key component for understanding energy transfer models between the soil-vegetation and the atmosphere [60]. In accordance with the methodology proposed by Avdan and Jovanovska [61] to estimate LST, we used data of the thermal infrared region of each sensor through Semi-Automatic Classification Plugin Version 7.0.11-Metera [50], in the free and open-source software QGIS Version 3.12 "Bucures , ti" [51] (Table 1). First, the spectral radiance (Lλ) of the top of the atmosphere (TOA) was calculated by introducing band 10 and equations taken from the USGS [62]: where ML represents a multiplicative rescaling factor of the specific band, Q cal is the value of the thermal band image, AL is the rescaling specific factor of the band, and O i is the correction of the thermal band. The next step was to calculate the brightness temperature (BT) from the absolute radiance values.
where K 1 and K 2 represent the band-specific thermal conversion constants from metadata. Furthermore, the results were converted from Kelvin to Celsius degrees, adding the absolute zero. In addition, the proportion of vegetation in each pixel was calculated, expressed as P v from NDVI, according to Equation (4); from this, the land surface emissivity was calculated (Equation (5)). Finally, the real LST was estimated according to Equation (6) [61,63].
where 0.004 corresponds to the average emissivities value of bare soil, 0.986 the average emissivity values of the vegetated areas, and P v is the proportion of vegetation [63].

Statistical Analysis
To determine the influence of the change in coverage on the LST, the years 1995 and 2019 were selected. The LST of the vegetation cover class of the year 1995 was used to compare with the LST of all those pixels that changed cover in 2019; the analysis was carried out using the Student's t-test. Furthermore, to determine if there were statistically significant differences in the mean of the LST between the cover classes, an analysis of the unidirectional variance (ANOVA) by season was performed. The dependent variables were formally tested for normality and homoscedasticity, while the response variables were transformed as needed with 1/x, log10(x), log10(x + 1), and sqrt(x) to meet linearity assumptions. Finally, to establish the relationship between the LST and the NDVI, the values of temperature and the NDVI were extracted for each of the coverage classes. Subsequently, a regression analysis with cross-validation by station was carried out, and the adjustment was compared using the coefficient of determination (R 2 ). All analyses were performed with the glm2 package [64] in the R [65] software. The general schematic workflow for this study is shown in Figure 2.

Land-Cover Changes of Mérida
The total analyzed area was 105,008.3 ha; between 1995 and 2019, 31,679 ha of landscape (30%) changed to another cover class, while 73,209 ha (70%) remained unchanged ( Table 3). The urban class registered the greatest increase from 14,505 to 24,215 ha, followed by the scattered vegetation class, which showed an increase from 18,151 ha to 24,225. In contrast, the areas of vigorous vegetation showed a negative trend, as they went from 72,343 to 56,448 ha. The same classification method was carried out in all the Landsat images, and Table 4 shows the overall accuracy (OA), Kappa coefficient, and highest accuracy per cover class per year.   Figure 3A shows the historical spatial distribution of LST in Mérida, calculated from the Landsat images from 1995 to 2019 by dry and wet seasons. In general, an increase in the LST was observed from 1995 to 2019 in both seasons; however, the increase is more evident in the wet season with an increase of approximately 2 to 3  Table 5). The LST values showed significant differences between the coverage classes for dry and wet seasons (F = 8532, p = 0.0001; F = 50,636, p = 0.0001, respectively). In the dry season, a greater variation in the LST was observed in the scattered vegetation and vegetation cover classes, contrary to what was observed in the wet season. All classes showed an increase between approximately 8 and 10 • C between dry and wet season (Figure 4).  Table 5). The LST values showed significant differences between the coverage classes for dry and wet seasons (F = 8532, p = 0.0001; F = 50,636, p = 0.0001, respectively). In the dry season, a greater variation in the LST was observed in the scattered vegetation and vegetation cover classes, contrary to what was observed in the wet season. All classes showed an increase between approximately 8 and 10 °C between dry and wet season (Figure 4).      Figure 3B shows historical spatial distribution of NDVI values in Mérida. In general, a decreasing pattern in NDVI values is observed in the last 24 years, suggesting a decrease in vegetation vigor, as well as a consequence of the change in land use. In particular, it was observed that the lowest NDVI values are located in the area corresponding to the center of the city, as opposed to the periphery of the study area, where the highest values are located. The variation of the NDVI values in the whole study area is between -0.3 and 0.93 in the dry season and between -0.2 and 0.93 in the wet season. Table 5 shows the median, as well as the maximum and minimum values of the NDVI, and LST per season per year. In terms of coverage classes, significant differences were observed between them ( Figure 5). Statistically significant differences were observed in the NDVI values with respect to the coverage classes in the dry season (F = 268,025, p = 0.0001), as well as in the wet season (F = 395,166, p = 0.0001) ( Figure 5).   Figure 3B shows historical spatial distribution of NDVI values in Mérida. In general, a decreasing pattern in NDVI values is observed in the last 24 years, suggesting a decrease in vegetation vigor, as well as a consequence of the change in land use. In particular, it was observed that the lowest NDVI values are located in the area corresponding to the center of the city, as opposed to the periphery of the study area, where the highest values are located. The variation of the NDVI values in the whole study area is between -0.3 and 0.93 in the dry season and between -0.2 and 0.93 in the wet season. Table 5 shows the median, as well as the maximum and minimum values of the NDVI, and LST per season per year. In terms of coverage classes, significant differences were observed between them ( Figure 5). Statistically significant differences were observed in the NDVI values with respect to the coverage classes in the dry season (F = 268,025, p = 0.0001), as well as in the wet season (F = 395,166, p = 0.0001) ( Figure 5).  Figure 3B shows historical spatial distribution of NDVI values in Mérida. In general, a decreasing pattern in NDVI values is observed in the last 24 years, suggesting a decrease in vegetation vigor, as well as a consequence of the change in land use. In particular, it was observed that the lowest NDVI values are located in the area corresponding to the center of the city, as opposed to the periphery of the study area, where the highest values are located. The variation of the NDVI values in the whole study area is between -0.3 and 0.93 in the dry season and between -0.2 and 0.93 in the wet season. Table 5 shows the median, as well as the maximum and minimum values of the NDVI, and LST per season per year. In terms of coverage classes, significant differences were observed between them ( Figure 5). Statistically significant differences were observed in the NDVI values with respect to the coverage classes in the dry season (F = 268,025, p = 0.0001), as well as in the wet season (F = 395,166, p = 0.0001) ( Figure 5).

Relationship between LST, Land Use Change (LUC), and NDVI
Significant differences were found in the mean LST between 1995 (29.36 • C) and 2019 (33.89 • C) (t, 201,651 = 677.02, p = 0.0001), for both the dry and rainy seasons (29.36 • C and 33.89 • C, respectively; (t, 211,288 = 1041.3, p = 0.0001) reported for the city of Mérida, Yucatán. Furthermore, statistically significant differences were found in the mean LST in those pixels that went from vegetation cover class in 1995 to urban in 2019 for dry (t, 20,245 = 6.99, p = 0.0001) and for northern (t, 16,294 = 41.42, p = 0.0001). Regression analysis showed a positive relationship between LST and NDVI using 80% of the data for training and 20% for validation. The model that presented the best fit for estimating the LST from the NDVI was that of the dry season with an R 2 = 0.2511 (RMSE = 2.32 • C); while for the northern season the fit was slightly lower, R 2 = 0.1476 (RMSE = 1.57 • C).

Discussion
Mérida City is historically recognized as a high cultural value city due to its link with the Mayan civilization. In recent years it has received several recognitions as one of the cities with the best quality of life in the world. Consequently, the real estate supply has intensified in the last decade, being one of the most profitable and productive economic sectors at the regional and national level [66][67][68].
Since its foundation (1953), Mérida had a configuration that favored green and treed areas; however, in the last 10 years it has experienced accelerated horizontal growth [66] towards the periphery, mainly in the eastern and western zones, under a construction model that bets on a spatial arrangement that seeks to increase the number of homes, sacrificing green spaces [67,68].
The main consequences of urban growth are the imminent loss or degradation of natural vegetation cover. In the case of Mérida, the horizontal growth of the city has implied the absorption of small settlements, transformation of areas with scattered vegetation, and reduction of the natural vegetation cover of 15,895 ha in the last 24 years. According to Sullivan and Ward [69], this type of growth in urban areas causes significant loss of green areas, agricultural land, in addition to increased levels of consumerism and larger traveling distances in people's daily lives [67,70].
NDVI values obtained suggest a global decrease in vegetation vigor in the region, as a consequence of the land use change in the last 24 years. Particularly for Mérida and its metropolitan area, the NDVI values estimated are between 0.2-0.5. Furthermore, the values observed within the city suggest that the current vegetation is not vigorous enough to maintain ecosystem services and that the observed mix between urban areas and vegetation suggests areas of environmental stress caused by edge effects. Regarding the vegetation outside the urban zone, different density and coverage conditions were observed (NDVI: 0.6-1) in a clear gradient of vegetation degradation that decreases with distance to the urban zone, similar to an edge effect. Finally, the correlation observed between cover classes and NDVI values in this study is in agreement with that found by Ferreira and Duarte [71] in shrub and bush species, suggesting that most of the vegetation is deciduous [46].
Degradation, or loss of vegetation, alter the balance of solar radiation that falls on the Earth's surface, that modifies the distribution of humidity and heat in the atmosphere, resulting in a local increase in temperature [5,71,72]. For the study area, in the last two decades, 15,895 ha of natural vegetation have been lost; 9710 ha of these have been replaced by artificial materials, such as concrete, cement, or pavement, with a significant impact on the modification of the local temperature [68].
The results show that the highest temperatures are focused into the Mérida urban zone, and they have increased over time, contrary to what was observed in peripheral areas that maintain vegetation cover. These results suggests that climate change processes occur at different spatial scales. At local scales, the increase in urban infrastructure increases albedo, decreases evapotranspiration, and generates the UHI phenomenon [5,71,72], while at regional scale, forest cover keeps the temperature relatively constant, and the climatic variations could be associated with regional or global climatic events or phenomena (El Niño and the Southern Oscillation, or global climate change [73][74][75]). This is similar to other tropical cities as Beijing, Changchun, and Hangzhou, China [1,8,11], Cameron Highlands, Malaysia [76], or Shiraz, Iran [77]. In all cases, UHI phenomena are associated with the extension of the urban area, type of structures and material, density housing, and presence of areas with continuous vegetation within the city.
In the case of Mérida, the difference in LST between urban areas and the peripheral areas with vegetation (scattered or dense) is on average 3.5 • C at the beginning and end of the study season, and is maintained in both seasons, similar to that observed by Adeyeri et al. [78], Ferreira and Duarte [71], and Dang et al. [79], among other authors.
Particularly in the urban zone, high temperatures were observed in the periphery (east and west) of the city, this is associated with concrete, pavement surfaces, and bare soil on the city's edge. In contrast, the north of the city, where residential houses and greater treed areas are, presents a cooler condition, in accordance with that reported by Pérez-Medina and López-Falfán [68] for the study area, and by Ferreira and Duarte [71], Martín-Vide and Moreno-García [7], and Wang et al. [80], among others, for similar studies.
Likewise, specific sites with extreme temperatures were identified, associated with the cement industry, the sanitary landfill, fires, and burned areas as part of agricultural activities, in the same way it was observed that, in recent years, the temperature between these hot spots and the rest of city is increasingly similar, which shows the occurrence of the UHI. In big cities where the UHI phenomenon has been recorded, the heat island is concentrated in the center of the cities [22,24,78,79,81]. However, for Mérida's urban zone, the heat island is not located in the center, but in the periphery, where the new developments are located, as pointed out by Pérez-Medina and López-Falfán [68]. With regard to the lower temperatures, they are distributed on the city outskirts, in areas with the presence of water (artificial lakes) and areas with dense, natural, or permanently irrigated vegetation, such as agricultural fields or the golf club.
The obtained results confirm the importance of green and treed areas, natural or artificial, to help reduce the temperature of the surface and the air [46,78], an issue of great relevance for a region with dry subhumid climate, such as Mérida [40], where the average annual temperature is 46 • C in the dry season, and 34 • C in the wet season.
In dry tropical cities such as Mérida, the climatic conditions accentuate the perception of the UHI effect, which reduces the comfort level of the inhabitants, increases the dependence on cooling, favors the presence of pollutants and health problems, among other effects [5,68]. In cities with projections of rapid expansion and extreme hot weather, such as Mérida, it is critical to consider, in the planning of new housing developments, an adequate distribution and density of vegetation as well as promoting the creation of green spaces within the current settlements, to mitigate and prevent the UHI phenomenon [26,82,83], among countless other environmental, social, and economic benefits [8,67,68].

Conclusions
The recognition of Mérida as one of the cities with the best quality of life in the world has been a driving force behind the demand for residential spaces in recent decades. The degradation and loss of vegetation in Mérida City and its surroundings has led to the modification of the microclimate, accentuating the high temperatures in the region within the urban area.
The extensive areas with vegetation, natural or irrigated, inside and outside the city, contribute to mitigating the warm natural temperatures in this dry tropical city, mainly in the dry season. The hotspot identification highlights the need to promote the conservation of areas with natural vegetation, and to integrate green spaces in the planning of new housing developments, to mitigate the UHI phenomenon, reduce energy consumption and greenhouse gas emissions, raise the comfort level of the inhabitants, and contribute to the sustainability of the city.