Assessing the Relationship of LST, NDVI and EVI with Land Cover Changes in the Lagos Lagoon Environment

Medium resolution Landsat imageries for 4 years were acquired from the United States Geological Surveys Earth Explorer portal (USGS 2020). On the Earth Explorer webpage, the appropriate search criteria were set and the imageries were downloaded. Following the approach of Ferrelli et al. (2018), imageries with excessive cloud cover were discarded. To maximise the coverage of the available scenes, a combination of Landsat scenes (191/055 and 191/056) was done to form a complete coverage of the study area at four average periods: 1984 (Thematic Mapper (TM)), 2002 (Enhanced Thematic Mapper (ETM+)), 2013 and 2019 (Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS)). These imageries are already geometrically corrected by USGS and ortho-rectified to Level 1 (Tatem et al. 2006, Gutman et al. 2013). A further confirmation of the geometric/positional suitability of the Landsat scenes was done with reference to the Lagos State boundary along the Lagos shoreline through overlay and the alignment was satisfactory. The scenes used for land cover classification covered the following periods: TM (1984–1986), ETM+ (2001/2002) and OLI (December 2013 and January 2019); while scenes used for LST determination covered the following: TM (December 1984), ETM+ (February 2000) and TIRS (December 2013 and January 2018). For NDVI and EVI extraction, the imagery dates are as follows: TM (December 1984), ETM+ (December 2002) and OLI (December 2013 and January 2019). Generally, the imageries cover dry season months. The Landsat imageries have a pixel spacing of 30 m and are referenced to a Universal Transverse Mercator (UTM) coordinate system defined on the WGS84 datum (WGS84 UTM Zone 31N).

Following the approach of Ullah et al. (2019), the first step involved a preliminary interpretation of the Landsat imageries aided by Anderson's (1971) Level I Classification Scheme and the use of Google Earth imagery for training sites development. This interpretation categorised the study area land cover into four classes, such as bare land, built-up area, vegetation and water body (Table 1). Next, training sites were created for each class and the imageries were then subjected to supervised classification by the Maximum Likelihood Classifier (MLC) in ENVI 5.3 image processing software. MLC incorporates the variance–covariance within the class distributions and for data that has a normal distribution; it gives a better performance than other known parametric classifiers. It performs well over various types of land cover, satellite systems and conditions (Bolstad, Lillesand 1991). It is a preferred choice for many applications and its good performance has been proven (Sharma et al. 2017, Brendel et al. 2019). After classification, the feature classes were transferred to ArcGIS for editing, refinement and assessment of classification accuracy. According to Hasmadi et al. (2009), accuracy assessment is important for checking the reliability of the classification results. The main idea behind any image classification process is to obtain the highest accuracy possible (Sharma et al. 2011). The minimum level of interpretation accuracy in the identification of land cover classes from remotely sensed data should be at least 85% (Anderson 1971). Using the Landsat imagery as a reference, 90 points were used to compare features interpreted on the imageries and their corresponding output in the classification. The point comparison detected the correctly classified pixels, pixels assigned to a certain class that did not belong to it (commission errors) and pixels that belong to one class but are included into other classes (omission errors) (Tran et al. 2017, Aboelnour, Engel 2018). The overall accuracy was then computed using the following expression: (1) Overall accuracy (%) = Total number of correct samplesTotal number of samples×100 {\rm{Overall}}\,{\rm{accuracy}}\,(\%) \, = \,{{{\rm{Total}}\,{\rm{number}}\,{\rm{of}}\,{\rm{correct}}\,{\rm{samples}}} \over {{\rm{Total}}\,{\rm{number}}\,{\rm{of}}\,{\rm{samples}}}} \times 100

Traditionally, the NDVI is calculated using a combined operation between the Red band and Near-Infrared (NIR) band (Qiu et al. 2018, Ferrelli et al. 2018, Ullah et al. 2019, Guha et al. 2020) as follows: (2) NDVI = (RNIR−RRed)/(RNIR+RRed) {\rm{NDVI}}\, = \,\left({{R_{{\rm{NIR}}}} - {R_{{\rm{Red}}}}} \right)/\left({{R_{{\rm{NIR}}}} + {R_{{\rm{Red}}}}} \right) where: –

