The Gontougo Region is one of the 31 regions of Ivory Coast. Since its establishment in 2011, it has been one of two regions in Zanzan District. The seat of the region is Bondoukou, and its area is 16,100 km². The population in the 2021 census was 917,828. Gontougo is currently divided into five departments: Bondoukou, Koun-Fao, Sandégué, Tanda, and Transua (RGPH-CI, 2021). It is located between 4°00′ and 2°50′ W longitudes and 6°30′ and 8°50′ N latitudes (Figure 1).

The soils of the region are generally ferritic, more or less desaturated, deep sandy-clay, and very suitable for cash crops (cashew, cocoa, coffee, and rubber) and food crops (e.g., yams, cassava, rice, maize, and tomatoes) (République de Côte d’Ivoire, 2019). According to the climatic division of Côte d'Ivoire, there are three climatic zones (Kouame, 2021; Kouame et al., 2020): northern, central, and southern zones. The central zone, of which the study area is a part, presents a bimodal seasonal cycle with less pronounced rainfall maximums in June and September. It has four seasons: a long rainy season with frequent rainfall and numerous thunderstorms from April to July, a small dry season in which the sky can remain overcast from August to September, a short rainy season with some small rainfall from September to November, and a long dry season from December to March.

The present study focused on satellite images and field training reference data collected between November and December 2022 in the Gontougo Region to examine LULC change and assess the correlation between NDVI, LULC, and LST. Satellite images include Landsat-8 Oli/Tirs and Sentinel-2 Msi.

