Application of Spatial Modeling of Biophysical Variables in an Urbanized Area in the Amazon: The Case of the Metropolitan Area of Belém-Pará

The type of land use and land cover plays a decisive role in land surface temperature (LST). As cities are composed of varied covers, including vegetation, built-up areas, buildings, roads and areas without vegetation, understanding LST patterns in complex urban spaces is becoming increasingly important. The present study investigated the relationship between LST and albedo, NDVI, NDWI, NDBI and NDBaI in the period between 1994 and 2017. Images from Thematic Mapper (TM) and Thermal Infrared Sensor (TIRS) onboard the Landsat 5 and 8 satellites, respectively, were used in the study. The images were processed, resampled (spatial resolution of 120 m) in the environment of the QGIS 3.0 software and, finally, centroids were extracted resulting in a total of 1252 points. A classical regression (CR) model was applied to the variables, followed by spatial autoregressive (SARM) and spatial error (SEM) models, and the results were compared using accuracy indices. The results showed that the highest correlation coefficient was found between albedo and NDBaI (r = 0.88). The relationship between albedo and LST (r = 0.7) was also positive and significant at р < 0.05. The global Moran's I index showed spatial dependence and non-stationarity of the LST (I = 0.44). The SEM presented the best accuracy metrics (AIC = 3307.15 and R2 = 0.65) for the metropolitan region of Belém, explaining considerably more variations in the relationship between explanatory factors and LST when compared to conventional CR models.


Introduction
In several places around the world, cities are the main areas for the development of human activities and interactions. Such environments have faced extensive changes in land use and land cover (LULC) (Li et al., 2011;Silva et al., 2014;Silva et al., 2015;Silva et al., 2016;Silva et al., 2018;de Oliveira Serrão et al., 2020;Margalho et al., 2020) directly associated with population growth and economic development. The transformation of natural landscapes into human settlements has increased significantly over the past few decades (Grimm et al., 2000;de Oliveira Serrão et al., 2020), resulting in manmade environments as large urban centers.
Changes generated in the environment due to the expansion of urban areas represent the main indicator of increased air and surface temperatures, which consequently form urban heat islands (UHI). Li et al. (2018) emphasize that LULC are directly associated with the disorganized growth of urban centers, especially in developing countries. This phenomenon is generated because the materials used in the construction of urbanized areas absorb and retain great amounts of heat. Thus, UHI consist of an atmospheric phenomenon observed in urban environments, caused by anthropogenic action in the change and use of the soil. Such geographical changes caused by man have as one of their consequences the increase in air and surface temperature (Guha et al., 2018;Li et al., 2018).
Urban heat islands can be investigated at various scales. Microscale methods (∼ 30 m) consider localized urban spaces and other meteorological factors. In turn, mesoscale methods (∼ 20 km) are generally applied in large geographic areas to assess the spatial distribution of UHI in terms of Earth's surface temperature (LST) and thus to see how the thermal gradient of the surface can impact the dynamics of tubular flows, and result in possible changes in wind circulation (intensity and direction) motivated by the presence of buildings. Still and according to (Jan et al. 2018;Pacifici et al. 2019) atmospheric stability can be altered which is closely related to mesoscale weather conditions. This is possible through the formation of satellite image mosaics for large urban centers, as noted by Wanderley et al., (2019) in Brazil;Arulbalaji et al., (2020) in India and Choudhury et al., (2019) in Egypt. Land surface temperature maps are created from radiation data obtained by orbital remote sensing (RS). Most studies have applied RS techniques and used several satellite images, such as MODIS, TM -Landsat 5, ETM + -Landsat 7 and OLI/ TIRS -Landsat 8, to provide LST data and compare them with other factors, such as sealed surfaces (Li et al., 2011), vegetation cover (Yuan and Bauer, 2007;Duncan et al., 2019), and composition and configuration of various LULC types (Chen et al., 2006;Li et al., 2011).
The determination of vegetation indices from multispectral satellite sensors and RS techniques help to detect changes in the landscape at multiple scales. Spectral indices, such as the Normalized Difference Vegetation index (NDVI), detect vegetated (values > 0) and non-vegetated (values < 0) surfaces (Rouse et al., 1974;Guha et al., 2018). NDVI allows monitoring the health of the vegetation (Huete et al., 2002) and can also be used to delimit wet and dry periods , diagnose LULC (Fu and Weng, 2016;Guha et al., 2018), and evaluate changes in spectral patterns of preserved environments (Mancino et al., 2014). For the analysis of urban development based on built-up areas, several researchers (Varhney, 2013, García andPérez, 2016;Guha et al., 2018) have used the Normalized Difference Built-up Index (NDBI), which clearly distinguishes urbanized areas from gardens and xerophilic vegetation across built-up urban environments.
Several studies, among them Di Leo et al. (2016) and Luo and Peng (2016), have discussed the relationship between LST and biophysical variables. Some of these studies including Zhou and Wang (2011) and Zhao et al. (2018) explored such relationships using geographically weighted regression (GWR) models. Many used one or more methods where biophysical variables are the independent variables, to examine their relationship with LST. Fotheringham et al. (2015) stated that GWR allows correlating the response variables to a set of independent variables and is a statistical method for spatial modeling of heterogeneous processes. It is an approach for heterogeneous modeling of spatially distributed processes, and due to its greater analytical capacity and accuracy, it provides estimates of high precision and efficiency (Ahamdi et al., 2018). A fundamental component of GWR is the spatial weight by which the spatial relationships are constructed. Usually spatial weights are defined by spatial nuclear functions such as Gaussian function (Fotheringham et al., 2015), in which weights are related to the closer observations. Therefore, the objectives of this study follow two lines: (1) Investigate the relationships between LST, Albedo, NDVI NDWI, NDBI and NDBaI and their spatial characteristics, (2) quantify the patterns of distribution and spatial dependence LST through the Moran I index for Belém-PA, seeking to verify how vegetated areas on urban soil they can lower the ambient temperature and thus alleviate the effects of UHI.

