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

Abstract

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 modifications 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 find 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 confirm an increase in the LST as a consequence of the loss of vegetation cover, which favors the urban heat island phenomenon.

1. 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 main causes of urban heat island: (1) Anthropogenic heat production (combustion, etc.). (2) Reduction of latent heat loss (low degree of evapotranspiration) due to the use of impermeable surfaces (cement, concrete, pavements, etc.), efficient drainage systems, and scarce vegetation. (3) Low sensible heat loss due to reduced wind speed. (4) Greater absorption of solar radiation trapped by urban canyons. (5) Decrease in outgoing long-wave radiation due to a low sky vision factor. (6) Greater heat storage in building materials during the day due to their high capacity to retain heat and release it at night. (7) Re-emission of long-wave atmospheric radiation to the ground resulting from air pollution (greenhouse effect).

Several studies have shown the importance of landscape composition and UHI phenomenon, for example, a study in Washington showed average differences of about 10.8 °C between urban and rural areas [9]. Also, it has been observed that UHI phenomenon tends to have severe effects towards the interior of urban sites, dissipating towards rural areas [10,11,12]. Because of the decrease in evapotranspiration due to the construction of infrastructure, UHI phenomenon has consequences on climate variables on a local or regional scale. For example, variations in rainfall patterns [13] and surface temperature have been observed modifying the environmental and regional climate by remodeling the boundary layer and Earth–atmosphere circulation [14]. However, few studies have been carried out in tropical cities to understand the effect of land use change on land surface temperature (LST) and heat islands, particularly in ecosystems with marked seasonality. Seasonality represents one of the major sources of natural variation in the climate and is known to drive cycles of plant productivity in many ecosystems [15,16,17]. It has been widely observed that interannual variations in precipitation have direct implications on evapotranspiration, as well as on soil moisture [18] and solar radiation [19]. In addition, it has implications for the phenology of vegetation, affecting the roughness of the Earth’s surface, the albedo, and the interception of water [20], which leads to variations in LST.

The UHI phenomenon can be addressed from two approaches, atmospheric and surface [21,22]. Atmospheric UHI is assessed using data from meteorological stations, while surface UHI is studied using remote sensing thermal infrared (TIR) data [23,24]. In recent years, remote sensing techniques have been successfully applied to detect UHI phenomenon in several regions, analyze urban surface radiation, estimate surface heat balance, study evapotranspiration processes, etc. [25,26]. Different sensors provide data in the infrared and thermal, for example, among the most used for their spatial and temporal resolution, are the Moderate Resolution Imaging Spectroradiometer (MODIS) [27,28,29], and Landsat Thematic Mapper (TM)/Enhanced Thematic Mapper Plus (ETM+)/ Operational Land Imager (OLI)/ Thermal Infrared Sensor (TIRS ) [22,28,30,31], which allow the study of large areas with spatial continuity, as well as being open and free data [8,32]. In addition, using remote sensors, it is possible to construct indices, such as the Normalized Difference Vegetation Index (NDVI) and Normalized Difference Water Index (NDWI), among others, which make it possible to evaluate the condition of vegetation over large areas [33,34,35], as well as to know the distribution patterns of the UHI from surface temperature [30,36,37]. Moreover, few studies have linked the LST and the NDVI in karst regions with such conspicuous variations in seasonality. Therefore, it is essential to contribute to the understanding of the relationship between NDVI as a predictor of LST [38].

In Mexico, one of the regions with the highest population and urban growth has been the Yucatán Peninsula, particularly Mérida City. In the last three decades, Mérida has quadrupled its population—in 1970 there were 269,000, in 2020 there are estimated to be 1,161,000 people—which has significantly increased the urban area. In addition, in recent years Mérida has been awarded, by various means, as one of the best cities in the world to live (https://www.a-nah.com/article/why-is-merida-the-best-city-in-the-world accessed on 2 June 2020; https://www.eluniversal.com.mx/english/merida-declared-best-city-world accessed on 2 June 2020). Its high ecological, scenic, and cultural value associated with the presence of the Mayan culture, as well as its proximity to the coastal area, archaeological sites, cenotes, and caves, make it an attractive housing destination, 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.

2. Materials and Methods

2.1. 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].

2.2. 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. List of Landsat images used to carry out the analysis.

ID	Image Name	Path	Row	Year	Month	Day
	Dry Season
