A novel AI approach for modeling land surface temperature of Freetown, Sierra Leone, based on land-cover changes

ABSTRACT Land use/land cover (LULC) indices can be considered while developing land surface temperature (LST) models. The relationship between LST and LULC indices must be established to accurately estimate the impacts of LST changes. This study developed novel machine learning models for predicting LST using multispectral Landsat images data of Freetown city in Sierra-Leon. Artificial neural network (ANN) and gene expression programming (GEP) were employed to develop LST prediction models. Images of multispectral bands were obtained from Landsat 4-5 and 8 satellites to develop the proposed models. The extracted data of LULC indices, such as normal difference vegetation index (NDVI), normal difference built-up index (NDBI), urban index (UI), and normal difference water index (NDWI), were utilized as attributes to model LST. The results show that the root-mean-square error (RMSE) of the ANN and GEP models were 0.91oC and 1.08 oC, respectively. The GEP model was used to yield a relationship between LULC indices and LST in the form of a mathematical equation, which can be conveniently used to test new data regarding the thematic area. The sensitivity analysis revealed that UI is the most influential parameter followed by NDBI, NDVI, and NDWI towards contributing LST.


Introduction
The climate change encompassing global warming leads to forest fire, heavy rain, and land drought environmental disasters throughout the world (Pinto et al. 2020;Otgonbayar et al. 2019). With recent advancements in technology, remote sensing through satellite systems integrated with geographic information systems (GIS) is emerging as an advanced tool that can be employed to detect environmental changes. It carries numerous applications owing to its acceptable accuracy, the satellite range of spectral solution, frequent time, and the availability of long-term pledge of data continuity (Guzinski and Nieto 2019). Among those applications are geographical, including land use and land cover (LULC), and environmental, including temperature change, are worth mentioning Abdi 2020;Li et al. 2020). Yang, 2015). Tran et al. (2017) applied kernel ridge regression (KRR) to estimate the LST based on LULC indices, and the results presented high performance for this method in estimating LST values. Bartkowiak et al. (2019) applied a random forest (RF) method to enhance LST estimation from lowaccuracy satellite images. Their results showed that the proposed method improved the extracted LST. Huntingford et al. (2019) summarized the machine learning concepts for climate change studies. It was concluded that machine learning would expose climate system qualities and improve the predictions across time scales. Herein, advanced machine techniques, including artificial neural network (ANN) and gene expression programming (GEP), are employed and compared to extract a new model for LST prediction in the study area. These methods are widely applied to predict different engineering applications (Jalal et al. 2021a,b;Kaloop et al. 2022;El-Diasty 2019). Furthermore, ANN was used in Brazil based on positional variables, temperature, and average air relative humidity to estimate the LST and the accuracy of the proposed model was shown high, equaling R 2 of 0.93 (Veronez et al. 2013). Similarly, Wang et al. (2021) proposed an ANN model with a reliable accuracy for modeling the LST in China by using the sensor's thermal infrared (TIR) bands. Moreover, traditional ANN and deep neural network accuracy was shown close in modeling LST . The relevant studies employing ANN for predicting LST can be found in (Yoo et al. 2020;Tan et al. 2019;Mustafa et al. 2020a;Haghbin et al. 2021). On the other side, Genetic programming (GP) was applied to estimate the global temperature depicting a better performance than the linear approaches (Stanislawska et al. 2012). The prediction models for atmospheric temperature revealed the superior performance of GEP over ANN (Azamathulla et al. 2018). Although GEP is widely applied in different engineering applications related to climate change due to its unique advantage of furnishing simple mathematical non-linear equations (Stanislawska et al. 2012;Azamathulla et al. 2018;Shirani Faradonbeh et al. 2017), its applications in predicting LST are still limited. The relative literature about the subject matter suggests a non-linear relationship between LULC and LST, a desideratum to evaluate. Moreover, the impact of contributing parameters towards LST is necessary to investigate in order to plan the land use of the study area and future case studies.
Although the literature contributed extensive studies in LST prediction and studied the impact of LULC on LST, limited studies have been implemented to design a formula between them. Therefore, developing a formula for predicting LST from the LULC data extracted through satellite images presents the novelty of the current study. Herein, the main objectives of the present study are directed to: (1) find the non-linear relationship between LULC indices and LST, (2) estimate a formula for predicting LST that correlates different factors involved in LULC, (3) develop and compare ANN and GEP models as machine learning algorithms to estimate accurate LST values, and (4) analyze the sensitivity and parametric study to see the effect of contributing input attributes in modeling LST. For this purpose, Freetown city in Sierra-Leon was selected as a case study to assess the influence of LULC changes. According to Landsat images, the growth of this city over the selected time period, 1998, 2000, 2010, 2018, and 2019, has been significant. This city was also selected because the urban heat island (UHI) impact shows high on the air and land temperature changes . Therefore, the relation between LST and LULC indices was modeled and evaluated. The impact of LULC change on LST was assessed and predicted to provide feedback to policymakers and urban planners in developing LST mitigation strategies for similar areas.