Study area
The object of study is the urban area of Belém Metropolis (Fig. 1) located in the state of Pará, North region of the Brazilian territory, between latitudes 1°28'43.1"-1°16'40.3" S and longitudes 48°30'24.4"-48°22'28.8" W. According to IBGE, its extension is approximately 1,060 km 2 , and its population was estimated in 1,452,275 inhabitants and demographic density in 1,315.26 inhab./km 2 in 2017. The present study focused on the urban area of Belém, covering 184 km 2 . The rural areas and river islands within the territory of Belém metropolis were excluded from this research.

Data
Products from the Landsat-5 (TM sensor) and Landsat-8 (OLI/TIRS sensors) purchased from the National Institute for Space Research (INPE) for the years 1994, 2008 and 2017 were used in this study, these dates were chosen due to the availability of images. The images obtained were those of the months of June and July, when conditions of cloudiness were favorable, since the occurrence of clear skies is more frequent at this time of the year. The information and specifications of each image are shown in Table 1.
Spectral bands were resampled to a spatial resolution of 120 meters. Then, representative centroids of the census sectors were extracted, resulting in 1,252 points distribu-ted in the study area, in which each sample represents the pixel value corresponding to the (raster) layer of the image. All processing was performed within a GIS environment in the QGIS 3.0 software. For statistical analysis, scripts were used in the R environment (Development Core Team, 2018) and modules thereof with specific functions for analysis.