1	LT05_L1TP_020045_19950423_20170109_01_T1	20	45	1995	4	23
2	LT05_L1TP_020046_19950423_20170108_01_T1	20	46	1995	4	23
3	LE07_L1TP_020045_20000428_20170212_01_T1	20	45	2000	4	28
4	LE07_L1TP_020046_20000428_20170212_01_T1	20	46	2000	4	28
5	LE07_L1TP_020045_20020317_20170131_01_T1	20	45	2002	3	17
6	LE07_L1TP_020046_20020317_20170131_01_T1	20	46	2002	3	17
7	LC08_L1TP_020045_20140427_20170307_01_T1	20	45	2014	4	27
8	LC08_L1TP_020046_20140427_20170307_01_T1	20	46	2014	4	27
9	LC08_L1TP_020045_20150414_20170228_01_T1	20	45	2015	4	14
10	LC08_L1TP_020046_20150414_20170228_01_T1	20	46	2015	4	14
11	LC08_L1TP_020045_20160416_20170326_01_T1	20	45	2016	3	26
12	LC08_L1TP_020046_20160416_20170326_01_T1	20	46	2016	3	26
13	LC08_L1TP_020045_20180406_20180417_01_T1	20	45	2018	4	6
14	LC08_L1TP_020046_20180406_20180417_01_T1	20	46	2018	4	6
15	LC08_L1TP_020045_20190511_20190521_01_T1	20	45	2019	5	11
16	LC08_L1TP_020046_20190511_20190521_01_T1	20	46	2019	5	11
	Wet Season
17	LT05_L1TP_020045_19951203_20170105_01_T1	20	45	1995	12	3
18	LT05_L1TP_020046_19951203_20170106_01_T1	20	46	1995	12	3
19	LE07_L1TP_020045_20001106_20170209_01_T1	20	45	2000	11	6
20	LE07_L1TP_020046_20001106_20170209_01_T1	20	46	2000	11	6
21	LE07_L1TP_020045_20021230_20170127_01_T1	20	45	2002	12	30
22	LE07_L1TP_020046_20021230_20170127_01_T1	20	46	2002	12	30
23	LC08_L1TP_020045_20140105_20170307_01_T1	20	45	2014	1	5
24	LC08_L1TP_020046_20140105_20170307_01_T1	20	46	2014	1	5
25	LC08_L1TP_020045_20151108_20170225_01_T1	20	45	2015	11	8
26	LC08_L1TP_020046_20151108_20170225_01_T1	20	46	2015	11	8
27	LC08_L1TP_020045_20161228_20170314_01_T1	20	45	2016	12	18
28	LC08_L1TP_020046_20161228_20170314_01_T1	20	46	2016	12	18
29	LC08_L1TP_020045_20181202_20181211_01_T1	20	45	2018	12	2
30	LC08_L1TP_020046_20181202_20181211_01_T1	20	46	2018	12	2
31	LC08_L1TP_020045_20190119_20190201_01_T1	20	45	2019	1	19
32	LC08_L1TP_020046_20190119_20190201_01_T1	20	46	2019	1	19

). 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. Characteristics of remote sensors (https://landsat.gsfc.nasa.gov/appendix/references accessed on 2 June 2020).

Sensors	Characteristics of Sensors
	Bands		Spatial Resolution (m)	Wavelength (µm)
Landsat 5MT * and Landsat 7 ETM+ **	1	Blue	30	0.45–0.515
	2	Green	30	0.525–0.605
	3	Red	30	0.63–0.69
	4	NIR	30	0.775–0.90
	5	SWIR 1	30	1.55–1.75
	6	Thermal *	120	10.41–12.5
	6	Thermal **	60	10.4–12.5
Landsat 8 OLI	2	Blue	30	0.450–0.515
	3	Green	30	0.525–0.600
	4	Red	30	0.630–0.680
	5	NIR	30	0.845–0.885
	6	SWIR 1	30	1.560–1.660
	10	Thermal	100	10.6–11.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.

2.3. 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.

2.4. 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.

$$NDVI=\frac{\left(NIR-RED\right)}{\left(NIR+RED\right)}$$

(1)

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 “Bucureș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]:

$$L\lambda =ML\ast Q_{cal}+AL-O_i$$

(2)

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.

$$BT=\left(\frac{K_2}{\ln \left(\frac{K_1}{TOA}\right)+1}\right)-273.15$$

(3)

where K₁ and K₂ 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].

$$P_V=\left[\left(NDVI-NDV{I}_{min}\right)\right]/{\left[\left(NDV{I}_{max}-NDV{I}_{min}\right)\right]}^2$$

(4)

$$\epsilon =0.004\ast P_{V }+0.986$$

(5)

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].

$$LST=\left[BT/\left(1+\lambda \ast \frac{BT}{\varrho \ast \ln (\epsilon )}\right)\right]$$

(6)

