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

Aplicação da Modelagem Espacial de Variáveis Biofísicas em uma área Urbanizada na Amazônia: O Caso da Região Metropolitana de Belém-Pará

Abstract

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.

Keywords
Land use; urban heat island; remote sensing; geoprocessing

Resumo

O tipo e uso cobertura do solo desempenham papel decisivo para a temperatura da superfície continental (LST). Como as cidades são compostas por coberturas variadas, incluindo vegetação, áreas construídas, prédios, estradas e áreas desprovidas de vegetação, compreender os padrões da LST no complexo espaço urbano se faz cada vez mais importante. O presente estudo investigou a relação da LST com o albedo, NDVI, NDWI, NDBI e NDBaI no período de 1994 e 2017. Foram utilizadas imagens dos sensores Thematic Mapper (TM) e Thermal Infrared Sensor (TIRS) a bordo dos satélites Landsat 5 e 8, respectivamente. As imagens foram processadas, reamostradas (resolução espacial de 120 m) no ambiente do software QGIS 3.0 e, por fim, foram extraídos centroides com um total de 1252 pontos. O modelo de regressão clássica (RC) foi aplicado às variáveis seguidas por modelos espaciais autoregressivos e de erro espacial (MEAR e MEE) e os resultados foram comparados a partir de índices de acurácia. Os resultados mostram que o maior coeficiente de correlação existe foi verificado entre albedo e NDBaI (r = 0,88). A relação entre o albedo e a LST, com r = 0,7, também é positiva e significativa ao nível de (р < 0,05). A partir do índice I de Moran Global verificou-se a dependência espacial e não estacionariedade da LST (I = 0,44). O modelo MEE apresentou as melhores métricas de acurácia (AIC = 3307,15 e R² = 0,65) explicando consideravelmente mais variações na relação dos fatores explicativos e da LST para a região metropolitana de Belém, quando comparados aos modelos convencionais de RC.

Palavras-chave
Uso do solo; Ilha de calor urbana; sensoriamento remoto; geoprocessamento

1. 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 man-made 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 (Zhao et al., 2018), 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 and Pé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.

2. Material and Methods

2.1. 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², and its population was estimated in 1,452,275 inhabitants and demographic density in 1,315.26 inhab./km² in 2017. The present study focused on the urban area of Belém, covering 184 km². The rural areas and river islands within the territory of Belém metropolis were excluded from this research.

Figure 1
Geographic location of the study area.

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

Table 1
Information provided by image metadata.

Dates	Orbit	Point	Central time (UTC)	Local time (Brasília)	Solar elevation angle
21/06/1994	223	61	12:41	09:41	46.2
13/07/2008	223	61	13:09	10:09	51.8
06/07/2017	223	61	13:22	10:22	53.3

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

2.3. Data processing

2.3.1. Landsat-5

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 (1987), that is Eq. (1):

(1) $$L_{\lambda i}=a+\left(\frac{b-a}{255}\right)\times ND$$

where: L_λi (W m⁻² sr⁻¹ µm⁻¹) 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):

(2) $$\rho =\frac{L_{\lambda }\times \pi }{K_{\lambda }\times \cos Z\times d_r}$$

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⁻²µm⁻¹); 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:

(3) $${\alpha }_{toa}=\left(0.293\times {\rho }_1\right)+\left(0.274\times {\rho }_2\right)+\left(0.233\times {\rho }_3\right)+\left(0.157\times {\rho }_4\right)+\left(0.033\times {\rho }_5\right)+\left(0.011\times {\rho }_7\right)$$

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

(4) $${\alpha }_{sup}=\frac{{\alpha }_{toa}-{\alpha }_{atm}}{{\tau }_{sw}^2}$$

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

(5) $${\tau }_{SW}=0.75+2\times {10}^{-5}\times z$$

where: z represents the altitude (m) of each pixel.

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

(6) $${\rho }_{\lambda ,b}=\frac{H_{\rho }\times Q_{CAL}+A_{\rho }}{\cos \left({\theta }_{SE}\right)}$$

where ρ_λ,b is the reflectance of band b; H_ρ is the specific multiplicative scaling of each band (2⁻⁵, 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):

(7) $${\theta }_{SE}=90-{\theta }_{SZ}$$

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