Data processing
2.3.1.  In image processing, radiometric calibration is processed by converting the gray level into spectral radiance of each band (L λi ), as proposed by Markham & Baker where: L λi (W m −2 sr −1 µm −1 ) is the spectral radiance of each of the bands; a and b are the minimum and maximum spectral radiance; ND is the pixel intensity, which ranges from 0 to 255; and i corresponds to the Landsat 5 bands.
The reflectance or reflectivity of each band (ρ λi ) is defined by the ratio between the reflected and the incident solar radiation flux, which, according to Allen et al. (2002), is given by Eq. (2): where: L λ is the spectral radiance of each band; K λ is the spectral solar irradiance of each band at the top of the atmosphere (Wm −2 µm −1 ); Z is the solar zenith angle; and d r is the square of the ratio between the Earth-Sun distance on the day when each image was obtained and the mean Earth-Sun distance. For determining the surface albedo Eq. (3), it is necessary to obtain the albedo without performing atmospheric correction, that is, the planetary albedo (α toa ), which is obtained from the monochromatic reflectances of bands 1 to 7, except band 6, in the case of TM -Landsat-5, as follows: where the values of each constant (weights) are determined by the ratio between the irradiance of each band and the sum of the irradiances. Then, the surface albedo is obtained, as in Tasumi et al. (2008), according to Eq. (4): where α atm is the reflectance of the atmosphere itself in the entire domain of solar radiation and which, and according to Bastiaanssen et al. (2000), the value 0.03 can be used; τ sw is the atmospheric transmissivity for clear skies, which according to Allen et al. (2002) can be estimated by the Eq. (5): where: z represents the altitude (m) of each pixel.

Landsat-8
To obtain the reflectance of the spectral bands of the OLI -Landsat-8 sensor, the digital number (DN) only needs to be converted into reflectance according to Eq. (6): where ρ λ,b is the reflectance of band b; H ρ is the specific multiplicative scaling of each band (2 −5 , constant); A ρ is the additive factor (−0.1); Q CAL is the digital value or digital number (DN); and θ SE is the solar azimuth angle estimated by the Eq. (7): where θ SZ is the sun elevation angle, available in the metadata of each image.
To determine the surface albedo, it is necessary to calculate the planetary albedo (α toa ) by combining the reflectances of each band (ρ λb ) with their respective weights (ϖ λ,b ), according to Eq. (8): where each weight (ϖ b ) is calculated from the ratio between the specific solar constant of each band (ESUN λ, b ) and the sum of all radiation constants ( P ESUN λ,b ). Therefore, the surface albedo or corrected albedo (α sup ) is calculated using the same expressions abovementioned for Landsat-5 (Eq. (4)).

Estimated land surface temperature (LST)
Images from the thermal band of TM -Landsat-5 (band 6) and TIRS-Landsat-8 (band 10) were used to estimate LST by Eq. (9). The processes and equations to estimate the surface temperature are similar and they will be described below. In both Landsat-5 and Landsat-8, the first step is to determine the spectral radiance of the thermal band. While for Landsat-5, the equation used to determine the radiance has already been described (Eq. (1)), for Landsat-8 the equation used to obtain the radiance of the thermal band is: where: L λ10 is the radiance of the thermal band; M L is the band-specific multiplicative rescaling factor, equal to 3,342 × 10 −4 Wm −2 sr −1 µm −1 ; Q cal is the pixel value; and A L is the band-specific additive rescaling factor, equal to 0.1. The second step is to calculate the emissivity (ε NB ). Then, after estimating the values of these two processes, the surface temperature can be estimated. However, to estimate emissivity, it was decided to adopt the leaf area index (LAI) as independent variable. To calculate LAI, it is necessary to obtain the soil-adjusted vegetation index (SAVI) which aims to mitigate the effects of the soil through Eq. (10): where L is a soil adjustment factor, varying between 0.25 in case of presence of dense vegetation, 0.5 in case of intermediate vegetation density, and 1 in case of low vegetation density (Huete, 1988). The value used in this study was 0.5. LAI is obtained by the ratio between the total area of leaves contained in the vegetation or in a given pixel and the area used by that vegetation or the area of the pixel. This index was calculated because it is present in the expression of emissivity. LAI is calculated according to Eq. (11): Emissivity (ε NB ) is introduced in the spectrum of the thermal band because each pixel does not emit electromagnetic radiation like a black body. According to Allen et al. (2002), emissivity can be obtained, in case NDVI < 0 and LAI < 3, by the Eq. (12): where ε NB is the emissivity in the domain of the thermal band; LAI is the leaf area index. For pixels with LAI ≥ 3, ε NB = 0.98; and for urban areas with NDVI < 0.2, ε NB = 0.92 (Nichol, 2009). LST was obtained using the Eq. (13): where: the K 1 and K 2 constants of the thermal band of Landsat-5 are 607.76 Wm −2 sr −1 µm −1 and 1260.56 K, and of Landsat-8 are 774.89 Wm −2 sr −1 µm −1 and 1321.08 K, respectively.

