Monitoring and forecasting spatio-temporal LULC for Akure rainforest habitat in Nigeria

The producer georeferenced Landsat 7 ETM+ imagery for the year 2009 was downloaded with strip lines. The strip lines are as a result of the failed scan line corrector attributed to Landsat 7 ETM+ datasets obtained after 31 May 2003. The Landsat (LE71900552009326ASN00) data attribute indicated a scan gap interpolation of two for the 2009 imagery. A total of 28 strip lines were identified, accounting for 11.7% of the demarcated study area. The strip lines were gap-filled using the ArcGIS 10.3 software. This was achieved with the Landsat toolbox ‘fix Landsat 7 scan line error’. The procedure was performed for each band 1–8 before the atmospheric correction procedure.

The atmospheric correction for the Landsat dataset was performed using the following steps. The digital number (DN) of each band was first converted to radiance at the level of top-of-atmosphere (TOA) radiance. The TOA radiance was then converted to TOA reflectance. The procedure was achieved using the raster calculator in the ArcGIS environment. The relevant metadata required for the procedure were obtained from the downloaded MTL text file described by Chavez (1996). For possible geometric distortions, the Landsat scenes utilised for this study were within the prescribed image-to-image tolerances of ≤12 m geometric root mean square error (RMSE) (Table 1).

Supervised classification was performed by developing a spectral signature of the verified land-use classes based on the Anderson (1976) classification regime. An accuracy assessment was performed to evaluate the conformity between the benchmark expected to be correct and the classified image pixels of unknown quality. Landsat 7 ETM+ and Landsat 8 OLI data were used to derive the LULC types for the study period using the Erdas Imagine 9.2 software. Before the supervised image classification, a false colour composite (FCC) analysis was performed to augment the ground-truth LULC information. This was achieved using the bands 4, 3, 2 in Landsat 7 ETM+ and bands 5, 4, 3 in Landsat 8 OLI. The five classes utilised are as follows: bare land, built-up area, cultivated land, forest cover and grassland. Bareland is any area with exposed land surface or rock outcrops. Built-up area is any area with asphalt or concrete roads. Cultivated land refers to an area under cultivation and intensive irrigation of crops pavements, buildings or houses. Forest cover is an area of dense tree cover with a thick closed canopy. Grassland is any non-vegetated, uncultivated farmland and open space. The rule-based feature extraction method combined with ground-truthed information was adopted for this purpose. The MLA was adopted for the supervised classification. It adopts the Bayesian decision rule based on the probability that a pixel belongs to a particular class. It is one of the most adopted and validated methods of classification in remote sensing (Hossen et al., 2018; Sahana et al., 2018).

The accuracy assessment is performed to evaluate the degree of misclassification among LULCs to the user classified image from collected ground (training) data with the producer information. A total of 181 training data samples were obtained as reference samples over the study area for the individual Landsat images. Two methods of assessment were utilised to evaluate the accuracy of the classification procedure; the error (confusion) matrix and Kappa statistics (Table 2).

Change detection entails a timely and accurate assessment of the changes in the LULC variables at periodic intervals. In this study, change detection analysis was performed for the five LULC classes at intervals, 1999–2009, 2009–2019 and 1999–2019. Three change maps were produced to display the change of the classified images. The change detection statistics were performed using the IDRISI software.

The LAC and LCR was determined using equations (1) and (2) (Yeates and Garner, 1976). To compute the LCR/LAC, the following variables are required. They include the land area of the city (built-up) in hectares and the corresponding population size, for the respective years of interest. The estimated built area of the city for the epoch of interest was obtained using the classification distribution (Table 3). (1) LCR=CityExtentinHectares(A)Population(P), {\rm{LCR}} = {{{\rm{CityExtentinHectares}}(A)} \over {{\rm{Population}}(P)}}, (2) LAC=A2−A1P2−P1, {\rm{LAC}} = {{{A_2} - {A_1}} \over {{P_2} - {P_1}}}, where A₁ is the extent (in hectares) for the early year, A₂ is the extent (in hectares) for the later year, P₁ is the population figure for the early year; P₂ is the population figure for the later year.

The population size was estimated (Parker, 2002) using the projected population equation (equation (3)) with the establishment census data of 2006 and the 3.0% growth rate for the study area.

The built-up area class is the key determinant for evaluating urban area growth. To formulate future urban area growth for the study area, the CA Markov model was adopted to capture the inherent spatio-temporal relationship based on the classified Landsat images. The CA Markov model analysis adopts the preceding LULC to project future spatial distribution (Liping et al., 2018).