(8) $${\alpha }_{toa}={\sum }^{\text{}}\left({\varpi }_b\times {\rho }_{\lambda ,b}\right)$$

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 (∑ESUN_λ,b).

Therefore, the surface albedo or corrected albedo (α_sup) is calculated using the same expressions abovementioned for Landsat-5 (Eq. (4)).

2.3.3. 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:

(9) $$L_{\lambda 10}=M_L\times Q_{cal}+A_L$$

where: L_λ10 is the radiance of the thermal band; M_L is the band-specific multiplicative rescaling factor, equal to 3,342 × 10⁻⁴ Wm⁻²sr⁻¹µm⁻¹; 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):

(10) $$SAVI=\frac{\left(1+L)\times ({\rho }_{IV}-{\rho }_V\right)}{\left(L+{\rho }_{IV}+{\rho }_V\right)}$$

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

(11) $$LAI=-\frac{ln\left(\frac{0.69-SAVI}{0.59}\right)}{0.91}$$

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

(12) $${\epsilon }_{NB}=0.97+0.0033\times LAI$$

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

(13) $$LST=\frac{K_2}{ln\left(\frac{{\epsilon }_{NB}{K}_1}{L_{\lambda thermal}}+1\right)}$$

where: the K₁ and K₂ constants of the thermal band of Landsat-5 are 607.76 Wm⁻²sr⁻¹µm⁻¹ and 1260.56 K, and of Landsat-8 are 774.89 Wm⁻²sr⁻¹µm⁻¹ and 1321.08 K, respectively.

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

(14) $$NDVI=\frac{\left({\rho }_{IVP}-{\rho }_{VRM}\right)}{\left({\rho }_{IVP}+{\rho }_{VRM}\right)}$$

where ρ represents the reflectance of each band, IVP (near infrared), VRM (red) and IVM (medium infrared).

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.

(15) $$NDBI=\frac{\left(d_{IVM}-d_{IVP}\right)}{\left(d_{IVM}+d_{IVP}\right)}$$

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

(16) $$NDBaI=\frac{\left(d_{IVM}-d_{IVT}\right)}{\left(d_{IVM}+d_{IVT}\right)}$$

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.

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

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

(17) $$I=\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n{w}_{ij}\left(y_i-\overline{y}\right)\left(y_j-\overline{y}\right)}{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

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; $\overline{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.

2.4.2. Geographically weighted regression models

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

(18) $$Y_i={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\dots +{\beta }_n{X}_n+{\epsilon }_i$$

where Y_i is an observation of the dependent variable; X₁, X₂,…, X_n are the independent variables; β = (β₀, β₁, β₂,…,β_n) are referred to as corresponding regression coefficients, and ${\epsilon }_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).

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

(19) $$Y_i=\rho {\displaystyle \sum }_{j=1}^n{W}_{ij}{y}_i+{\displaystyle \sum }_{q=1}^Q{X}_{iq}{\beta }_q+{\epsilon }_i$$

where ${\epsilon }_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 ∑_j W_ijy_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.

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

(20) $${\epsilon }_i=\rho {\displaystyle \sum }_{j=1}^n{W}_{ij}{\epsilon }_j+u_i$$

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.

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

(21) $$AIC=D\left(y;\widehat{\mu }\right)+2p$$

where $D\left(y;\widehat{\mu }\right)$ 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 $\widehat{\beta }$. 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.

3. Result and Discussion

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

Figure 2
Scatterplot between average values (1994, 2008 and 2017) of Albedo, NDVI, NDWI, NDBaI and LST for the urban area of Belém Metropolis.

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), Liu and Zhang (2011) and Ranagalage 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).

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

Figure 3
Scatterplot of Moran's I index for the average (1994, 2008 and 2017) pixel value.

Several researches have used the Moran's I index (Deilami et al., 2018; Mavrakou et al., 2018 and Yang and Santamouris, 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 represented 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.

Figure 4
Land surface temperature cluster map.

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.

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

Table 2
Coefficients and Moran's I residuals of the regression models fitted for the LST.