R_NIR represents the spectral reflectance in the NIR band,

The EVI is computed (Hoek van Dijke et al. 2019, Semeraro et al. 2020) as follows: (3) EVI=G×((RNIR−RRed)/(RNIR+C1×RRed−C2×RBlue+L)) {\rm{EVI}} = {\rm{G}} \times \left({\left({{R_{{\rm{NIR}}}} - {R_{{\rm{Red}}}}} \right)/\left({{R_{{\rm{NIR}}}} + {C_1} \times {R_{{\rm{Red}}}} - {C_2} \times {R_{{\rm{Blue}}}} + L} \right)} \right) where: –

G is the Gain factor,

These enhancements incorporated into EVI reduce the background noise, atmospheric noise and saturation in most cases.

Landsat Level 2 NDVI and EVI products were ordered and downloaded from the USGS Earth Resources Observation and Science (EROS) Centre Earth Science Processing Architecture on-demand interface. The ESPA is an incubation environment that provides users with an on-demand interface to process and customise Landsat science products (USGS 2020). The coefficients adopted on ESPA in the EVI algorithm for Landsat 4–7 and Landsat 8 are G = 2.5, C₁ = 6, C₂ = 7.5 and L = 1. The Level 1 scene lists were uploaded on the interface and processing was done for the surface reflectance, NDVI and EVI Level 2 products. The surface reflectance enhances comparison between multiple Landsat imageries over the study area by considering atmospheric effects such as aerosol scattering and thin clouds, which can help in the detection and characterisation of surface changes. According to USGS (2019), surface reflectance is generated from Level-1 inputs that meet the <76°solar zenith angle constraint and include the required auxiliary data inputs to generate a scientifically viable product. The NDVI and EVI products were delivered from ESPA as single bands with 16-bit integers and values ranging from −10,000 to +10,000 (USGS 2020). However, using the given scale factor of 0.0001, this was scaled down to a range from −1 to +1.

The retrieval of LST followed the single-channel method (Oguz 2013, Ferrelli et al. 2015, Obiefuna et al. 2018). Landsat TM thermal Band 6, ETM Band 6_1 and TIRS Band 10 were used for the retrieval. The key stages of the methodology are discussed below.

The formula for converting DN in Landsat 5 and 7 to spectral radiance is given by Zareie et al. (2016): (4) Lλ=(LMAX−LMINQCALMAX−QCALMIN)×(QCAL−QCALMIN)+LMIN {{\rm{L}}_{\rm{\lambda}}} = \left({{{{\rm{LMAX}} - {\rm{LMIN}}} \over {{\rm{QCALMAX}} - {\rm{QCALMIN}}}}} \right) \times \left({{\rm{QCAL}} - {\rm{QCALMIN}}} \right) + {\rm{LMIN}} where: –

L_λ is spectral radiance at the sensor's aperture (W m⁻² sr⁻¹ μm⁻¹),

The formula to derive the spectral radiance for Landsat 8 is given by USGS (2015): (5) Lλ−ML×QCAL+AL {L_{\rm{\lambda}}} - {M_L} \times {\rm{QCAL}} + {A_L} where: –

M_L is Radiance multiplicative scaling factor for the band,

The values for LMIN, LMAX, QCALMIN, QCALMAX, M_L and A_L are derived from the Landsat metadata file.

After calculating the spectral radiance (L_λ), the next step was the calculation of the TOA brightness temperature. The TOA approximation formula is given by Dewan and Corner (2012), Zareie et al. (2016), Hamoodi et al. (2019) and Guha et al. (2020): (6) T=K2/log(1+K1/Lλ) T = {K_2}/\log \left({1 + {K_1}/{L_{\rm{\lambda}}}} \right) where: –