Determination of biophysical parameters
In the study, some indices were investigated for their possible relationship with surface temperature. A methodology validated in several studies was used to obtain these indices. NDVI (Rouse et al., 1973) and NDBI (Zha et al., 2003) are used to characterize land use and study the relationship between land use and surface temperature. NDVI, estimated from (Eq. (14)), stands for normalized difference vegetation index and is used to express the vegetation density (Purevdorj et al., 1998). The values range from −1 to +1; soils with vegetation cover present values between 0 to 1, and water bodies (lakes, river beds, etc.) present negative values (< 0).
In turn, NDBI, expressed by (Eq. (15)), was adopted in the present study for it is sensitive to built-up areas (Zha et al., 2003). For this index, the range varies from −1 to +1; positive values represent built-up and/or waterproofed areas and negative values, areas with vegetation cover.
According to Zhao and Chen (2005), NDBI allows the differentiation between primary and secondary bare land. This index is based on significant differences in the spectral signature in the near infrared between uncovered soil and the background. In turn, NDBaI is calculated from (Eq. (16)) as follows (Zhao and Chen, 2005): where d represents the digital numbers (DN) of the Landsat-TM for each band, IVP (medium infrared) and IVT (thermal infrared). NDBaI was proposed to analyze the spectral characteristics of different classes of land cover. Fully exposed soils (beach, bare land, developing land, among others) are distinguished in the images with NDBaI > 0, which also refers to primary bare soil. According to Chen et al. (2006), these indices can be used to classify different classes of land cover (vegetation, built area, among others). Index values vary for different types of land cover, and they are not constant.

Statistical analysis
Values of the spectral indices used in the spatial modeling are representative for each centroid of the census sectors in the urban area of Belém-PA and correspond to the average of the pixels of the three images used in the present research (1994, 2008 and 2017). All statistical analysis was done in the R software through the GeoR Package.

Moran I index (I)
The global Moran's I index (I) is a statistical mechanism for verifying the spatial dependence of a given variable and is one of the most used for this purpose. It is estimated by the following expression, as in Almeida et al. (2009) Eq. (17): where n is the number of observations; w ij is the element in the neighborhood matrix for the pair i and j; W is the sum of the matrix weights; y i and y j are deviations from the mean; y is the mean. This index measures the spatial autocorrelation based on the product of deviations from the mean. Thus, the Moran's I index is a global measure of spatial autocorrelation, as it indicates the degree of spatial association present in a data set. It varies in an interval (−1, 1). In case of spatial independence, its value is zero (0); positive values (between 0 and +1) indicate direct correlation, that is, perfect association with spatial dependence; and negative values (between 0 and −1) indicate an inverse correlation, that is, perfect dispersion.