Covariables	CR (p-value)	SARM (p-value)	SEM (p-value)
Intercept	25.28 (< 0.001)	18.50 (< 0.001)	26.85 (< 0.001)
Albedo	4.84 (< 0.001)	–	5.05 (< 0.001)
NDVI	–1.99 (< 0.001)	–1.97 (< 0.001)	–3.83 (< 0.001)
NDBI	6.71 (< 0.001)	5.91 (< 0.001)	–
NDWI	–	–	–3.81 (< 0.001)
R2	0.43	0.51	0.65
AIC	3758.94	3599.47	3307.15
BIC	3779.47	3620	3327.68
ρ	–	0.30 (< 0.001)	0.62 (< 0.001)
Test residual	0.45 (< 0.001)	0.22 (< 0.001)	–0.050 (< 0.001)

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 coefficient of determination (R² = 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² 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. According 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.

Figure 5
Residual maps for the classical regression (CR), spatial autoregressive (SARM), and spatial error (SEM) models for surface temperature in the urban mesh of the city of Belém.

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.

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

Acknowledgments

The authors would like to thank the Academic Unit of Atmospheric Sciences of the Federal University of Campina Grande (UFCG) for the structure provided and also the Coordination for the Improvement of Higher Education Personnel (CAPES), for the scholarship granted to the first author and for the financing of the Research Project entitled “Analysis and Prediction of Intense Hydrometeorological Phenomena in the East of Northeast Brazil” - Pro-Alertas Notice n° 24/2014 (Case No. 88887.091737/2014-01). The authors also thank the National Council for Scientific and Technological Development (CNPq) for financing Research Projects under numbers 446172 / 2015-4 and 409499 / 2018-8, as well as the Research Productivity Scholarship (Process No. 304493 / 2019-8).

Reference

• ALLEN, R.G.; TASUMI, M.; TREZZA, R.; BASTIAANSSEN, W. SEBAL (Surface Energy Balance Algorithms for Land), Idaho Implementation: Advanced Training and User's Manual. NASA EOSDIS/Raytheon Company/Idaho Department of Water Resources, 2002, 97p.
• ALMEIDA, A.S.D.; MEDRONHO, R.D.A.; VALENCIA, L.I.O. Análise espacial da dengue e o contexto socioeconômico no município do Rio de Janeiro, RJ. Revista de Saúde Pública, v. 43, p. 666-673, 2009.
• ARULBALAJI, P.; PADMALAL, D.; MAYA, K. Impact of urbanization and land surface temperature changes in a coastal town in Kerala, India. Environ. Earth Sci, v. 79, p. 400, 2020.
• BARRON, L.; RUGGIERI, D.; BRANAS, C. Assessing Vulnerability to Heat: A Geospatial Analysis for the City of Philadelphia. Urban Science, v. 2, n. 2, p. 38, 2018.
• BASTIAANSSEN, W.G.M. SEBAL-based sensible and latent heat fluxes in the irrigated Gediz Basin. Turkey. Journal of Hidrology, v. 229, p. 87-100, 2000.
• BIVAND, R.; PIRAS, G. Comparing implementations of estimation methods for spatial econometrics. Journal of Statistical Software, v. 63, n. 18, p. 1-36, 2015.
• CâMARA, G.; CARVALHO, M.S.; CRUZ, O.G.; CORREA, V. Análise Espacial de áreas. Análise Espacial de Dados Geográficos. Brasília: Embrapa, p. 157-8, 2004.
• CHEN, X.L.; ZHAO, H.M.; LI, P.X. YIN, Z.Y. Remote sensing image-based analysis of the relationship between urban heat island and land use/cover changes. Remote Sensing of Environment, v. 104, n. 2, p. 133-146, 2006.
• CHOUDHURY, D., DAS, K., DAS, A. Assessment of land use land cover changes and its impact on variations of land surface temperature in Asansol-Durgapur development region. Egypt J. Remote Sens. Sp. Sci, v. 22, p. 203-218, 2019.
• CORRAR, L.J.; PAULO, E.; DIAS FILHO, J.M. Análise Multivariada: Para os Cursos de Administração, Ciências Contábeis e Economia São Paulo: Atlas, 2007.
• COUTTS, A.M., HARRIS, R.J., PHAN, T., LIVESLEY, S.J., WILLIAMS, N.S., TAPPER, N.J. Thermal infrared remote sensing of urban heat: Hotspots, vegetation, and an assessment of techniques for use in urban planning. Remote Sens. Environ, v. 186, p. 637-651, 2016.
• DEILAMI, K.; KAMRUZZAMAN, M.; LIU, Y. Urban heat island effect: a systematic review of spatio-temporal factors, data, methods, and mitigation measures. Int. J. Appl. Earth Obs. Geoinf, v. 67, p. 30-42, 2018
• DE OLIVEIRA SERRãO, E.A.; SILVA, M.T.; FERREIRA, T.R.; DE ATAIDE, L.C.P.; WANZELER, R.T.S.; DA SILVA, V.P.R., DE LIMA; A.M.M.; DE SOUSA F.D.A.S. Land use change scenarios and their effects on hydropower energy in the Amazon. Science of The Total Environment, v. 744, p. 140981, 2020.
• DI LEO, N.; ESCOBEDO, F.J.; DUBBELING, M. The role of urban green infrastructure in mitigating land surface temperature in Bobo-Dioulasso, Burkina Faso. Environ. Dev. Sustain, v. 18, n 2, p. 373-392, 2016.
• DOBSON, A.J.; BARNETT, A. An Introduction to Generalized Linear Models Texts in Statistical Science, Boca Raton, FL: Chapman & Hall/CRC Press, 2011.
• DUNCAN, J.M.A., BORUFF, B., SAUNDERS, A., SUN, Q., HURLEY, J., AMATI, M. Turning down the heat: an enhanced understanding of the relationship between urban vegetation and surface temperature at the city scale. Sci. Total Environ, v. 656, p. 118-128, 2019.
• FISCHER, M.M.; LESAGE, J. Spatial Econometric Methods For Modeling Origin Destination Flows. Handbook of Applied Spatial Analysis Springer, Berlin/Heidelberg, p. 409-432, 2010
• FOTHERINGHAM, A.S.; CRESPO, R.; YAO, J. Geographical and Temporal Weighted Regression (GTWR). Geogr Anal, v. 47, p. 431-452, 2015. doi:10.1111/gean.12071
  » https://doi.org/10.1111/gean.12071
• FU, P.; WENG, Q. A time series analysis of urbanization induced land use and land cover change and its impact on land surface temperature with Landsat imagery. Remote Sens. Environ, v. 184, p. 175-187, 2016.
• GARCIA, P.; PEREZ, E. Mapping of soil sealing by vegetation indexes and built-up index: a case study in Madrid (Spain). Geoderma, v. 268, p. 100-107, 2016.
• GRIMM, N.B.; GROVE, J.G.; PICKETT, S.T.A.; REDMAN, C.L. Integrated Approaches to Long-Term Studies of Urban Ecological SystemsUrban ecological systems present multiple challenges to ecologists-pervasive human impact and extreme heterogeneity of cities, and the need to integrate social and ecological approaches, concepts, and theory. Bioscience, v. 50, p. 571-584, 2000.
• GUHA, S.; GOVIL, H.; DEY, A.; GILL, N. Analytical study of land surface temperature with NDVI and NDBI using Landsat 8 OLI and TIRS data in Florence and Naples city, Italy. Eur. J. Remote Sens, v. 51, p. 667-678, 2018.
• GUO, G.; WU, Z.; XIAO, R.; CHEN, Y.; LIU, X.; ZHANG, X. Impacts of urban biophysical composition on land surface temperature in urban heat island clusters. Landsc. Urban Plan, v. 135, p. 1-10, 2015.
• JAN, G.; MICHAL, L.; STEVAN, S.; Dragan, M. Modelled spatiotemporal variability of outdoor thermal comfort in local climate zones of the city of Brno, Czech Republic. Sci. Total Environ, v. 624, p. 385-395, 2018.
• HO, H.; KNUDBY, A.; HUANG,W. A spatial framework tomap heat health risks at multiple scales. Int. J. Environ. Res. Public Health, v. 12, p. 15046, 2015.
• HUETE, A., DIDAN, K., MIURA, T., RODRIGUEZ, E.P., GAO, X., FERREIRA, L.G. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ, v. 83, p. 195-213, 2002.
• LI, J.; SONG, C.; CAO, L.; ZHU, F.; MENG, X.; WU, J. Impacts of landscape structure on surface urban heat islands: a case study of Shanghai, China. Remote Sens. Environ, v. 115, p. 3249-3263, 2011.
• LI, G.; ZHANG, X.; MIRZAEI, P.A.; ZHANG, J.; ZHAO, Z. Urban heat island effect of a typical valley city in China: responds to the global warming and rapid urbanization. Sustain. Cities Soc, v. 38, p. 736-745, 2018.
• LI, W.; CAO, Q.; LANG, K.; WU, J. Linking potential heat source and sink to urban heat island: Heterogeneous effects of landscape pattern on land surface temperature. Sci. Total Environ, v. 586, p. 457-465, 2017.
• LIU, C.; LIU, J.; HU, Y.; WANG, H.; ZHENG, C. Airborne thermal remote sensing for estimation of groundwater discharge to a river. Groundwater, v. 54, p. 363-373, 2016.
• LIU, L.; ZHANG, Y. Urban heat island analysis using the Landsat TM data and ASTER data: A case study in Hong Kong. Remote Sensing, v. 3, n. 7, p. 1535, 2011.
• LUO, X.; PENG, Y. Scale effects of the relationships between urban heat Islands and impact factors based on a geographicallyweighted regression model. Remote Sensing, v. 8, n. 9, p. 760, 2016.
• MANCINO, G.; NOLE, A.; RIPULLONE, F.; FERRARA, A. Landsat TM imagery and NDVI differencing to detect vegetation change: assessing natural forest expansion in Basilicata, southern Italy. iForest-Biogeosciences and Forestry, v. 7, p. 75-84, 2014.
• MARGALHO, E.S. SILVA, M.T.; CARDOSO, L.K.S.; OLINDA, R.A.; MENEZES, J.F.G. Influência da mudança do uso e cobertura do solo sobre a temperatura da superfície continental na área urbana de Belém-PA. Anuário do Instituto de Geociências - UFRJ, v. 43, p. 7-19, 2020.
• MARKHAM, B.L.; BARKER, J.L. Thematic mapper band pass solar exoatmospherical irradiances. International Journal of Remote Sensing, v. 8, n. 3, p. 517-523, 1987.
• MAVRAKOU, T.; POLYDOROS, A.; CARTALIS, C.; SANTAMOURIS, M. Recognition of thermal hot and cold spots in urban areas in support of mitigation plans to counteract overheating: application for Athens. Climate, v. 6, p. 16, 2018.
• NICHOL, J. An emissivity modulation method for spatial enhancement of thermal satellite images in urban heat island analysis. Photogrammetric Engineering & Remote Sensing, v. 75, n. 5, p. 547-556, 2009.
• PACIFICI, M.; RAMA, F.; MARINS, K.R.D.C. Analysis of temperature variability within outdoor urban spaces at multiple scales. Urban Clim, v. 27, p. 90-104, 2019.
• PUREVDORJ, T.S.; TATEISHI, R.; ISHIYAMA, T.; HONDA, Y. Relationships between percent vegetation cover and vegetation indices. International Journal of Remote Sensing, v. 19, p. 3519-3535, 1998.
• RANAGALAGE, M.; ESTOQUE, R.C.; MURAYAMA, Y. An Urban Heat Island Study of the Colombo Metropolitan Area, Sri Lanka, Based on Landsat Data (1997-2017). ISPRS Int. J. Geo-Inf., v. 6, p. 189, 2017.
• ROUSE, J.; HAAS, R.; SCHELL, J.; DEERING, D.; HARLAN, J. Monitoring the Vernal Advancement of Retrogradation of Natural Vegetation: Final Report [online] [Viewed 6 February 2016], Available from: NASA/GSFC, Greenbelt, pp. 371, 1974.
• ROUSE, J.W.; HAAS, R.H.; SCHELL, J.A. & DEERING, D.W. Monitoring vegetation systems in the Great Plains with ERTS. In: 3rd Earth Resources Technology Satellite Symposium, p. 309-317, 1973.
• SILVA, M.T.; SILVA, V.P.R.; SOUZA, E.P.; ARAúJO, A.L. Aplicação do modelo SWAT na estimativa da vazão na bacia hidrográfica do submédio rio São Francisco. Revista Brasileira de Geografia Física, v. 8, n. 6, p. 1615-1627, 2015.
• SILVA, M.T.; SILVA, V.P.R.; SILVA, M.M.M.A.; SILVA, H.C.D.; OLIVEIRA, N.F. Space time variability of surface temperature in the semi-arid Pernambuco based image TM/Landsat. Journal of Hyperspectral Remote Sensing, v. 4, p. 111-120, 2014.
• SILVA, V.P.R.; SILVA, M.T.; SOUZA, E.P. Influence of land use change on sediment yield: a case study of the sub-middle of the São Francisco river basin. Engenharia Agrícola, v. 36, n. 6, p. 1005-1015, 2016.
• SILVA, V.P.R.; SILVA, M.T.; BRAGA, C.C.; SINGH, V.P.; SOUZA, E.P.; SOUSA, F.A.S.; HOLANDA, R.M.; ALMEIDA, R.S.R.; BRAGA, A.C.R. Simulation of stream flow and hydrological response to land-cover changes in a tropical river basin. Catena, v. 162, p. 166-176, 2018.
• TASUMI, M.; ALLEN, R.G.; TREZZA, R.; WRIGHT, J.L. Satellite-based energy balance to assess within-population variance of crop coefficient curve. Journal of Irrigation and Drainage Engineering, v. 131, n. 1, p. 95-108, 2008.
• TYRALIS, H.; MAMASSIS, N.; PHOTIS, Y.N. Spatial analysis of the electrical energy demand in Greece. Energy Policy, v. 102, p. 340-352, 2017.
• VARHNEY, A. Improved NDBI differencing algorithm for built-up regions change detection from remote-sensing data: an automated approach. Remote. Sens. Lett, v. 4, p. 504-512, 2013.
• WANG, J.F.; HU, Y. Environmental health risk detection with GeogDetector. Environ. Model Softw, v. 33, p. 114-115, 2012.
• WANDERLEY, L.N.; DOMINGUES, R.M.L.; JOLY, A.; DA ROCHA, C.R.H. Relationship between land surface temperature and fraction of anthropized area in the Atlantic forest region, Brazil. PLOS ONE, v. 14, n. 12, e0225443, 2019.
• YANG, J.; SANTAMOURIS, M. Urban heat island and mitigation technologies in Asian and Australian cities - impact and mitigation. Urban Sci, v. 2, p. 74, 2018.
• YUAN, F.; BAUER, M.E. Comparison of impervious surface area and normalized difference vegetation index as indicators of surface urban heat island effects in Landsat imagery. Remote Sens. Environ, v. 106, p. 375-386, 2007.
• YUE, W.; XU, J.; TAN, W.; XU, L. The relationship between land surface temperature and ndvi with remote sensing: Application to shanghai Landsat 7 ETM+ data. Int. J. Remote Sens, v. 28, p. 3205-3226, 2007.
• ZHA, Y.; GAO, J.N.I.S. Use of normalized difference built-up index in automatically mapping urban areas from TM imagery. International Journal of Remote Sensing, v. 24, p. 583-594, 2003.
• ZENG, Y.; HUANG, W.; ZHAN, F.B. Study on the urban heat island effects and its relationship with surface biophysical characteristics using MODIS imageries. Geo-spat. Inf. Sci., v. 13, p. 1-7, 2010.
• ZHANG, X.; ESTOQUE, R.C.; MURAYAMA, Y. An urban heat island study in Nanchang City, China based on land surface temperature and social-ecological variables. Sustain. Cities Soc, v. 32, p. 557-568, 2017.
• ZHAO, H.M.; CHEN, X.L. Use of normalized difference bareness index in quickly mapping bare areas from TM/ETM+. Geoscience and Remote Sensing Symposium, v. 3, p. 1666-1668, 2005.
• ZHAO, H.; REN, Z.; TAN, J. The spatial patterns of land surface temperature and its impact factors: spatial non-stationarity and scale effects based on a geographicallyweighted regression model. Sustain. For, v. 10, p. 2242, 2018.
• ZHOU, X.; WANG, Y.C. Dynamics of land surface temperature in response to land-use/cover change. Geographical Research, v. 49, n. 1, p. 23-36, 2011.

Internet Resources

• GeoR Package, https://cran.r-project.org/web/packages/geoR/index.html
  » https://cran.r-project.org/web/packages/geoR/index.html

Publication Dates

• Publication in this collection
  23 Apr 2021
• Date of issue
  Apr-Jun 2021

History

• Received
  22 May 2020
• Accepted
  13 Jan 2021