T is TOA brightness temperature (K),

Values for K₁ and K₂ for Landsat TM and ETM+ are shown in Table 2, whereas Table 3 shows the values for Landsat 8 TIRS.

The brightness temperature was subsequently converted to LST using Eq. (7) (Hamoodi et al. 2019, Ullah et al. 2019): (7) ST=T1+(λ×T/ρ)logε {S_T} = {T \over {1 + \left({{\rm{\lambda}} \times T/\rho} \right)\log \varepsilon}} where: –

S_T is LST (K),

Finally, the LST in Kelvin was converted to degree Celsius by subtracting from 273.15.

Statistical analyses are commonly applied to establish relationships between parameters, for example, Land cover-LST (Dewan, Corner 2012), LST-NDVI (Ferrelli et al. 2018, Hamoodi et al. 2019, Malik et al. 2019, Ullah et al. 2019, Guha et al. 2020, Mukherjee, Singh 2020), NDVI-EVI (Chen et al. 2006b, Matsushita et al. 2007, Li et al. 2010, Uyeda et al. 2017) and LST-EVI (Phompila et al. 2015). In this study, an inventory of the parameters was created for the study area. To do this, the LST, NDVI and EVI raster maps were converted to grids of point shape files with the associated parameter values at pixel centres in the attribute table. The grids were then overlaid on the land cover, LST and NDVI maps within ArcGIS. The values on the land cover maps coinciding with the points were extracted. To explore the relationship between the parameters, descriptive statistics were computed using the Statistical Package for the Social Sciences (SPSS) software v20. In line with Guha et al. (2020), the minimum, maximum, mean and standard deviations (SDs) of the parameters were generated. Pearson's correlation analysis was used to analyse the correlation between variables. Linear regression was also carried out to determine the temporal relationship between the LST, NDVI and EVI in different land cover classes. The effect of land cover on the LST was calculated using the Contribution Index (CI). The CI quantifies the warming or cooling extent of a land cover type and this is related to the proportion of the total land area it occupies. It links spatial structure and long-term changes in land cover to LST intensities. The equation for calculating the CI is given by Eq. (8) (Odindi et al. 2015, Odindi et al. 2017, Tarawally et al. 2018). (8) CI=Dt×S CI = {D_t} \times S where: –

D_t is the average LST of study area minus average LST of land cover class type,

If CI of a land cover type is positive, it indicates that it contributes to raising the LST and vice versa.

Going further, a trend analysis was carried out using the Mann-Kendall test and Sen's estimator of slope. The Mann-Kendall's test is a non-parametric test for trend analysis (Gilbert 1987). It is particularly useful as the data need not conform to any particular distribution (Hamed 2008). The use of mean values of multiple observations in a period was recommended by Gilbert (1987) for Mann-Kendall trend analysis. If the Mann-Kendall's statistic (S) is a large positive number, measurements taken later in time tend to be larger than those taken earlier. Similarly, if S is a large negative number, measurements taken later in time tend to be smaller. The first scenario is to test the null hypothesis (H_o) of no trend against the alternative hypothesis (H_A) of an upward trend. H_o is rejected in favour of H_A if S is positive and if the probability value is less than the significance level (p-value) of the test. Similarly, to test H_o against the alternative hypothesis H_A of a downward trend, Ho is rejected and H_A accepted if S is negative and if the probability value is less than the a priori specified p-value. For this study, a trend analysis of the LST, NDVI and EVI was done for the years 1984, 2002, 2013 and 2019. The Sen's slope estimator, Q_med (Sen 1968, Gocic, Trajkovic 2013), is the magnitude of the upward trend per year. The variables with the upward trend or linear trend were further probed for the annual increase or slope value. The Q_med sign reflects data trend reflection, while its value indicates the steepness of the trend.