Classical regression (CR)
A regression model (Eq. (18)) is based on the interest in evaluating the relationship of a variable (Y) with independent or covariate variables (X), that is, the relationship between two or more variables so that one of them can be explained or predicted by other variables (Corrar et al., 2007). In the case of spatial data, when there is spatial autocorrelation, the generated model must incorporate the spatial structure because the dependence between the observations affects the explanatory capacity of the model (Câmara et al., 2004).
where Y i is an observation of the dependent variable; X 1 , X 2 ,…, X n are the independent variables; β = (β 0 , β 1 , β 2 ,…, β n ) are referred to as corresponding regression coefficients, and ε i is the error associated with the observations of the dependent variable.
The assumption that the observations are independent simplifies the model, but in the context of data of area, this simplification is unlikely because of the possibility of spatial dependence between the error terms. An alternative is to use a mixed spatial autoregressive model (spatial lag model), which attributes the unknown spatial autocorrelation to the response variable y i . Another alternative is to apply a spatial error model, in which spatial effects are considered as noise, that is, as a factor to be removed (Fischer and Lesage, 2010).

Spatial autoregressive model (SARM)
The spatial autoregressive model allows the observations of the dependent variable y i in area i (i = 1, .., n) to depend on observations in neighboring areas with j ≠ i (Câmara et al., 2004), taking the form (Eq. (19)): where ε i is the error, W ij is the (i, j)-th element of the spatial matrix of n order (that is, n by n); the scalar ρ is a parameter (to be estimated) that will determine the intensity of the spatial autoregressive relationship between y i and P j W ij y j , whose interpretation corresponds to the average effect of the dependent variable in relation to the spatial neighborhood in the region analyzed; the Wy vector is known as space lag, the X matrix contains the observations of the independent variables, and the β vector has coefficients for the independent variables.

Spatial error model (SEM)
The spatial error model occurs when the spatial dependence is obtained through the error process, in which the errors of the different areas can present spatial covariance (Bivand and Piras, 2015). It is determined by the expression Eq. (20): where ρ is the auto-regressive parameter that indicates the intensity of the spatial autocorrelation between the residuals of the observed equation, and which measures the average effect of the neighbors' errors in relation to the residual of the region in question; and u i is the random error term, typically assumed to be i.i.d. The spatial autocorrelation in these models appears in the error terms.

Model selection
According to Dobson and Barnett (2011) the adjustment algorithm must be applied not only to one model, but to several models in a very wide set that must be relevant to the nature of the observations to be analyzed. If the process is applied to a single model, without taking into account possible alternative models, there is a risk of not obtaining the more adequate models to the data. To the usual regression models, the Akaike criterion (AIC) must help finding a sub-model for which the following quantity is minimized Eq. (21): where D y;μ ð Þ is the distance between the log-likelihood function of the saturated model (q parameters) and of the model under investigation (p parameters) evaluated in the maximum likelihood estimateβ. A small value for the deviation function indicates that, for a smaller number of parameters, the adjustment is as good as that obtained with the saturated model.