where BT = brightness temperature; λ = wavelength of the radiance emitted in each band; ρ = constant of Boltzmann (1.38·10–23 J/K), and ε = surface emissivity) [61].

2.5. 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²). 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.

3. Results

3.1. 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. Changes in surface of the city and class cover.

	1995			2019
LSC	Area (ha)	Occupancy %	# Polygons	Area (ha)	Occupancy %	# Polygons
Urban	14,505.5	14.0	8577	24,215.66	23.1	5715
Scattered vegetation	18,151.1	17.1	6493	24,225.47	23.1	6509
Vegetation	72,343.9	68.9	3701	56,448.76	53.8	2821
Total	105,000.48	99.99	18,771	104,889.89	99.9	15,045

). 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. Kappa details.

Year	Class Cover	Total Precision [%]	Kappa Coefficient	Producer Accuracy [%]	User Accuracy [%]	Overall Accuracy [%]	Kappa Hat
1995	Urban	95.21	0.96	82.44	89.16	89.16	0.88
	Scatt. Veg.			89.08	89.11	89.11	0.87
	Veg.			98.61	97.49	97.49	0.92
2000	Urban	97.86	0.96	89.29	98.14	98.14	0.98
	Scatt. Veg.			95.28	95.03	95.03	0.94
	Veg.			98.54	98.85	98.85	0.96
2002	Urban	95.62	0.92	91.08	92.80	92.8	0.92
	Scatt. Veg.			89.61	95.14	95.14	0.94
	Veg.			98.17	95.98	95.98	0.9
2014	Urban	96.13	0.92	90.36	82.54	82.54	0.82
	Scatt. Veg.			85.34	94.59	94.59	0.94
	Veg.			99.48	97.92	97.92	0.93
2015	Urban	82.43	0.71	93.93	95.88	95.88	0.95
	Scatt. Veg.			98.27	79.25	79.25	0.56
	Veg.			30.33	74.16	74.16	0.68
2016	Urban	94.1	0.9	85.71	88.74	88.74	0.86
	Scatt. Veg.			97.99	97.71	97.71	0.97
	Veg.			100	99.41	99.41	0.99
2018	Urban	79.95	0.69	80.99	91.78	91.78	0.91
	Scatt. Veg.			55.86	92.53	92.53	0.89
	Veg.			99.88	69.77	69.77	0.49
2019	Urban	81.2	0.69	93.56	92.81	92.81	0.92
	Scatt. Veg.			61.91	63.78	63.78	0.52
	Veg.			84.76	84.38	84.38	0.66

shows the overall accuracy (OA), Kappa coefficient, and highest accuracy per cover class per year.

3.2. Historical LST Results

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 °C. Particularly, in the dry season, the highest temperatures were observed in 2014 and 2019 with an average of LST of 32.90 °C (18.3–47.5 °C) and 35.45 (23.1–47.8 °C), respectively. In the wet season, the year 2000 showed the highest average of LST with 28.25 °C (16–40.5 °C) (

Table 5. Changes in Normalized Difference Vegetation Index (NDVI) and land surface temperature (LST) of the metropolitan city of Mérida.

	NDVI						Temperature °C
	Dry		Wet				Dry		Wet
Year	Min	Ave	Max	Min	Ave	Max	Min	Ave	Max	Min	Ave	Max
1995	−0.35	0.29	0.93	−0.38	0.28	0.93	23.2	30.60	38.0	20.2	26.35	32.5
2000	−0.3	0.18	0.66	−0.53	0.14	0.81	19.7	33.55	47.4	16.0	28.25	40.5
2002	−0.35	0.21	0.77	−0.42	0.17	0.76	13.8	30.05	46.3	13.3	26.90	40.5
2014	−0.65	0.07	0.79	−0.88	0.01	0.9	18.3	32.90	47.5	18.0	25.85	33.7
2015	−0.67	0.07	0.81	−0.81	0.03	0.87	18.7	30.05	41.4	5.7	19.30	32.9
2016	−0.53	0.145	0.82	−0.86	0.01	0.87	21.6	30.90	40.2	16.7	25.00	33.3
2018	−0.69	0.045	0.78	−0.2	0.18	0.56	21.1	31.25	41.4	19.4	25.45	31.5
2019	−0.73	0.04	0.81	−0.74	0.03	0.8	23.1	35.45	47.8	14.5	14.50	33.3

). 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).

3.3. Historical NDVI Results

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).

3.4. 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² = 0.2511 (RMSE = 2.32 °C); while for the northern season the fit was slightly lower, R² = 0.1476 (RMSE = 1.57 °C).

4. 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].

5. 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.