Table 4 shows the areal distribution of land cover in the study area for the years 1984, 2002, 2013 and 2019. The analysis shows that bare land increased by 613.7% from 7.08 km² in 1984 to 50.53 km² in 2002 (at the rate of 2.41 km² a⁻¹). Thereafter, it decreased by 72.99% to 13.65 km² between 2002 and 2013 (at the rate of 3.35 km² a⁻¹) and finally decreased by 82.05% to 24.85 km² between 2013 and 2019 (at the rate of 1.87 km² a⁻¹). Built-up areas that occupied only 40.77 km² in 1984 expanded significantly by 219.7% to 130.34 km² between 1984 and 2002 (at the rate of 4.98 km² a⁻¹) and subsequently increased by 219.37% to 416.27 km² between 2002 and 2013 (at the rate of 25.99 km² a⁻¹). Between 2013 and 2019, it increased slightly by 6.06% to 441.49 km² (at the rate of 4.20 km² a⁻¹). There were depletions in the coverage of vegetation within the area. Between 1984 and 2002, 9.25% of the vegetation cover (122.28 km²) was lost at the rate of 6.79 km² a⁻¹. In the period between 2002 and 2013, 20.24% of vegetation cover representing 242.82 km² was lost at the rate of 22.07 km² a⁻¹. Finally, 2.39% of vegetation cover (22.85 km²) was lost between 2013 and 2019 (at the rate of 3.81 km² a⁻¹). In relation to other land cover classes, the loss in water body was moderate over the entire period. Water body decreased from 415.37 km² in 1984 to 404.63 km² in 2002, 398.40 km² in 2013 and finally dropped to 384.83 km² in 2019.

The accuracy assessment of the classification yielded the following overall accuracies: 96.67% (1984), 81.11% (2002), 93.33% (2013) and 87.78% (2019).

Tables 5–8 show descriptive statistics of the LST, NDVI and EVI per land cover excluding water body for the years 1984, 2002, 2013 and 2019, respectively. The mean LSTs for the four periods are 22.68°C (1984), 24.34°C (2002), 26.46°C (2013) and 28.40°C (2019), respectively. Comparison was done with historical monthly temperatures of Lagos from WWO (2020). The online temperature data showed that the state had average temperatures of 28°C in both December 2013 and January 2018. These monthly averages are in the same range of the average LSTs from Landsat shown in Tables 7 and 8. The trend reveals a gradual rise in surface temperatures and warmness within the environment surrounding the Lagos Lagoon over the years. The lowest temperature recorded was 17.63°C in 1984 whereas the highest temperature was 36.90°C in 2019. The built-up area class had the highest mean LST values while vegetation had the lowest mean values. The NDVI values ranged from 0.0 to 0.5 (1984), −0.2 to 0.70 (2002), −0.12 to 0.85 (2013) and −0.06 to 0.80 (2019). For the EVI, the ranges are as follows: from 0.0 to 0.41 (1984), −0.08 to 0.56 (2002), −0.05 to 0.74 (2013) and −0.04 to 0.60 (2019). In Lagos State, many patches of bare land are sandwiched between built-up areas. Hence, this mixed set-up could explain the little differentiated values between the NDVIs of bare land and built-up area. It has been discovered that the values of EVI are generally lower than the NDVI values. This is due to the atmospheric correction and reduction of soil reflectance influence, which has been factored in the EVI model. Figure 3 is the histogram representation of the range of values of LST, NDVI and EVI, respectively, from 1984 to 2019. The histogram shows the dominance of LST values within the range of 23–24°C, the NDVI values within the ranges, 0.35–0.42 and 0.58–0.62 dominate the Lagos Lagoon environment over the years, the EVI values within the range of 0.30–0.40 dominate the Lagos Lagoon environment over the years. These values of VI show the dominance of greenness in the Lagos Lagoon environment despite the urbanisation. Figures 4–7 present maps of land cover, LST, NDVI and EVI for the periods under study. The progressive increase in built-up areas and the decrease in vegetation cover are evident in Figure 4. The LST map also shows a progressive rise in temperatures emanating from the centre of the metropolis. Juxtaposing the four maps, it can be deduced that the lagoon environment has varied vegetation where vegetation decline over the 35 years follows a pattern of increasing disappearance from the west to the east. This also indicates that urban development is very high in the west and it increases yearly towards the north-western and eastern regions of the lagoon environment. Generally, the vegetation-covered areas have the highest NDVI and EVI values. Areas with abundant vegetation cover display low LST values (Guha et al. 2020). The trend of overall LST ranging from higher to lower for built-up areas, bare land and vegetation conforms with the findings of Ullah et al. (2019).