The CA Markov model is one of the commonly used models used for monitoring the stability of land-use developments via future prediction. The model integrates CA with the Markov chain to predict the LULC trends and characteristics over time (Hamad et al., 2018). Equations (4) and (5) illustrates the computation of predicted land use (Kumar et al., 2014). (4) s(t,t+1)=Aij×s(t) s(t,t + 1) = {{\bf{A}}_{{\bf{ij}}}} \times s(t) where S(t) is the system status at the time t, S(t+1) is the system status at the time of t+1; A_ij is the transition probability matrix computed using equation (5). (5) ‖Aij‖=[A1,1A1,2…A1,NA2,1A2,2….A2,N…….….….AN,1AN,2….AN,N] \left\| {{{\bf{A}}_{{\bf{ij}}}}} \right\| = \left[ {\matrix{ {{A_{1,1}}} & {{A_{1,2}}} & \ldots & {{A_{1,{\rm{N}}}}} \cr {{A_{2,1}}} & {{A_{2,2}}} & { \ldots .} & {{A_{2,{\rm{N}}}}} \cr \ldots & { \ldots .} & { \ldots .} & { \ldots .} \cr {{A_{{\rm{N}},1}}} & {{A_{{\rm{N}},2}}} & { \ldots .} & {{A_{{\rm{N}},{\rm{N}}}}} \cr } } \right] where 0 ≤ A_ij ≤ 1 determines the probability of a pixel changing to another LULC or maintaining its original LULC.

The CA Markov module embedded on the IDRISI Software (Sang et al., 2011) was utilised to generate the LULC prediction map for the year 2039.

The gains/losses, net changes and contributors are presented in Figure 4. In Figures 4a–4c, the green bar indicates the gain per class in km², while the loss (to the left) of each class is displayed the purple.

Figure 4

In Figure 4a, it can be observed that for the period 1999–2009, bare land class lost 329.63 km² and gained 196.29 km², built-up area class entirely gained 143.73 km², cultivated land class lost 1.24 km² and gained 12.13 km², forest cover class lost 62.29 km² and gained 94.87 km² and grassland cover class lost 238.49 km² and gained 184.64 km². In Figure 4b, the gain and loss of LULC from 2009 to 2019 we derived to be bare land lost 351.15 km² and gained 58.47 km², the built-up area lost 19.95 km² and gained 345.74 km², cultivated land entirely gained 89.65 km², forest cover lost 103.59 km² and gained 112.21 km² and grassland cover lost 265.66 km² and gained 134.28 km². The loss and gain from 1999–2019 (Figure 4c) showed that bare land lost 474.21 km² and gained 48.19 km², the built-up area gained 469.51 km², cultivated land entirely gained 100.53 km², forest cover gained 167.54 km² but lost 126.34 km². Grassland cover lost 278.15 km² and gained 92.93 km². The steady increase in the forest cover can be attributed to climatic variables as well as the population's cultural belief against forest intrusion.

Figures 4d–4f illustrates the net changes in terms of loss (negative) and gain (positive) along the x-axis for the five LULCs of interest between 1999 and 2009 (Figure 4d), 2009 and 2019 (Figure 4e) and 1999 and 2019 (Figure 4f). For the respective intervals, ‘bare land’ recorded −133.34 km², −292.68 km² and −426.02 km²; ‘built-up area’ recorded +143.73 km², +325.79 km² and +469.51 km²; ‘cultivated land’ recorded +10.89 km², +89.65 km² and +100.53 km²; ‘forest cover’ recorded +32.58 km², +8.62 km² and +41.20km² and ‘grassland’ recorded −53.85 km², −131.38 km² and −185.22 km². Figures 4g–4i identified bare land and grassland as the main contributors to the consistent increase in the built-up area class across the study intervals.

The projected population estimates for the year 1999, 2009 and 2019 were computed as 394,176, 529,752 and 711,942, respectively (National Population Commission, 2010). Adopting the projected population size for the LCR/LAC analysis, it can be observed that the LCRs for the year 1999, 2009 and 2019 were computed (in hectares) to be 0.00505, 0.00362 and 0.0687, respectively. These values revealed there is a steady rise of the land consumed as a result of projected population growth. The LAC was also determined. With equation (2), the computed result revealed an LAC of 0.1635 from the years 2009 to 2019. Also, the absorption coefficient for 1999–2019 was computed to be 0.127. These LAC values obtained were higher than the coefficients reported by Salghuna et al. (2018) for Andhra Pradesh in India.

Before the prediction, a genetic optimisation from the year 1999 to the year 2019 was performed. The 2009 image was used to validate the prediction of the CA Markov by comparing the predicted LULC in 2009 and the actual supervised classification for the year 2009. The overall Kappa coefficient was derived to be 0.864. This value was within the acceptable limits (Kappa index of agreement); thus, the prediction simulation was embarked upon. The CA Markovian analysis predicted the changes that are likely to occur over the next two decades (2039). It projected that about 675.011km² of the study will be built-up area (Figure 5). The predicted built-up area increase (+37.92%) for this study was within limit when compared with the +27.31% increase for Damman coastal city, Saudi Arabia (Rahman et al., 2017). By contrast, the predicted built-up area increase was far less than the +111.34% for the Tarai region, Nepal (Rimal et al., 2020), and the +300% for Dhaka, Bangladesh (Ahmed et al., 2013).

Figure 5

From Figure 5, it can be deduced that the rising uncoordinated urban growth in Akure is very likely to translate to declined biodiversity and increased urban air pollution from increased, unregulated anthropogenic activities. This is, therefore, likely to expose the residents to respiratory health challenges as well as directly influence the city's urban heat island (UHI) analysis.