Correlation between biophysical parameters
The correlation coefficients (upper diagonal correlogram) between LST and biophysical variables are shown as a dispersion matrix (lower diagonal correlogram) in Fig. 2 and the three images used in the present study were constructed from the average of the pixel values (1994, 2008 and 2017). The highest correlation coefficient was observed between albedo and NDBaI (r = 0.89). This is due to the increased absorption of sunlight by elements that characterize buildings and urbanization. It is observed that the albedo stands out with a higher density of points in the range of 0.13 to 0.16. The relationship between albedo and LST (r = 0.67) was also positive and significant (р < 0.05). The relationship between albedo, NDVI and NDWI was negative because of the lower reflectance rate of vegetation and water bodies; these surfaces absorb more energy and use it in the evapotranspiration process.
According to Liu and Zhang (2011), vegetation has a great influence on the reduction of surface temperature. Thus, there is a negative relationship between the accumulation of green areas and UHI. On the other hand, thermal inertia is higher at impermeable levels, resulting in increased LST (Zhang et al., 2017). There was also a negative relationship between NDBaI and NDWI (r = −0.84) which possibly resulted from water forces characterized by NDWI values > 0, which indicate the presence of moisture and surface water in the urban perimeter of Belém-PA.
According to the histogram (Fig. 2), NDVI presented a minimum value of −0.35 and a maximum of 0.81, corresponding to water bodies and very vigorous vegetation, respectively. In contrast, LST presented an average value of 23.9°C for the period 1994-2017, with a minimum of 19.9°C and maximum of 27.7°C.
In studies such as those by Chen et al. (2006), Zeng et al. (2010, Zhang (2011) andRanagalage et al. (2017), the relationships between LST and biophysical variables showed that the surface temperature is positively correlated with indices for detection of built-up areas (NDBaI) and negatively with vegetation indices (NDVI).

Moran's I Index
In order to verify whether LST is dependent on space, it was necessary to apply the statistical index called Moran's I index, whose aim is to ascertain whether a particular variable has fluctuations in space. The scatterplot of the Moran's I index (Fig. 3) showed that the areas located in the upper right corner (first quadrant) showed positive spatial autocorrelation, that is, these areas formed clusters of similar values between neighbors. areas. The areas in the first quadrant that stand out in the graph are worrisome, because in addition to presenting an intense elevation of LST (High-High), they are also surrounded by areas with the same characteristic, that is, hot areas with warm surroundings. And they are represented by neighborhoods that have a high density of construction.
Several researches have used the Moran's I index (Deilami et al., 2018;Mavrakou et al., 2018 andYang andSantamouris, 2018) to analyze the adverse effects of UHI and under the perspective of appropiate urban development. It is noteworthy that the main object of the technique is based on the spatial heterogeneity of the distribution of LST (for example, high or low clustering pattern and dispersed spatial distribution).
The line that passes through the first and second quadrant corresponds to the slope of the simple regression; this value is equal to the global Moran's I index (0.44), showing that the variable has a positive and direct relationship with spatiality, that is, LST is dependent on space. Coutts et al. (2016) proposed the readjustment of green elements (urban reforestation) to reduce thermal hotspots. The study by Zhao et al. (2018) showed that spatial patterns of LST are statistically significant in hotspot zones in the center of the study area and partially extended to western and southern industrial areas, which indicated that the intensity of UHI had a notable spatial cluster in the city.
In the scatterplot of Moran's I index, it is observed that LST has spatial variability. Thus, in areas of the urban mesh of Belém (Fig. 4), preserved green zones are repre-sented by quadrant 2 (Q2) in the scatterplot of Moran's I index, presenting low LST values and surrounded by other low LST values (Low-Low). These vegetation zones are located in areas of an ecological park (Utinga Park), military areas in which the international airport of Belém is located, and areas still in the process of urban expansion, with few buildings, no asphalt or vertical constructions, and with the presence of vegetation.
The delimitation of spatial clusters becomes a strategic tool in the urban zoning of large cities. Furthermore, such analyses are extremely important, as highlighted by Wang and Hu (2012), due to their possible applications in environmental health (Ho et al., 2015), heat risk (Liu et al., 2016), water balance analysis (Tyralis et al., 2017), and energy generation.
In areas characterized by the quadrant 1 (Q1), it is possible to see higher LST values when the neighborhood has high LST values (High-High), demonstrating that areas with higher density of buildings and pavements present higher LST values. In these areas are some of the most populous neighborhoods in the city of Belém as well as the areas with lower vegetation density, and they also represent the downtown area. Yet, another quadrant (Q4) can be seen in this central region, represented by areas with lower ST values and neighboring areas with higher ST values (Low-High), representing zones that still have green areas (squares and forests), less verticalization compared to Q1 and greater presence of houses.
By means of a geospatial analysis of the vulnerability of formation of UHI in the city of Philadelphia, USA, Barron et al. (2018) demonstrated that areas of the city that had more urban characteristics, that is, areas with soil sealing and no vegetation, presented a greater risk of heat absorption and storage, causing greater vulnerability to formation of UHI. The authors also observed that regions with greater urbanization presented clusters (High -High) for streets with greater vulnerability for the formation of UHI. This study carried out in Philadelphia ratifies the findings of the present research, in which it is possible to see that regions of greater urbanization, greater population density, and decreased vegetation cover form clusters that facilitate the formation of an urban heat island due to increased heat storage.

Spatial model
The spatial variation of LST was related to the spectral indexes albedo, NDVI, NDWI and NDBI. Table 2 summarizes the accuracy of the spatial regression models used in the present study. The models generated for analysis were the classical regression model (CR), spatial autoregressive model (SARM), and spatial error model (SEM) ( Table 2).
The latter being the combination of the first two. According to Table 2, it was possible to verify that the SEM was better adjusted, as it presented a higher coeffi- cient of determination (R 2 = 0.65) and a lower Akaike criterion (AIC) value equal to 3307.15. In addition, it is observed that the significance of covariables in all models is adjusted for the analysis of the UHI in the city of Belém.
The maximum R 2 for the adjusted values between the dependent variable (LST) and the independent variables (biophysical parameters) were, respectively, 0.43 and 0.51, for the CR and ARSM models, respectively. Accord- ing to the spatial analysis, a spatial non-stationarity pattern (ρ > 0) was detected in the distribution of LST as well as in its associations with the explanatory variables.
In the result of the residuals of each model (Fig. 5), it was observed that the SEM presented the best fit, mainly in areas of the urban center of Belém-PA, with the highest urban density. For all models (CR, SARM, SEM), the transition areas from vegetation to buildings showed an overestimation of LST values, while underestimations were observed in the vegetated areas. However, both overestimation and underestimation values are relatively low.
In general, regression analyses revealed that the composition of LULC and the morphology of the terrain are closely related to the effects of UHI for the vast majority of the pixels used in the study. Additionally, the composition of the biophysical indices proved to be an essential factor directly influencing the LST pattern of the urban area of Belém-PA. The present research confirmed previous findings that an increase in building density (NDBI) tends to exacerbate the UHI effect, while an increase in the intensity of vegetation cover tends to mitigate such effect on the microclimate (Guo et al., 2015;Li et al., 2017;Stock et al., 2017).
As a natural process in large urban centers, the increase in LST exhibits high spatial heterogeneity that is difficult to be characterized with conventional regression methods. However, most of the previous studies derived the spatial relationship, focusing especially on cities in Asia (Yue et al., 2007). Such studies still highlight that the degree of relationship between LST and biophysical variables may be associated with climatic characteristics in a region where LST is less affected by land cover. Strong  associations between LST and biophysical variables are likely to be found in cities with more homogeneous characteristics of the land.

Conclusions
The results of the present spatial distribution analysis based on the global Moran's I index indicated the rejection of the null hypothesis, due to the observation of spatial clusters of LST. In short, high and low values tended to agglomerate and create a sort of cluster. Also, the clustering pattern expanded towards the north and southwest during the study period.
The results are useful and satisfactory because, even with the great cloud cover existing in the Amazon, it was possible to detect the spatial characteristics of the LST for the city of Belém. It was also evident that the disordered occupation of land use, caused the hot spots found. Future perspectives on this topic in the city of Belém, should address the evolution of land use and coverage and its relationship with air temperature, as well as the application of a trend test in these series. Our work can be useful to decision makers in the strategic planning of land use and environmental monitoring in Belém.
The LST mapping revealed that coldspots were more dominant than hotspots. However, the trends of spatial formation are directly associated with the land use and occupation dynamics. Regarding the spatial distribution analysis based on the global Moran's I index indicated the rejection of the null hypothesis, due to the observation of spatial clusters of LST. In short, high and low values tended to agglomerate and create a sort of cluster. Also, the clustering pattern expanded towards the north and southwest during the study period. Future studies can use the results of this study to estimate LST at a resolution of 30 x 30 m taking the bands that determined the biophysical indices.