According to America’s space agency [Nasa Earth Observation (https://earthobservatory.nasa.gov/)], LST is how hot the “surface” of the Earth would feel to a touch in a particular location, and NDVI is used to quantify vegetation greenness and density and is useful in understanding plant health. Calculations of NDVI for a given pixel always result in a number that ranges from minus one (−1) to plus one (+1); however, no green leaves give a value close to zero. A zero means no vegetation, and close to +1 (0.8–0.9) indicates the highest possible density of green leaves.

The overall change rate (Tg) is used to estimate the overall progress (the proportion of gain or loss) of the areas of LULC units (Dahan et al., 2021; Aka et al., 2022). It is obtained from the following mathematical formula: T g = S 2   –   S 1 /   S 1 * 100 , (7) where Tg is the overall rate of change (%); S1 represents the area of the class on date t₁ (initial date); S2 is the area of the class on the date t₂ (final date), and t₂ ˃ t₁.

We used the Pearson correlation coefficient, also known as the Pearson R statistical test, to measure the strength between the different variables and their relationships (Eq. (8)). Here, we measured a statistical test between two variables (NDVI and LST). Therefore, according to previous studies (Zheng et al., 2014; Peng et al., 2020), the calculation of the correlation coefficient value enables comprehension of the strength of the relationship between the two variables. Pearson correlation coefficient can range from +1 to −1, where +1 indicates the perfect positive relationship between the variables considered, −1 indicates the perfect negative relationship between the variables considered, and 0 indicates no relationship between the variables considered (Zheng et al., 2014; Njoku and Tenenbaum, 2022): r x y = n ∑ x i y i − ∑ x i ∑ y i n ∑ x i 2 − ∑ x i 2 n ∑ y i 2 − ∑ y i 2 , (8) where r _xy is Pearson r correlation coefficient between x and y, n represents the number of observations, x _i means the value of x (for ith observation), and y _i indicates the value of y (for ith observation).

Through qualitative and quantitative comparison of the change, the LULC types were compared with the LST values over the classification years based on their respective changes (Moisa et al., 2022a).

Over the study period, between 2015 and 2022, the NDVI of the Gontougo Region was determined according to the NIR and Red multispectral bands of the Landsat images. We noticed that the year 2021 presented the highest NDVI mean value with 0.246 and the year 2016 presented the lowest with 0.205 (Table 2). Each NDVI value per year of study was statistically presented in a histogram. These histograms show the distribution of vegetation cover according to the channel in which the NDVI value is high, medium, or low (Figure 3). The distribution of NDVI revolves around the mean values where the highest peaks are observed.

Qualitatively, the results show that the northeast, northwest, and west parts are much more devoid of vegetation than the other parts of the study area. The central and southern areas have some greenery. Overall, the values vary from −0.108 to 0.566 (Figure 4). The mean values of each NDVI were extracted to appreciate the trend of the vegetation (Figure 5), and the years 2019–2022 showed higher average vegetation than the years 2015–2018, in which the trend was lower.

The LST is used to monitor the temperature and surface processing of land features. In the present study, LST is calculated using thermal bands from Landsat images from 2015 to 2022. Table 5 reveals that the maximum temperature peaked at 50.11°C in 2021, whereas the minimum temperature was recorded at 9.21°C in 2019. Overall, LSTs are around 31°C when considering mean values. Similarly, at the mean level, there has been a temperature decrease of 0.19°C from 2015 to 2022 (Table 3).

To better appreciate these values, we generated histograms for each year of LST. These histograms show the distribution of LST values. The horizontal axis shows the minimum values at the beginning and middle, and the mean and maximum values at the end. On the vertical axis, we can see the evolution, especially the concentration of the LST values from the diagonal line (Figure 6). For all years, the distribution of LST is more or less concentrated on the mean in which the highest peaks are observed.

The findings of the study revealed that the northeastern and some western parts of the study area had high LST except in the year 2019, when the western part was not high, whereas the southern and center-eastern parts had low LST (Figure 7).

Based on the mean values, we represented the trend of the LST. The year 2017 recorded the lowest temperature at 30.47°C, and the year 2016 had the highest temperature at 32.73°C. The overall trend shows a constant and slight evolution along the LST trend line (Figure 8).

All the results in this section were obtained from the classification of Sentinel-2 images using GEE. Therefore, we used the following abbreviations: UA, user accuracy; PA, producer accuracy; OA, overall accuracy; Acrops, annual crops; Pcrops, perennial crops; Sav, savannah; Forest, forest; BarSet, bareland and settlement; Wb, waterbody.

The relationship between LST and NDVI was examined using Pearson correlation analysis. Pearson correlation is a measure of the linear relationship between two variables. It ranges from −1 to 1, where −1 indicates a perfectly negative correlation, 0 indicates no correlation, and one indicates a perfectly positive correlation. It is used in statistics to examine the strength and direction of the relationship between two variables (Zheng et al., 2014; Abulibdeh, 2021). We performed statistical analysis by comparing the LST and NDVI values after reclassification in which both variables had the same number of classes. We had eight images on each side; therefore, eight matches were made per pair (LST and NDVI together). The LST values ranged from maximum (50.11°C) to minimum (9.21°C) temperatures, whereas the NDVI values ranged from 0.567 to −0.108, from maximum to minimum. The results show that LST and NDVI have a strong negative correlation each year (Table 10). Following Table 12, which shows the type of correlation between NDVI and LST, we investigated the plots generated by these tables, which show how the regression line declines at each correlation (Figure 12). This indicates a negative correlation between NDVI and LST. Hence, the Pearson correlation coefficient (r) recorded at each correlation indicates negative values from 2015 to 2022. These values are −0.992, −0.982, −0.975, −0.977, −0.966, −0.994, −0.990, and −0.978, respectively, from 2015 to 2022.

The distributions of spatial LST values for the different LULC classes were reclassified and mapped for the two LULC years, 2015 and 2022, to depict temperature differences on the maps and show the LULC areas that are more or less heated (Figure 13). These LST values reclassify maps to show that the hottest areas are the savannah and the bareland and settlement. The results are similar for the low-temperature areas, which are more or less forest areas. In addition, there is a significant spread of LULC classes and temperatures that vary immensely over the whole area on both maps. Furthermore, the surface temperatures varied between 17.33°C and 45.22°C in 2015 and 16.53°C and 45.81°C in 2022. These values can be observed on the map through the legend.

Subsequently, the mean LST values were matched to the LULC classes. In 2015, the bareland and settlement class recorded the highest surface temperature (35.18°C), followed by the Savannah class (34.60°C). The forest class recorded the lowest temperature (28.69°C) (Table 11). In 2022, still with the mean LST values, the correlation showed again that the bareland and settlement class recorded the highest LST at 35.41°C, and the lowest LST was recorded for the forest class (28.46°C). The differences in mean temperature are a decrease in LST of 0.23 for the forest class and an increase for the bareland and waterbody classes.

The spatiotemporal analysis of the Gontougo LULC change made it possible to categorize six LULC classes each year: annual crops, perennial crops, savannah, forest, bareland and settlement, and waterbody. The OA obtained is 93.77% with a kappa coefficient of 0.87 for the year 2015 map and 96.01% OA and 0.92 kappa coefficient for the year 2022 map. These results are close to those obtained by Kadio et al. (2022), who mapped the Sentinel-2 image and compared it with Landsat-8 Oli in the Azagny site Ramsar in Côte d’Ivoire. The OA was 90.63%, and the kappa was 0.94. In addition, Zhang et al. (2023) used GEE and RF with three Sentinel-2 images in the crop mapping type for agronomic purposes to quantify cultivated areas and determine yields in China. They obtained 95.72%, 97.21%, and 98.13%, respectively, as OA and 0.94, 0.96, and 0.98 as kappa for the years 2017, 2018, and 2020. This latter study is somewhat similar to ours in the sense that it discusses the use of ML and RFA and highlights all the cartographic precision obtained in this study. In remote sensing, a classification is deemed acceptable when the kappa coefficient value is greater than 75% (Jensen J. R., 1986; Congalton, 1991; Zhang et al., 2023). It confirms that the results of this analysis are statistically acceptable. Currently, RFA is widely used in satellite image classification because it gives better results (Teluguntla et al., 2018; Phalke et al., 2020; Yan et al., 2022; Zhao et al., 2023). We are increasingly turning to ML technology, such as GEE, to implement LULC change. These authors also performed the same in their LULC study using Landsat collection images (Teluguntla et al., 2018; Phalke et al., 2020; Suryono et al., 2021; Yan et al., 2022; Ashane et al., 2023; Zhao et al., 2023). The particularity of our study here is to use Sentinel-2, which is one of the studies that has so far mapped the LULC change in the entire Gontougo Region and assessed the LST using ML technology (GEE) and RFA. However, it should be noted that Sentinel-2A was launched in June 2015, and the first images were available on 23 June 2015. These early images of the Sentinel-2 mission over our study area showed quite a bit of artifact and constituted bits of images with continuous haze. The mosaicking and speckle cleaning gave us a homogeneous image from which the 2015 classification was made. This justified the low precision obtained compared to 2022, when the filter interval search was conducted over a full year.

The other striking aspect of our research is the drastic loss of forests and the dramatic decrease in annual crop areas. Indeed, from 2015 to 2022, forests lost 62.19% of their area, and annual crops lost 86.6% of their area. These various losses benefit perennial crops by up to 48.13% for the forest and 65.14% for annual crops. Indeed, according to a study conducted in the region on the re-profiling of roads by the State of Cote d'Ivoire, and according to the latest census of the population and habitat, the Gontougo Region’s potential economy is based primarily on agriculture, and the reputation of this region is based on the famous variety of yam called “Kponan” (République de Côte d’Ivoire, 2019; RGPH-CI, 2021). The current situation shows a landscape strongly dominated by perennial crops, with cashew nuts in the lead. The region is one of the largest cashew producers in Côte d’Ivoire, with 65,018 tons in 2022, according to the Cashew Cotton Council (CCA). In detail, at the level of each Gontougo department, it is 49,735 t in Bondoukou, 3,315 t in Koun Fao, 7,725 t in Sandégué, 4,094 t in Tanda, and 194 t in Transua. These statistics are given by the Agence Ivoirienne de Presse (AIP) in an article published on 3 February 2023 at 8:01 GMT entitled «Côte d'Ivoire—AIP/Noix de cajou: L'Indénié-djuablin et le Gontougo ont réalisé une production annuelle de 68,697 tonnes en 2022» and edited by Memel Franck Niagnely¹. This necessitates the availability of a large production area and explains why the vegetation cover and yearly crop area have been drastically reduced. As for the yam crop, in most cases, it is used to prepare or lay out the soil for the next cashew season. According to the yam farmers themselves, they transplant small cashew seedlings into their producing yam fields. After the yam harvest, the cashew field takes over completely, and new annual crop areas are sought for the following year. This is one of the factors contributing to the annual decline in the area planted for annual crops. The fact that annual crops are not grown in isolation and that it is important to leave trees in the fields of annual crops, particularly yam, to encourage the growth of yam buds, also serves to justify the low precision of mapping technically and discrimination of the annual crops class ¹ .

Land surface heating is the cumulative effect of LULC change and global atmospheric warming (Moisa et al., 2022b; Moisa et al., 2022c). The decline in vegetation cover driven by agricultural expansion in the study area substantially increases the LST. Indeed, the results show that the lowest observed temperatures are in areas where forests are present. Respectively, the LULC forest areas of 2015 and 2022 recorded the lowest mean LST of 28.69°C and 28.46°C. This translates into a slight decrease in temperature over the year. This trend is verified within NDVI calculations. According to some authors (Carlson and Riziley, 1997; Davi et al., 2006; Morawitz et al., 2006; Du et al., 2017; Wu et al., 2019; Syawalina et al., 2022), healthy vegetation absorbs most of the visible light that hits it and reflects a large portion of the near-infrared light. Unhealthy or sparse vegetation reflects more visible and less near-infrared light (Davi et al., 2006; Dibi N’da et al., 2008; Khan et al., 2010; Phan et al., 2021). Therefore, NDVI has a huge influence on LST calculation and the surface temperature content. Therefore, we used it to correlate with the LST. Similar to the LST, NDVI was reclassified into four classes and was well-corrected with the LST (Table 12). This NDVI classification corresponds to the vegetation standards from sub-Saharan humid tropical countries because the definition of the forest depends on the forestry code of each country (Dibi N’da et al., 2008; RCI-Code Forestier, 2019b; RCI-Code Forestier, 2019a). Accordingly, the negative NDVI values corresponding to the “non-vegetation land” description recorded the highest temperatures. This can also be perceived with the Pearson correlation coefficient (r) for each year (Figure 12) and confirms that the absence of vegetation induces a high surface temperature compared to the vegetated land, which records a low temperature (Zheng et al., 2014). However, the negative Pearson correlation coefficient values found in the relationship between the LST and NDVI indicate that the presence of green areas mitigates the surface temperature effect because the vegetation is wet and contains water almost throughout most months of the year through the leaves compared to other types of LULC. Therefore, it moistens the surface and absorbs more and reflects less solar radiation (Davi et al., 2006; Morawitz et al., 2006).

In contrast, water surfaces have often recorded high temperatures (33.74°C in 2015 and 34.54°C in 2022 for the mean temperature). In fact, the waterbody class is not a large river or big stream. It comprises puddles, backwaters, small lakes, ponds, and swamps. As noticed, it is the class with the smallest extent and represents less than 1% of the total area of the classes in 2015 and exactly 0.34% in 2022. This means that water dries up, especially in the dry season, and becomes bareland. As a result, the temperatures will be very high, like in the bareland and savannah areas. Indeed, our study has shown that the hottest surfaces with the highest temperatures are the bareland/settlement and savannah. The mean temperatures of these LULC classes are 35.18°C and 35.41°C for bareland and settlement, respectively, in 2015 and 2022 and 34.60°C and 33.44°C for savannah, respectively, in 2015 and 2022. The bareland and settlements are devoid of vegetation, and the savannah is dominated by grasses and shrubs, which cannot strongly absorb sun rays. Therefore, reflections are strong, which justifies the high temperatures recorded at these LULC classes. A recent study (Moisa et al., 2022a; Moisa et al., 2022b; Moisa et al., 2022c) on the wettest parts of Ethiopia, covered in the present study, confirmed that the increasing trend of temperature is driven by weak NDVI values and LULC conversion. One of their outputs concerns the great variation between daytime and nighttime temperatures, which can influence the LST. As highlighted by the same authors, the LST increases during the daytime more than at night. Akomolafe and Rosazlina. (2022) found a moderate negative relationship between NDVI and LST in Penang Island, Peninsular Malaysia. This shows that higher LSTs correspond to lower vegetation cover and vice versa, as our results also show. Through the urban heat island characteristics and mitigation strategy analysis of eight arid and semi-arid gulf region cities, Abulibdeh (2021) stated that vegetation, particularly trees, can prevent direct surface heat as a result of solar radiation by providing shade, cooling the air by generating cool island effects, evapotranspiration and emissivity processes, and reducing the wind speed under the canopies. In the Pearson correlation coefficient analysis, the same author (Abulibdeh, 2021) stated that the LST has a significant relationship with NDVI but in different directions. The relationship between the LST and NDVI is negative, indicating that the presence of green areas mitigates the LST preponderance effect. All these claims by different authors confirm our results and demonstrate once again the importance of LULC changes and their effects on the LST modification.