Fig. 3

Fig. 4

Fig. 5

Fig. 6

Fig. 7

Figures 8–10 present graphical illustrations of the mean LSTs, NDVIs and EVIs associated with land cover types in the eight LGAs surrounding the Lagos Lagoon with error bars representing ±1 SD. The eight LGAs (Lagos Island, Somolu, Ibeju Lekki, Ikorodu, Somolu, Epe, Eti-Osa and Kosofe) with their spatial locations show varying degrees of urbanisation. The highest mean LST values are observed in Lagos Island, Lagos Mainland and Somolu LGAs and this is possible because of the high population growth rate and urbanisation, which invariably causes a high degree of human-induced heat. Hence, the increase in urbanisation must have resulted to high mean values of LST in each of the years of study. This is one of the effects of UHI. Notwithstanding, areas such as the extreme eastern part of Ikorodu, Ibeju-Lekki and Epe exhibit the highest means of NDVI and EVI values (Figs 9 and 10). This might be attributed to the predominance of vegetation cover. In agreement with the findings of Wilson et al. (2003) and Yue et al. (2007), green vegetation cover in the areas suggests relatively a higher rates of evapotranspiration and favouring of latent exchange between surface and atmosphere as compared with built-up areas. In the study by Ferrelli et al. (2018) on the Spatio-temporal relationship between LST and NDVI in Monte Hermoso City, Argentina, minimum NDVI was observed in summer whereas the LST showed the opposite behaviour. Similarly, in the present study which was conducted with Landsat imageries acquired in the dry season, it is observed that mean LSTs are in opposite trend with the mean NDVIs and EVIs. Also, where high LSTs are observed, NDVI diminishes due to variations in vegetation state or abundance.

Fig. 8

Fig. 9

Fig. 10

Table 9 presents the coefficients of correlation derived from the analysis. In all the pairs of comparisons, levels of correlation were detected at the 0.01 level (2-tailed). Generally, a strong inverse relationship was observed between the NDVI and LST and between EVI and LST. An inverse relationship between NDVI and LST was also noted in studies by Li et al. (2010), Ferrelli et al. (2018), Mukherjee and Singh (2020) and Guha et al. (2020). The strongest positive correlations were observed between NDVI and EVI over the 35-year data span. This means that the LST increased as the NDVI and EVI decreased over time. However, the value of the coefficient of correlation (r) decreases from 1984 to 2019 (r values shown in bold). This implies that although high positive r was obtained, the VI around the lagoon environment is decreasing yearly, inferred from the gentle gradient of r from 1984 to 2019.

The regression analysis between the parameters and the time period presented in Table 10 shows a negative linear relationship between the LST and the VIs over the years. Several authors have also noted this inverse relationship, for example, Dewan and Corner (2012), Ferrelli et al. (2015), Guha et al. (2020) and Mukherjee and Singh (2020). This means that the LST increased as the NDVI and EVI decreased over time. The linear models show the existence of significant relationships between the VIs and the LST. In consonance with the submission of Ferreira and Duarte (2019), the R² shows that the NDVI has a better fit and relationship with LST. The relationship between the LST, NDVI and EVI in vegetation cover shows a better fit than the bare land and built-up areas with R² values as follows: vegetation: (0.729, 0.770), bare land (0.542, 0.557) and built-up area (0.410, 0.421). This shows the homogeneity in the relationship between the LST and the VIs in the vegetation class. Figure 11 shows the graphical representation of the relationship between the LST, NDVI and EVI across the years. Figure 11a shows the heterogeneous distribution of the LST over the years with respect to the VIs (NDVI and EVI) in bare land. Between the years 2014 and 2019, high LST values 27.5–28.5°C have moderate NDVI and EVI values between −0.2 and 0.6. Figure 11b shows a less heterogeneous distribution of the LST over the years with respect to the VIs (NDVI and EVI) in built-up area. Between the years 2014 and 2019, high LST values (29–30°C) have moderate NDVI and EVI values between −0.2 and 0.2. Figure 11c shows a homogenous distribution of the LST over the years with respect to the VIs (NDVI and EVI) in vegetation. Between the years 2014 and 2019, high LST values (29–30°C) have moderate NDVI and EVI values between −0.2 and 0.2.