Study area
This study was implemented in Freetown (capital city), a western coastal area in Sierra-Leon ( Figure 1). The development of this city affects land use and the air and surface temperature changes (Tarawally, Wenbo et al. 2018;USAID 2016;BTI 1990;Tarawally, Xu, et al. 2018). Freetown's position is on the Atlantic Ocean, and it is the major port city in Sierra-Leon. It lies between 08 o 05 ' N and 08 o 30 ' N latitude and 12 o 30 ' W and 13 o 20 ' W longitude. The city has a total area of 357 km 2 , and its population accounts for 15.53% of the entire Sierra-Leon population (c). Recently, Tarawally et al. ( , 2019, evaluated the urban growth of Freetown. It was found that the built-up showed a continuous increase in the city; in addition, the dense vegetation and agricultural land reduced by 3807 and 9145 ha, respectively, between 2000 and 2015. Furthermore, the average surface temperature of the city increased from 23.7°C to 25.5°C between 1988 and 2015 . The topography of Freetown fluctuates, ranging between 100 and 700 m of elevation change. The seasonal variation is that May to October is the rainy season, while November to April is the dry season of Freetown. The dry season was selected in this research to study the high impact of LST change in Freetown, and it was also the best period for collecting cloud-free optical satellite images. The annual average minimum and maximum temperatures of Freetown are about 23.8°C and 29.9°C, respectively . The prevailing winds are the southwest monsoon during the wet season and the northeastern harmattan, a dust-laden wind from the Sahara Desert during the dry season. The pink bordered areas ( Figure 1) are rapidly developing. In this study, the Landsat images were captured in the hot and dry seasons to evaluate the effectiveness of LULC changes on the LST.