Fig. 11

Table 11 presents the CI calculated for bare land, built-up area and vegetation. The analysis suggests that the increasing LSTs are largely attributable to the built-up areas due to the highest CIs of 0.08 in 1984, 0.14 in 2002, 0.67 in 2013 and 0.35 in 2019. Conversely, vegetation had a cooling influence and negative contribution to LST as evident in the CIs of −0.10 in 1984, −0.19 in 2002, −0.68 in 2013 and −0.36 in 2019. The cooling influence of the vegetation is increasing over time while the warming influence of built-up areas is also increasing.

In 2013, vegetation recorded its highest cooling influence. However, between 1984 and 2002 there was no significant cooling influence of vegetation on the LST. This is confirmed by the depletion of vegetation in the lagoon environment within the same period. Remarkably, while the CI value for vegetation was declining from 2003, that of built-up areas was increasing (Fig. 12). This indicates that the cooling influence of vegetation reached its maximum in 2013 before it started declining. However, that of the built-up area followed the opposite trend. This implies that the cooling effect of the vegetation was very high between 2003 and 2013 in the vegetation-covered areas but very low in the built-up areas. Possibly this could be a result of the tree greening programmes enforced by the Lagos State Government to create more green areas by planting trees and expanding the green environment in Lagos State generally (Adegboye 2013, Soladoye, Oromakinde 2013). The 2013 tree planting campaign was the fifth anniversary of tree planting initiated by the state government since 2008. A possible explanation for the reduction in built-up area CI between 2013 and 2019 is the influence of sea breeze which conditions the climates of coastal urban areas, leading to incidence of low temperatures in coastal areas (Ferrelli et al. 2018). Interestingly, a point of intersection was observed (Fig. 12) in 2005 with CI value of approximately −0.3 and 0.2 for vegetation and built-up area, respectively. This proposes a critical point of decision because in 2019 the same value is almost becoming the intersection point. This needs to be further investigated to initiate a decision plan for positive action in the lagoon environment.

Fig. 12

Conversely, the relationship between the CI value of vegetation and bare land shows that bare land had almost no significant contribution to the cooling effect in the study area (Fig. 13). However, in 2003 while the cooling effect by the vegetation began to increase, that of bare land was decreasing very slowly until it reached a peak of zero in 2019.

Fig. 13

Table 12 presents the mean LST, NDVI and EVI across the years, which were variables in the trend analysis. At 90% confidence interval, i.e. p-value = 0.10, there is a significant upward trend in the LST values over the years across the study area with the Mann-Kendall statistic (S) of 6 and significance value lower than 0.10. Further trend analysis was carried out for individual land cover classes. The parameters LST, NDVI and EVI appear to have an upward trend in the bare land class with S value of 6 and significance value 0.042. The LST, NDVI and EVI do not follow an upward trend across the years for built-up areas. The LST follows an upward trend across the years in the vegetated area but not the NDVI and EVI. Sen's estimate of slope (Q_med) shows that the mean LST increases by 0.178°C annually in the Lagos Lagoon environment. When considered on a land cover class basis, the bare land mean LST increases by 0.115°C a⁻¹, the mean NDVI and mean EVI by 0.006 a⁻¹ and 0.003 a⁻¹, respectively.