Image processing and datasets
3.1.1. Image datasets The data for the selected periods 1998, 2000, 2010, 2018, and 2019 were collected from the Earth Resources Observation and Science (EROS) center through the United States Geological Survey (USGS) Global Visualization Viewer. Herein, cloud-free and geometrically corrected Landsat images were used. The path/row 202/54 was used with different bands to estimate an accurate LST and LULC, as presented in Table 1. The date of images and sensors used are shown in   Table 1. The selected period is designed to study the low and high variation of the data changes in the modeling process, considering the visibility and quality of available images. The satellite images were projected by the Universal Transverse Mercator (UTM) within Zone 29 North-Datum World Geodetic System-1984 (WGS84). GPS reference points for image adjustment and the data of population and temperature from Statistics and metrological Departments were obtained and used. High-resolution contemporary satellite imagery (GEOEYE-1 and Google Earth historical image of 2015), administrative spatial datasets from the National Tourist Board and Environmental Protection Agency, and ancillary secondary maps were also used as ground truth data for accuracy assessment. Various software packages were utilized because each one has strength in some operations needed for this study.
3.1.2. Image processing and LST estimation TM and OLI sensors were used with third-order polynomial fitting model to fit and nearest neighbor resampling techniques. Digital numbers (DN) of Landsat TM and OLI images were deposited as 8 and 16 bits, respectively. These numbers were transformed to the top of the atmosphere (TOA) spectral radiance (L), which can be calculated as follows (Larnicol et al., 2018;Tarawally et al. 2019): where, A and B = coefficients depend on the atmospheric and geometric circumstances but not on the surface; r and r e = the pixel and normal surfaces reflectance's, respectively, for the pixel and surrounding areas. S is the round albedo of the air, and L a is the radiance backscattered by the atmosphere.
In Eq. (1), to reduce the atmospheric effects, the band interleaved by line (BIL) format was obtained from the reflective bands' radiance. The FLAASH method is applied to provide the corrections of the surface reflectance imagery as follows (Tarawally et al. 2019;Larnicol et al., 2018): In order to calculate the LULC, the multi-spectral optical data were stacked into a multi-band layer required as input for classification. Table 2 presents the vegetation and urban indices used in the current study. Based on the previous studies (Mustafa et al. 2021a,b;Mustafa et al. 2020b;Tarawally et al. 2019Gbanie 2015;Tarawally, Xu, et al. 2018), these indices were observed to affect the LST changes highly.
In order to estimate LST, Weng et al. (2007) summarized the LST determination as follows: (i) transformation of DN of Landsat images to radiance, (ii) transformation of radiance to the temperature of blackbody, and then (iii) the blackbody temperature due to brightness transforms to LST (emissivity correction). These steps can be presented as follows: Eq. (3) presents the transformation process of DN into spectral radiation. The data of DN for the TIR bands of images is transformed. Eq. (4) is applied to obtain thermal infrared images of Landsat 8, based on the standard of USGS.
Where, L l denotes spectral radiance in W/(m 2 srmm) gained by the sensor found from each pixel from the Landsat imagery. Band-specific multiplicative is denoted by M L and rescaling factors which are gathered from MT L image file is found to be additive and is denoted using the term by A L . Q cal and QCAL max , respectively represent DN of each image and maximum DN. The value of Q cal and QCAL max for the 16-bit Landsat 8 is 65535 and for other Landsat missions is 255, respectively. QCAL min is the least DN (0). Radiances from top of the atmosphere (TOA) are denoted as L max and L min that are scaled to QCAL max and QCAL min in W/(m 2 srmm), respectively.
The radiant images can be converted to temperature detected in blackbody as follows: Here, the brightness temperature in the Kelvin unit is measured using sensor, which is denoted with T b , spectral radiance in W/(m 2 srmm) is denoted with L l and through image MT L file prelaunch calibration constants in Kelvin unit is obtained. These are represented as K 1 and K 2 . Error in the measurement of surface temperature can arise during the execution process due to detection of earth as a blackbody because of its brightness (Weng 2001). To minimize these error values, emissivity correction is performed to attain LST from T b using Eq. (4) (Nichol 1994).
Meanwhile, a threshold method was suggested by Sobrino et al. (1990) to achieve land surface emissivity from NDVI. If NDVI < 0.2 then pixels are acknowledged as barren land. The emissivity of this barren land was derived from red spectral region. If NDVI > 0.5, pixels are acknowledged as lands thoroughly covered with plants or agricultural fields. Their emissivity value was estimated to be 0.99 (Sobrino et al. 2004). When the NDVI value lies between the range of 0.2 and 0.5 then pixels are acknowledged as partially covered with vegetation. The surface emissivity (1) can be estimated as follows: Where, 1 v and 1 s represent the emissivity of vegetation coverage and soil surface, respectively, P v is computed to obtain the proportion of vegetation as follows: Where NDVI s and NDVI v are the NDVI values of pure soil and pure vegetation, respectively. D1 in Eq. (6) indicates distribution of the land surfaces according to its geometry in addition to their internal reflection. The value of internal reflection is very low for ordinary and uniform surfaces (Sobrino et al. 2004). The value of internal reflection is, however found to be 2% for rough and heterogeneous surface. D1 can be estimated as follows: Where shape factor is denoted by F. The mean value for different land surfaces according to its distributions based on its geometry is assumed 0.5.
From Eqs. (7) and (8), the estimation of emissivity can be given as follows: Where, LST can be estimated by minimizing the errors of surface emissivity as follows (Tarawally, Xu et al. 2018): Where, l is the wavelength of emitted radiance (11.5 µm) (Tarawally, Xu et al. 2018), r = hc/s (mK), K is the Stefan-Boltzmann's constant (1.38 × 10 −23 JK −1 ), h is the Planck's constant (6.26 × 10 −34 Js), c is the velocity of light (2.998 × 10 8 ms −1 ), and 1 is the surface emissivity. Finally, the derived LST values were converted to the conventional Degree Celsius (°C) unit by adding the absolute zero, which is approximately minus 273.5°C.

LULC and LST extractions
In order to model the LST of Freetown city, 400 points over the city map were randomly generated and selected. Appendix A presents the extracted maps of LULC indices and LST of Freetown. Figure 2 presents the DEM of the study area and selected points. The selected points were generated to cover the whole area of the city to estimate an accurate LST map in the prediction stage. To design an accurate prediction model, the LULC indices are considered in this study. The geographic effects of land commonly affect the LST (Khandelwal et al. 2018); in this study, the land elevation has been evaluated; however, it is not significantly affecting the LST in Freetown, as presented in the result section. Therefore, NDVI, NDWI, NDBI, and UI are used to model the LST. Here, the data from 1998 and 2000 were used to design the proposed models (training stage), 2010 and 2018 datasets were used for testing and validating the models, respectively. In addition, as presented in Figure 2, three sectors' of 2019 datasets were used to validate and evaluate the accuracy performance of the proposed models. Table 3 presents the statistical evaluation, mean, range, and standard deviation (SD), extracted LULC indices, and calculated LST.
From Table 3, it can be seen that there are significant changes in the LULC indices from 1998 to 2018. The range change of NDVI increased till 2010 and then decreased, as the same other indices. The SD change of the index is small. However, the LST change enormously increased from 23.43°C to 26.53°C between 1998 and 2018, respectively. This infers that the relation between LULC indices and LST is non-linear. Figure 3 and Appendix B demonstrate the linear correlation coefficient between LST and input variables to further understand the relationship. From the linear correlation results, it can be seen that the individual correlation between LST and input variables can be detected. A high positive and negative correlation is observed between LULC indices and LST, with a correlation between 0.59 and 0.89. This also indicates that the LULC indices variables significantly affect the LST values in Freetown. Previously, different studies were implemented based on one or two variables for predicting LST (Mustafa et al. 2021a(Mustafa et al. ,b, 2020b. However, the non-linear contribution of these variables on LST cannot be estimated through multi-linear regression models. Thus, this study aims to find the contribution of these variables to LST estimation.

Prediction models design and assessment
3.2.1. Gene expression programming (GEP) method GEP is an advanced and mature version of genetic algorithms (GAs) and genetic programming (GP) (Sharadini Faradonbeh et al. 2017). GAs consists of individuals having linear chromosomes of fixed length and encompasses simplified expressions. GP comprising evolutionary algorithms (EAs), is an automated modeling technique based on Darwin theory of the 'survival of the fittest'. It consists of branched structure chromosomes of variable sizes and shapes (Eldrandaly 2010;Jalal, Xu, Iqbal, Javed et al. 2021a,b). GEP uses the advantage of both the GAs and GP to produce linear chromosomes and efficient computer programs (CPs) (Sharadini Faradonbeh et al.   Figure 4 describes the step-wise procedure of the GEP process. Initially, input and output attributes are inserted into the GEP model. The input data is bifurcated into the training and validation data. The number of chromosomes, genes, and other properties (head size) are randomly fed into the computer program, which are expressed in the form of mathematical expressions and ETs. The fitness of the chromosome is evaluated in the sequential step according to the fitness function. Chromosomes are selected and replicated to the next generation based on the fitness functions. Frequently used fitness functions includes: the root mean square error (RMSE) , the sum of squared error (SSE) (Jalal, Xu, Iqbal, Jamhiri et al. 2021a), Correlation coefficient (R) , and mean absolute error (MAE). Furthermore, transposition, mutation, inversion, and recombination are applied to the chromosome to produce the next generation. Table 4 manifests the setting parameter for the developed GEP models. A varying number of chromosomes executed several trials, the number of genes and head size to achieve an accurate model. It can be seen that three genes, 50 chromosomes, and ten head sizes yielded an accurate model.

Artificial neural network (ANN) model
ANN is an empirical modeling technique used for prediction following the desired output (Köroğlu 2019). The motivation of the neural network originated from the human brain (Prakash et al. 2008;Aleboyeh et al. 2008). In this method, the relationship between the input and output is established by the number of hidden layers; hence the model is termed 'black box model' (Elmolla et al. 2010). The neural network typically consists of three or more layers, namely an input layer, one or more hidden layers, and an output layer. (Köroğlu 2019). A multi-layer feed-forward network with backpropagation is robust enough to solve complex engineering problems due to high non-  Onyelowe et al. 2021 ;). The error backpropagation (EBP) algorithms, also known as the steepest descent algorithms, were the breakthrough in training neural networks. EBP algorithms are used commonly, but they have an associated demerit of slow convergence. This drawback was addressed by Gauss-Newton algorithms using second-order derivatives of error functions, and, as a result, they were observed to converge faster, provided the quadratic error function. Gauss-Newtonian algorithms are primarily divergent. Kenneth Levenberg and Donald Marquardt blended EBP and Gauss-Newtonian algorithms as Levenberg-Marquardt (LM) algorithms which can switch between the EBP and Gauss-Newtonian algorithms and converge fast (Yu and Bogdan 2010). The LM backpropagation algorithm has been widely used and a preferable type for network training in the recent past (Suratgar et al. 2005;Mansouri and Kisi 2015) and is thus employed in the current study. Figure 5 shows the preliminary architecture of ANN used in this study. The neural network randomly categorizes the data imported from the database as training, testing, and validation sets. The training set is used to train the network by pairing a set of inputs with the respective expected output. The feed-forward pass of LMBP algorithms transmits  Figure 5. Architecture of ANN model. a signal from the input vector to the neurons in the hidden layer. The neurons in the hidden layer map the external inputs and expected outcomes in some meaningful way. The presence of multihidden layers uses the output signal of the first hidden layer as input for the second hidden layer, and the process propagates through the layers. During the forward pass, outputs are generated, and weights are fixed. The error between the network generated and the desired output is calculated and transmitted as an error signal during the backward pass. All weights are adjusted during the backward pass according to the error correction rule. The major limitation of the ANN method includes its black box nature and susceptibility to overfitting (Tu 1996). The limitation of overfitting was addressed by selecting the model with minimum objective function (OBF) as given in Eq. (12). Table 5 lists the setting parameters for the ANN model.
The subscript T and V represent the training and validation data, and n is the total number of sample points. r represents the performance indicator.

Models' performance assessments
The correlation between observed and predicted values is evaluated using Correlation coefficient (R) and variance account for (VAF) statistical indices to assess the proposed models' performances. Root mean square error (RMSE), mean and maximum absolute errors (MAE and XAE, respectively), and percentage of model error (PE) is calculated to assess the error of the proposed models in estimating LST. The assessment indices can be calculated as follows: where, d and y represent the measured and predicted LST values, respectively;d mean , d max , d min are the average, maximum, and minimum of measured LST values, and n is the number of observed datasets. s d and s y are the standard deviation of measured and predicted LST values.
Moveover, the sensitivity analysis (SA) of input parameters towards LST was determined as follows: where, f max (v i ) and f min (v i ) corresponds to the maximum and minimum values of i th input domains.
Herein, the interpolation method of Kriging in GIS is applied to estimate the prediction map of LST for the study area. Kriging is a multistep geostatistics approach that uses the variogram as a basic tool, then selects appropriate methods and performs optimal linear unbiased interpolation estimates on parameter space structure and randomness (Shi 2014;Tran et al. 2021). It posits that the distance between sample sites shows a spatial correlation that can be utilized to explain surface variation. It weights the surrounding measured values to make a prediction for an unmeasured location. When there is a spatially associated distance or directional bias in the data, Kriging is the best option. Also, it is frequently utilized in LST modeling, forecasting, and mapping (Yang 2004;Wang et al. 2021;Bhattacharjee et al., 2020). The Mathematical formula for Kriging that can be used for the prediction of unknown LST at location (S 0 ) can be formulated as: Where LST(S i ) represents the measured LST at (i th ) location, and l i is the weight of the measured value. Figure 6 presents this study's conceptual framework of data and models processing.

Results and discussion
According to Appendix A, the LST changed significantly from 1998 to 2019, as illustrated in Figure  A.1. The LST changes in coastal areas are highly increased due to the development changes (Tarawally et al. 2019. The expansion of the built-up area ( Figure A.4) into agriculture and sparse vegetation ( Figure A.2), and wetlands ( Figure A.3), affected the LST changes (Tarawally et al. 2019. In addition, a high correlation between LST and UI can be extracted from Figures A.1 and A .5; this agrees with the obtained results in  and Mustafa et al. (2021aMustafa et al. ( ,b, 2020b. Appendices A and B, as well as Figure 3, show that the LULC indices significantly impacted LST values. Tarawally et al. (2019 and Mustafa et al. (2021aMustafa et al. ( ,b, 2020b evaluated the LULC and LST of Freetown and designed regression models to predict LST. Furthermore, they validated the extracted LST and LULC in the study area. Here, more advanced methods have been used to estimate the LST of Freetown. This section presents the detailed performance and validation of the developed models, followed by 2nd level validation carried out in three distinct sections. The data from 1998 to 2018 were utilized for considering the development change of the overall area of Freetown city. 2019 datasets were used to assess the performance of the proposed model with different concepts of data distribution of the study area and check the accuracy of the designed models. The developed models were also validated in the third stage by evaluating the parametric and sensitivity contribution of the input attributes.

Models performance and 1 st level validation
During the training of the models, several trials were run to obtain optimized parameters for the ANN and the GEP models. The models were trained using the best hyperparameters parameters. Table 6 shows the statistical performance of the proposed models in the training, testing, and validation stages. The training data for the ANN model shows a minimum value of R equaling 0.87 for the test data. The training and validation data revealed R values of 0.88 and 0.91. For the GEP model, the values of R were observed as 0.84, 0.85, and 0.89 for the training, test, and validation data, respectively. The value of R obtained in the test and validation data is either comparable to or greater than the value of R for the training data. This means that the models were trained such that no overfitting issue was created during training models. The performance of the ANN model in training was high in correlation assessment indices. Also, the model error was 0.08°C in RMSE terms. The model can be used to estimate LST with a 7.81% error. The performance of ANN was shown to be high too in the testing and validation stages. The correlation and VAF values were average at 0.89 and 0.79, respectively. Previously the literature revealed that the values of R greater than 0.8 represent a close agreement of observed and forecasted results (Iqbal et al. 2021a,b;Onyelowe et al. 2021;Tran et al. 2021). The values of R for ANN model in the testing and validation phase were almost equal or greater than the training values, hence manifesting no overfitting of the developed ANN model (Iqbal et al. 2021b). In addition, the average model RMSE, MAE, and XAE were 0.91, 0.77, and 4.47°C, respectively. Therefore, the ANN model could be used to estimate LST with an 8.6% error, on average. In comparison, the average performance of GEP in terms of R and VAF were 0.86 and 0.69, respectively. The performance of average RMSE, MAE, and XAE for the GEP model were 1.08, 0.95, and 4.59°C, respectively. Furthermore, the GEP model could predict the LST with the percentage of error approaching 10.25%. The magnitudes of correlation obtained for the testing and validation phase also avoided the issue of the overfitting model in the GEP model. Therefore, these results indicate that both models can be used in modeling the LST of Freetown. However, the ANN outperformed GEP in estimating LST values. Figure 7 presents the estimated maps of GEP, and ANN models compared to the observed LST from satellite images.
From Figure 7, it can be seen that there were small errors that can be detected. The GEP model's performance showed accuracy in modeling 1998 datasets, while the ANN model outperformed GEP in modeling LST for 2000, 2010, and 2018 datasets. In general, both models can be used to estimate the LST; however, ANN was superior to estimating the accurate LST for Freetown. Herein, the accuracy of ANN and GEP models was assessed at 95% confidence level. For the testing and validation stages, the empirical cumulative distribution function (CDF) of model errors was computed and is shown in Figure 8. It is shown that the majority of error values at 95% confidence of ANN model were within ±1.70°C and ±1.46°C for the testing and validation stages, respectively. For the GEP model, most model errors at 95% confidence interval were within ±2.08°C and ±2.10°C, respectively. Thus, ANN was better than GEP in modeling the LST of Freetown.
The equation was extracted from the developed MATLAB model of GEP, and the expression trees were generated as depicted in Figure 9 (Iqbal, Zhao, et al. 2021;Jalal, Xu, Iqbal, Javed, et al. 2021b). The corresponding variables d0, d1, d2, and d3 in Figure 9 are input attributes, namely NDVI, NDWI, NDBI, and UI, respectively. The constant 'c8' in Sub-ET1 was 1.76; c0, c7, c1, c4, and c3 in sub-ET2 were 1. 44, 1.15, 7.61, 2.17, and −9.23, respectively; and the constants in sub-ET3 were c5 and c7 equaling 7.63 and −7.85, respectively. From Eq. (23_ and Figure 10, it can be seen that the pattern impact of UI index significantly affected the LST change, followed by the NDWI index. Meanwhile, the pattern impact of NDBI and NDVI indices was respectively small. Herein, the geographic impact was evaluated using the GEP model. Two scenarios were considered by using an unused geographic impact, elevation height, in GEP modeling of LST. Considering the elevation height, the GEP model accuracy was not changed. Therefore, the effect of geographic change could be neglected in modeling the LST of Freetown, as presented in Figures 9  and 10.

2nd level validation
The three sections see Figure 2; datasets were applied in developing ANN and GEP models to validate the proposed models further. Different changes of LULC indices at different elevation heights were considered to assess the model's accuracy in the different design patterns. In addition, the data number effect was evaluated for each section. Table 7 presents the statistical analysis performance of the three sections. The predicted LST was calculated using the developed structure of ANN and Eq. 23 of the GEP model. Figure 11 illustrates the observed and predicted LST values for the three sections. From Table 7 and Figure 11, it can be seen that the correlation between observed and predicted LST was high for the three sections of both models; however, the ANN performance was overall better and more robust. The data number did not affect the performance of both models; the percentage of model's errors was high for sec-3, which had a data number greater than sec-2.
In conclusion, the ANN model was robust to use in predicting the LST of Freetown, and the GEP model could be used as a numerical model to estimate the future performance of LST values based on the range of input variables. Herein, Eq. 23 could be considered the non-linear relationship between LULC indices and LST for Freetown city.

3rd level validation
Due to the overfitting problem, we validated the ANN model generated in this study on a completely different data set while formulating AI models. For this purpose, a simulated dataset was created, as shown in Table 8. The effect of changing parameters was investigated while keeping the remaining variables constant. The results of parametric and sensitivity analyses are presented in Figures 12 and 13, respectively.
It has been investigated that an increase in NDVI and NDWI yielded a corresponding decrease in LST (Marković et al. 2021;Chen et al. 2013); hence the current parametric study (Figure 12) captured the same trend of LST with varying NDVI and NDWI. Keeping other impacting parameters constant, change in NDVI from −0.37 to 0.20 decreased LST by 3.68°C. Similarly, an increase in the NDWI decreased LST by 1.51°C. In contrast, the rise in NDBI yielded an increasing trend of LST. An increase of 3.9°C was observed for NDBI variation (−0.19-0.32). Finally, LST was observed to have a complex 3-degree polynomial variation with a rise in UI.
The sensitivity analysis reveals how sensitive a developed model was to a particular change in the parameters considered for the respective model since the parameters are ranked as per their relevance (Jalal, Xu, Iqbal, Javed et al. 2021b). The developed ANN model was used to test new unseen  data sets to observe the effect of changing parameter values in its range. The contribution of input parameters toward LST is presented in Figure 13. Here, the remaining input parameters were maintained at average values when determining the maximum and minimum values of i th input domains. It can be seen from Figure 13, that UI was the most influential parameter contributing 30.16%, followed by NDBI (27.5%), NDVI (24.73%), and NDWI (18.04%) towards LST, respectively. Therefore, it can be concluded that the urban changes (UI and NDBI) can be significantly used to evaluate and estimate the LST change, while the spectral vegetation indices (NDVI and NDWI), or changes in green spaces, cannot only use for analyzing the change of LST.

Conclusions
Land surface temperature (LST) significantly depends on land use and land cover (LULC) parameters. The mitigation of LST and proper urban planning require an accurate estimation of LST. This study investigated the relationship between LST and LULC. Two novel machine learning prediction models, namely ANN and GEP, are employed to develop LST prediction models for Figure 11. 2019 evaluation of sections 1,2, and 3 (upper to lower rows, respectively) for the GEP (upper) and ANN (lower) models. Freetown city, Sierra Leone. Moreover, the interpolation method of Kriging in GIS is applied to generate the prediction maps of LST. Following conclusions can be drawn from this study.
. The LULC indices achieved from the study area for 1998-2018 were used to train non-linear ANN and GEP models with the best hyperparameters obtained using trial and error analysis. The Landsat 4/5 and 8 images were used to collect the input and output variables. The evaluation of the proposed models showed that the developed ANN and GEP models could be used to estimate LST with 8.6% and 10.25%, on average, respectively. In addition, a close agreement is observed between the estimated and observed LST values for both models. The uncertainty models showed that most error values at 95% confidence of ANN and GEP models fall within ±1.58°C and ±2.09°C, respectively. . The developed models were validated by testing the entire three new datasets of LULC and LST, taken from three different sections of the study area in 2019, on the basis of trained GEP and ANN models. The validation data revealed an average MAE of 1.15 and 1.37°C for three sections of the ANN and GEP model, respectively. Besides, the average correlations equaling 0.833 and  0.78 were obtained for the ANN and GEP models. This analysis revealed that the ANN model is more robust than the GEP model. For this reason, third-level validation was conducted for ANN model in the form of parametric and sensitivity analysis. . The sensitivity of ANN model reflected the existence of the highest contribution of UI towards LST. ANN model suggests that NDVI, NDWI, and NDBI have appreciable contributions to the LST of Freetown; hence, they shall be considered by the town planners in the land use process. Meanwhile, the NDVI and NDWI indices cannot only be used to analyze the LST change. The parametric contribution of UI, NWDI, NDBI, and NDVI showed adherence to the literature, thus further validating the developed ANN model. The parametric study also supported the results obtained from the individual linear correlation depicting a high negative and positive coefficient between NDVI ∼ LST and UI ∼ LST, respectively. . The geographic impact i.e. the effect of datum value on LST was evaluated using the GEP model, and the results showed that the effect can be neglected in modeling the LST of Freetown. The evaluation of the developed models at different sectors with changes in LULC indices showed that the ANN model is more robust to predict Freetown's LST. However, the GEP model, can also be used as a numerical model to estimate the future performance of LST values of Freetown based on the range of input variables.

Data availability
The authors confirm that the data supporting the findings of this study are available within the article.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This work was supported by Korea Agency for Infrastructure Technology Advancement: [Grant Number 21CFRP-C163381-01].