Time Series Analysis of Landsat Images for Monitoring Flooded Areas in the Inner Niger Delta, Mali

The analysis of the environmental dynamics using remote sensing data has received considerable attention in Earth sciences during the last decades (Mertikas et al., 2021; Payra et al., 2023; Sultana et al., 2023). Satellite images are widely recognised as a valuable source of spatial information (Costel and Bariou, 1992; Lerat, 1987; Richards, 2013). To process these data effectively and accurately, a variety of algorithms for processing satellite images have been developed in remote sensing domain. A general approach of using multi-temporal satellite images for environmental monitoring is based on their pairwise comparison aimed at evaluating changes in land cover patterns (Hiernaux et al., 2021; Ogilvie et al., 2015). Time series analysis of Earth observation data has proven to be effective in the evaluation of landscape changes using several images covering the same area in various consecutive years. For instance, satellite-derived trends are used as monitoring methods in a wide variety of environmental applications: mapping and monitoring wetlands (Wu, 2018), assessment of deforestation and forest degradation (Haarpaintner and Hindberg, 2019; Mashhadi and Alganci, 2022; Masolele et al., 2021; Schneibel et al., 2017), monitoring wetland dynamics (Kovács et al., 2022; Xie et al., 2022), evaluation of vegetation cover fraction and soil depletion (Dube et al., 2017; Gallo et al., 2023), computing vegetation indices (Lemenkova and Debeir, 2022a; Liu et al., 2022; Venter et al., 2020), estimating variations in land surface temperature (Carrillo-Niquete et al., 2022), assessment of spatio-temporal variations in night lights emissions in urban studies (Rehman et al., 2021) and more.

Such case studies assume that evaluating land cover changes using time series analysis of satellite images is possible by a comparison of images taken with repetitive time gap of years during the same months. This is based on finding the differences in spectral reflectance of land cover types in various bands of a satellite image for the same spatial extent and processing images using various classification techniques (Gandhi and Sarkar, 2016; Merry et al., 2023; Shahi et al., 2023). Since spatial object depicted on the image scenes has different brightness reflected in distinct colours of pixels, a proper combination of spectral band followed by image classification enables to highlight target features that categorise land surface objects. Hence, one can compare the shifts in spatial geometry of the land patterns on the Earth using a sequence of satellite images taken in different years for the same region.

This assumption, however, requires effective methods of image processing for analysis of a time series rapidly and accurately. While many existing GIS software involve a variety of manual operations for processing each image (Epuh et al., 2022; Spiekermann et al., 2015; Tappan and McGahuey, 2007), the success in classification of the satellite images requires the development of more advanced approaches. The tasks of satellite image processing require machine-based programming algorithms for accurate classification and mapping. Thus, remote sensing data processing and analysis for mapping vegetation indices involve several steps of workflow including the analysis and composites of spectral bands, computing the indices using algorithms, originating from existing formulae and visualising and mapping the results, which makes it a cumbersome task by traditional GIS software. In contrast, automated workflow, which involves programming languages such as R, recognises specific features of images using machine-based algorithms and scripts and computer vision algorithms of image processing that can classify the images using spectral reflectance data or compute the vegetation indices accurately and rapidly where manual image processing takes more efforts and time.

We believe that the traditional GIS approaches do not meet the requirements of modern applications of image processing in time series. For example, the efficiency of time series analysis based on rapid and accurate processing of the image set is a fundamental requirement to be able to perform comparative analysis of the land cover types by classification of multiple scenes taken on different years. The satellite image collections are growing at significant rates. For instance, only the Landsat-8 and Landsat-9 are collecting daily 1,500 new scenes archived into the USGS public repositories (U.S. Geological Survey, 2015). Together with other satellites (Sentinel by ESA, MODIS by NOAA) including commercial ones (SPOT), this creates a precious pool of spatial information (Maini and Agrawal, 2006). The processing of such big data requires automated tools for rapid and effective handling to convert them from the raw source of information into practical knowledge.

In contrast, programming methods of image processing use machine-based algorithms for detecting object features and pattern recognition (Barma et al., 2020; Lemenkova and Debeir, 2022c; Lighezzolo et al., 2019; Mosca et al., 2020). In such cases, the automated discrimination of land cover types is supported by the embedded algorithms of image processing with minimised control and supervision from the user. Furthermore, while a set of programming libraries can be used to process the lines of the prewritten code in general tasks, the interactive nature of image processing and the algorithms of classification call for the ability to allow users to define land cover classes (Liu and Mason, 2009). Thus, besides the statistical analysis of R applied in various research and reported in relevant studies (Bivand, 2000; Grunsky, 2002; Lemenkova, 2019), it can also be effectively used in mapping for discriminating land cover types.

Land cover types automatically discriminated on the satellite images can be divided into diverse categories to be recognised and retrieved from the images using real situation of the local vegetation patterns. For instance, vegetation associations can include various types of forests, grasslands and savannah, wetlands and floodplains, bare soil, urban lands, cultivated agricultural spots, built-up areas and artificial objects (buildings, roads, etc.) (Lillesand et al., 1994). The automated approach classifying pixels for object interpretation on the satellite images uses the programming methods. This is based on the machine-based discrimination of similar spectral characteristics and assigning them to the target land use classes.

The k-means clustering technique was used for unsupervised image classification as a commonly accepted method for remote sensing data processing (Bovolo et al., 2018; Esche and Franklin, 2002; Hou et al., 2016; Paola Patricia et al., 2020). The k-means clustering in image processing by R presents a method of image partition that divides the raster matrix with n pixels (or cells of this matrix) into k groups. Here, each pixel of the image scene belongs to one cluster based on the value of its spectral reflectance. The criterium of this partition is the nearest mean where the cluster centroids present the centres of the groups. More specifically, the initial class means are the two important class parameters of the pixels on the satellite images that were defined using the characteristics of spectral reflectance of the recording of the pixels (i.e., cells of the image matrix). Such algorithm is implemented via the image analysis techniques in unsupervised classification by R. Spectral reflectance is a common parameter of the image’s pixels which represent the earth surface features.

The environmental vegetation studies use remote sensing data to analyse the phenology of plants by evaluating information collected from the satellite sensors (Karkauskaite et al., 2017; Schucknecht et al., 2017). Measured strength of light absorbed and reflected by green leaves of vegetation gives the information regarding the health and greenness of canopy. Such algorithm identified spectral reflectance and values of Digital Numbers (DN) of pixels by creating new data for detecting higher and lower values in vegetation indices which point at the amount of chlorophyll related to health and growth conditions of the plants (Alcaras and Parente, 2023; Lemenkova and Debeir, 2023a). Thus, the data gathered by satellite sensors measure strength of reflected light that indicates vegetation health through values in Red and NIR bands.

The core approach of the VI analysis consists in the properties of pigments in plant leaves which strongly absorb wavelengths of Red visible light and reflect those of near-infrared (NIR) light (Garzonio et al., 2017). Therefore, using a combination of Red/NIR band values gives information on the state and condition of health in leaves in plants. Detected wavelengths of light absorbed and reflected by green plants split the image scene through regions of VI values with high and lower values around canopy areas (Lemenkova and Debeir, 2023b; Marino and Alvino, 2019). The areas of the vegetation coverage are then defined and used to evaluate the level of greenness to discriminate the healthy plants on the images against the wilting and diseased plants. This process is based on colour values of pixels to identify vegetation vigour for canopies breaking the images into areas of healthy or diseased vegetation and bare land.

After the values of the VI are calculated for each pixel on the image scene, the image is converted to a map showing the distribution of these values across the scene with values of pixels evaluated for the VI using the algebra of band combination for each index, respectively. By tracking the image for these values, the segments of various VI values are traced and marked with assigned values of VI as parts constituting the matrix of pixels on each image scene. This resulted in a matrix of pixels representing the level of vegetation healthiness and greenness (Ma et al., 2019). For the NDVI index, for instance, higher values signify healthy vegetation, lower values – sparse of decease vegetation, while negative values – no vegetation, water areas, etc. The interpretation of the VI values may slightly differ in various indices.

In this paper, we present a case study of remote sensing data assessment by R language for environmental monitoring of Inner Niger Delta, Mali, West Africa (Figure 1). As largest inland wetland in West Africa, the habitats of Inner Niger Delta offers high potential for biodiversity of the flood-dependent ecosystems with endemic woody plants and endangered species. Moreover, rich water and land resources of the Inner Niger Delta create high agricultural potential and favourable conditions for irrigated agriculture and rice crop. As a result, the Inner Niger Delta is one of the most productive areas of rice in West Africa and environmentally vulnerable region of Mali.

Figure 1.

Despite the environmental importance of this region, there are relatively few published papers on the analysis of the land cover dynamics in Inner Niger Delta, Mali, using Earth observation data. Existing studies have focused on the meteorological observations from the gauge stations, fieldwork measurements and locally obtained agro-hydrological information (e.g. Dadson et al. (2010)) assessing ecosystem services (Rebelo et al., 2013), simulations of datasets by statistical analysis (Mascaro et al., 2015), hydrodynamic modelling based on streamflow data and water levels to estimate the inundation extent (Haque et al., 2021, 2020). In contrast, obtaining information about the object through the analysis of its characteristics acquired from the space-borne data collected at a distance present a more accurate, advanced and reliable method (Bergé-Nguyen and Crétaux, 2015; Maini and Agrawal, 2010).

To contribute to the environmental assessment of this geographically unique region by the advanced methods, we provide a general machine-based solution of satellite image processing. The objective is the evaluation of land cover dynamics using six Landsat 8–9 images of the surroundings of Mopti and Djenne, central Mali. Each image was taken on November from 2013 to 2022. The evaluation of the land cover types and vegetation patterns using remote sensing data processed by the programming methods enables to get better insights into landscape dynamics resulted from the complex land-atmosphere interactions.

Our approach is based on using libraries of R programming language for analysis of spatial and temporal structures of the satellite images on 2013, 2015, 2018, 2020, 2021 and 2022. Using algorithms of R, we processed the remote sensing data to yield estimates of landscape dynamics and visualise changes in major land cover types of the Inner Niger Delta that are temporally coherent. Further, we analysed fluctuations of spatial extent in fluvial wetlands caused by the hydrological processes of seasonal flooding. Finally, by allowing packages of R to support image processing, an approach to mapping several vegetation indices and visualising changes in distribution of different land cover classes over time is realised. We tested R approaches on three computed vegetation indices which included the Normalised Difference Vegetation Index (NDVI), the Soil Adjusted Vegetation Index (SAVI) and the Enhanced Vegetation Index 2 (EVI2).

Mali, a landlocked sub-Saharan West African country (Figure 1), is one of the hottest countries in the world with an average temperature varying between 24°C in January and 35°C in May (for Development Practitioners and Makers, 2022).

Subject to recent climate change, it is characterised by a significant variability in rainfalls, droughts and occasional floods. Rising temperatures, interannual to decadal variability in rainfalls and decreased precipitation lead to the desertification and deteriorated plant communities in the ecosystems of Mali. In turn, the decline of vegetation that earlier retained water in the soil aggravates soil dryness and depletion contributing to the desertification. For arid climate of Mali, such a vicious circle of interrelated processes has disastrous consequences on the vegetation growth and as a consequence, on food production and security (Davies, 1996). Since recent decades, the climate setting in West Africa changed with the most notable trends in decrease of precipitation and rainfall, rise of temperature and warming trends, increase in drought frequency and dryland expansion (Leal Filho et al., 2022; Oguntunde et al., 2017, 2006; Sidibe et al., 2019). Many papers are published on climate change issues and reported consequences for biodiversity and socio-economic systems (Huang et al., 2016; Le Houérou, 1996; Yu et al., 2017).

These processes triggered the desertification, escalated the salinity of soil and effects on ecosystems, e.g., loss of biodiversity (Heubes et al., 2013), changed geographic distribution of species (Coulibaly et al., 2023) and declined vegetation health and tree growth (Sanogo et al., 2022). Further effects include the decline in agricultural yields that are highly dependent on regional meteorological conditions and rainfall variability (Zare et al., 2017). For instance, the decrease in precipitation negatively affects rainfed crop yields in sub-Saharan Africa (Bambio et al., 2022; Raes et al., 2021) with reported cases on rice (Adjah et al., 2022; Akpoti et al., 2022), millet (Bado et al., 2022), maize (Srivastava et al., 2016) and sorghum (Adam et al., 2020).

In this study, we focus on monitoring the dynamics of the flooded flat plain savannah areas in the Inner Niger Delta located in central Mali (13°30′ N–15°30′N). Situated in the midst of the Malian Sahel, this seasonal floodplain comprised multiple meandering channels of Niger River, the third largest river in Africa. The Inner Niger Delta is one of the larger wetlands in Africa with 200 km of spatial extent and area of 36, 000km² (Kuper et al., 2003). The vegetational pattern of Inner Niger Delta is interspersed by the seasonally flooded savannah, plateaus with the desiccated sands and wooded areas in the north-Sudanian zone. The geomorphology of the floodplain in Inner Niger Delta presents almost flat relief with variations from 272 to 262 m (Gallais, 2003).

The alluvial depositions result in a series of micro deltas where the hydrography of Niger divides forming a local basin and then regroups in the downstream in a confluence section of the river (Dickens et al., 2018). The geologic setting of the Inner Niger Delta includes a wide Precambrian basement occupied by primary sandstones where the interior delta is formed. A quaternary deposits, accumulated by the winds from the Sahara, added the effect of a dune dam to the generally flat relief. The Niger River meanders in the conditions of the flat low-relief plains and divides into several branches of small channels, whose alluvial ridges enclose more depressed areas (Gourou, 1969).

High poverty of Malian population increases socio-economic vulnerability and low adaptive capacity towards natural hazards in changing hydro-climatic conditions. Specifically, the environmental vulnerability of Mali to climate change – one of the world’s poorest countries – includes factors such as hydroelectricity production and dam operations (Haque et al., 2019), navigation over Niger, seasonality of the irrigated and rainfed agriculture during dry and rainy periods and flooding in the Inner Niger Delta (Ghile et al., 2014). As a result, climate extremes contribute to low sustainability of food resources depending on fishing and agricultural farming that are highly sensitive to seasonality and weather constraints (Liersch et al., 2019; Morand et al., 2012). Changes in selected crop production patterns have been reported recently. For instance, Laris et al. (2015) reported the decline in cotton yields in Mali due to rapid intensification of agriculture or soil fertility loss, while rise in maize yields as less sensible towards soil conditions and less demanding crop in farming system.

The Malian climate has three distinct seasons: a dry season from March to June, a rainy season from June to September and a cold windy season from October to February. Such fluctuations of the sahelian climate correspond to the regular seasonal flooding of the Inner Niger Delta which results in a unique hydrologic phenomenon – seasonal floods (Brunet-Moret et al., 1986). The Inner Niger Delta becomes inundated during several months a year with floodplains and remaining a few islets of dry land (Courel and Chamard, 1994). Previous studies report (Diakite, 1986) the area of 20, 000km² regularly covered by the flood. For the remaining six months, the floodplain dries out again to a fertile plain with essentially grass vegetation, e.g., bourgou (Planhol, 1971). The repetitive regular pattern of climatic-hydrological shifts affects the distribution and geochemistry of soils and vegetal associations (Tabeaud, 1980). Water-rich alluvial soils, together with arid climate and high temperatures create perfect conditions for plantations (e.g., rice harvest) as a precious resource for Malian fishermen, pastoralists and farmers (Retaillé, 1984). However, harvest is at risk of lost during flash heavy rains or floods.

At the same time, the variability of inflow into the Inner Niger Delta, the period of the floods and the extent of the regularly flooded surface area have direct consequences on livelihoods, and food production (Liersch et al., 2013). Repetitive droughts on the one hand and dynamics of repetitive floods on the other hand strongly affect the ecology and biodiversity of Inner Niger Delta with direct consequences on vegetation system (Mariko, 2003). Besides, shifts in precipitation and temperature patterns in Inner Niger Delta result in depletion of the irrigated soils and salinisation as a consequence. Apart from the distribution of saline soils (solonchaks), this region is under pressure from the intensive irrigation system which uses poorly mineralised surface water. The decrease in precipitation, increase in annual temperature and high evapotranspiration amplify the climate-related effects and result in the deposition of the dissolved salts on the surface (Valenza et al., 2000).

In turn, small-scale farming systems create additional environmental challenges by overgrazing and uncontrolled livestock contributing from their part to land cover changes. Besides livelihood and food-related issues, stagnant water during floods contributes to health risks of the diseases of population through rapid distribution of parasites and insects, such as malaria (Cools et al., 2013; Klinkenberg et al., 2002). All these issues require to take measures on environmental monitoring of floodplain wetlands in Inner Niger Delta. Previous studies pointed the need for cartographic visualisation of the unique territories of the African continent which should be supported by maps for quantitative investigations (Cabot, 1968), e.g., in phytogeographical studies (Robequain, 1948). Insofar as we can map hydrological processes and their relationship with vegetation growth, we can get deeper insights into the complex links of climate-environmental interactions in Africa.

The time series analysis is based on the Landsat 8–9 satellite images collected from the USGS survey (Department of the Interior U.S. Geological Survey, 2022) on the surroundings of Mopti and Djenne, Mali. The metadata of the scenes are summarised in Table 1. The images were captured from the EarthExplorer repository (Figure 3). We used a series of the six Landsat 8–9 images to map the areas of flooded vegetation within the area of Inner Niger Delta during the years 2013, 2015, 2018, 2020, 2021 and 2022.

Figure 3.

The images were collected always on November for the following reason. Due to the seasonality of floods and repetitive cycle of inundation, a strong contrast in moisture and green vegetation reaches its peak in autumn period of intense rains, that is, from late October to mid January which enables to compare the contrasting flooded areas against the dry lands. While the period of rainfalls ends in the Inner Niger Delta by late October and soils and grassy vegetation become dry in the regions untouched by floods, the southern segment of the Inner Niger Delta demonstrates the highest levels of flooded areas. Such contrast in the northern, central and southern segments of the Inner Niger Delta favours visualisation of the flooded/not flooded areas and estimation of the intensity of floods in this particular year. The second target requirement was set up for all the imagery regarding the least possible cloudiness which is below 5% for all the scenes.

Therefore, we have collected the images on November for each evaluated year with the minimal cloud coverage in technical characteristics of Landsat. The characteristics of the scenes are summarised in Table 1. The map in Figure 1 has been prepared using the Generic Mapping Tools (Wessel et al., 2019) by applied existing cartographic workflow reported earlier (Lemenkova and Debeir, 2022b,d).

The cartographic technical characteristics common for each Landsat band include the coordinate reference system World Geodetic System (WGS) 84, the Universal Transverse Mercator (UTM) zone 30N (EPSG:32630) and the spatial resolution of most of the bands at 30 meters. Specifically, the OLI multispectral bands have pixel size 30 m. These include Bands 1 to 7: Band 1 (Coastal aerosol), Band 2 (Blue), Band 3 (Green), Band 4 (Red), Band 5 - Near Infrared (NIR), Band 6 (SWIR 1), Band 7 (SWIR 2) and Band 9 (Cirrus). The panchromatic Band 8 and Thermal Infrared (TIRS) 1 and 2 bands were not used in this study. Here, we used the multispectral bands which allow the precision of the maps to 30 m per pixel. The six Landsat 8–9 OLI images used in this study are visualised in RGB colours in Figure 4.

Figure 4.

The principal characteristics of the Landsat 8–9 images consist in multi-spectral bands which enable to detect difference in reflectance characteristics of various land cover types by image analysis. This is explained by different ratio in energy reflected and incident by these surfaces and recorded in wavelengths. Therefore, various land cover types (water, land, diverse vegetation types, asphalt-covered roads or buildings) have different spectral reflectance (Chapman and Gasparovic, 2022; Richards and Jia, 2006). Spectral reflectance is measured as a function of wavelength in the Landsat 8–9 bands. Such spectral characteristics enable to separate these pixels into classes on the images.

The concept of the entire process of image processing is shown in Figure 5. It summarises the main steps of the two separated processes – the calculation of the three vegetation indices (NDVI, SAVI and EVI) and the classification score using k-means clustering approach based on Landsat 8–9 images, followed by accuracy assessment. The flowchart was plotted using R library DiagrammeR. The dataset included the 6 Landsat images taken on various years. The selected bands were imported in R environment and processed using libraries. Two different separated steps were applied in this work – computing the vegetation indices (VI) and classification of land cover types in the Inner Niger Delta, Mali.

Figure 5.

For computing the VI, three different indices were calculated – NDVI, SAVI and EVI using equations for band combinations for each case. The selection of indicators is based on the diapason of the vegetation indices which ranges from −1 to 1 for NDVI, SAVI and EVI. After running the script in R, the results were obtained including maps and histograms of data distribution. Next, we performed clustering of the images based on the k-means algorithm. The results of clustering and computed VI were visualised on the maps. To confirm the advantages of the methods and evaluate the performance of the classification models, we assessed the accuracy by Kendall confusion matrix.

In this section, we briefly describe our GMT scripting method used to plot the topographic map of Mali. The spatial analysis aims to correctly identify the land properties and geographic setting. First, we define the image spatial subset using the grdcut module with target study area of Mali by argument ‘-R-13/5/9.5/25.5’. The map was visualised in Equidistant conic projection, a conformal equal-area projection, defined in module pscoast by the argument -JD-4/18/9.5/25.5/6.5i, where −4° and 18° stand for longitude (W) and latitude (N) of the projection centre, the 9.5°/25.5° are the two standard parallels depicting the southern and northern borders of the map of Mali in Figure 1, and the 6.5i defines the physical dimensions of the map as 6.5 inches.

Also, we consider plotting the isolines using the grdcontour module visualised with interval 250 m. The psclip module was applied for optional translucency of the background image of the neighbouring countries. Generating texts and annotations was done using pstext module. Once the geographic image was plotted, the psbasemap module was adjusted to visualise cartographic grid. We used the GEBCO/SRTM grid raster file for main image, and ETOPO1 for plotting the generalised contours. The textual annotations and country borders were adopted from the Digital Chart of the World (DCW) data. The final mapping procedure is given in Listings presented in GitHub.

The image processing was performed using R programming language (R Core Team, 2022) and selected libraries for graphical visualisation (Murrell, 2005). The key to successful image processing aimed at vegetation analysis and land cover types identification is correctly establishing band (or channels) correspondence. For images with present flooded areas and dominant savannah vegetation that have diverse composition and plant structure due to the effects of the climatic gradient, correspondence extraction is difficult since repeating structures creates mosaics of vegetation associations many of which are repetitive. Woody and herbaceous vegetation often intersperses by partitioning the landscapes into smaller patches.

To accurately classify savannah, grasslands and agricultural lands, we applied the automated method of k-means clustering and used well-established vegetation indices that indicate the presence of green healthy vegetation on the scenes. To this end, we computed several vegetation indices – NDVI, SAVI and EVI2 – for the region of Inner Niger Delta with methods and formulae described below for each index. The histograms showing frequency of values in pixels are computed for each vegetation index to indicate data distribution.

The Normalised Difference Vegetation Index (NDVI) was computed using a general formula of band combination for NDVI = (NIR − Red) / (NIR + Red) (Tarpley et al., 1984). In R syntax, the NDVI was defined using the code presented in GitHub. The SpatRaster was created from Bands 4 and 5 as a spatially referenced surface in RStudio (Figure 6). A raster object was made from the combination of Bands 4 and 5 and divided into the three-dimensional cells (rows, columns, and layers). The calculation and mapping of NDVI were followed by plotting the histograms that show the data distribution for each of the NDVI diapason of values. Here, each pixel from Band 5 (NIR) and Band 4 (Red) quantifies the vegetation biomass and plant vigour with a general range from −1 (min) to +1 (max). We used the R terra library and applied the code (a case for 2015, repeated likewise for images on 2013, 2018, 2020, 2021 and 2022).

Figure 6.

The calculation of the Soil Adjusted Vegetation Index (SAVI) (Huete, 1988) corrects the NDVI for the effects of soil brightness in desert regions, bare soil and other places where vegetative cover is low. The general formula used for calculation SAVI is as follows SAVI = ((Band 5 − Band 4) / (Band 5 + Band 4 + 0.5)) * (1.5) (Huete, 1988; Huete et al., 1994). Note that the numeration of these bands corresponds to 4 and 3 because the numeration in R syntax of RasterObject starts from 0 (zero). Therefore, in formula below, the NIR and Red bands are assigned numbers of 4 and 3 in the lapp function, as shown below. Using R syntax where NIR = Band 5 and red channel corresponds to Band 4, the computing of SAVI is performed using programming script.

The process of computing SAVI vegetation index is performed in R library terra using the code R (here: example for 2021, repeated likewise for all other images). The differences between soil colour and brightness around the regions near the Saharan desert, such as Mali, and affect the NDVI values. As a result, the lower boundary condition of identified plants may affect and bias land cover mapping in the Sahel region. To this end, we visualised the set of SAVI indices for the sequence of years from 2013 to 2022 and compared them to the previously computed NDVI visually and numerically.

Computing the Enhanced Vegetation Index (EVI) was based on using the equation of Landsat 4 and 5 bands: EVI2 = 2.4 * (Band 5 − Band 4) / (Band 5 + Band 4 + 1) (Huete et al., 2002). In R language, this formula was defined and converted using its syntax f or calculation and visualisation of the Enhanced Vegetation Index (EVI) as presented in Listing in the GitHub of the author.

The time series of EVI are used to recognise the agriculture lands and planting areas (e.g., detecting rice planting and harvest for millet, maize and sorghum) in contrast to the savannah, flooded areas and bare lands. Besides, variations of EVI over years obtained from the satellite images are useful for assessment of land degradation and desertification in semi-arid and arid regions of Sahel and southern Sahara of Mali.

The matrix of image cells is defined as a sequence of cells containing pixels that represent land cover types with 30-m resolution for Landsat 8–9 images. In k-means clustering technique, the consecutive massive of the pixels on the raster image with similar values is changed by the move to another land cover types which interrupts the mosaic of landscapes (e.g., water-land borders), and followed by the next landscape patch. The landscape mosaic is identified through iterative assignment of pixels to target land cover classes on the whole image in a repetitive way by the k-means algorithm of clustering. Thus, land cover types on the image are defined using search for the minimal distance between each pixel on the scene and the centroid of the k-means clusters, which are defined automatically by R. The assignment of the pixel to a given land cover class is based on the algorithm which evaluates the sequence of the cells in the image matrix based on the principle of the closest distance between the pixel and the cluster’s centre.

The k-means method was selected since it enables to allocate pixels into classes automatically, and thus to present a reliable correspondence between the groups of land cover classes and assigned pixels. Each iteration is recalculating the mean values of the class and re-assigning pixels to the new means until optimal values are reached using the embedded algorithms by R. For this purpose, the contrast of values in each land cover classes against the neighbouring classes has been tested using empirical trials of k-means algorithm which evaluates the distance of pixels to the centroids based on minimal distance principle. Finding correspondence between the value of pixels and land cover categories is a key operation in remote sensing data processing. Given several pixels identical up to spectral brightness, the task is to identify a set of similar pixels and group them into relevant classes well distinct from the others. This task is readily generalised by R library ‘RSToolbox’ using embedded algorithms of image processing and analysis of pixels on the raster scene. The k-means clustering technique is based on grouping the pixels on the raster image into clusters according to their spectral reflectance (Ose et al., 2016; Richards, 2022). The code of clustering algorithm is summarised in Listing in the GitHub.

Thus, the classification task is formulated as a pixel-based unsupervised clustering using k-means (Sreevalsan-Nair, 2020). The technical approach is performed by libraries terra’, ‘raster’ and ‘RSToolbox’ to efficiently discriminate land cover features and their changes over time. The target centroid of pixels was selected automatically with allocated neighbouring pixels assigned to groups based on the similarity of values. In such a way, the inter-annual variability of vegetation associations was detected in the inundated areas of the Inner Niger Delta of Mali. The code used for clustering is presented in the Appendix for a case of 2013 and applied with minor modifications for all the other images on years 2015, 2018, 2020, 2021 and 2022.

In the map of land cover types, we depicted 20 major dominating vegetation associations arranged from the most affected flooded areas to the relatively dry areas and untouched bare soils including the following 28 land cover types adopted from the existing studies and interpreted on the Landsat imagery Marie (2000): 1) Woodland; 2) Floodplain with aquatic grassland; 3) Savannah; 4) Thicket; 5) Wild rice; 6) Grassland; 7) Irrigated rice fields; 8) Savannah in shallow floodplain; 9) Bare soil (loams and clays); 10) Bare soil (sands); 11) Vegetation mosaics, main channel banks; 12) Lacustrine woodland; 13) Shallow lacustrine woodland; 14) Shrubby savannah in uplands; 15) Anthropic woody savannah; 16) Thicket in uplands; 17) Woody savannah in uplands; 18) Shrubby savannah in uplands; 19) Island palm grove; 20) Margins palm grove.

The correlation analysis was performed using plots of the non-parametric Kendall correlation method (Kendall, 1938) which is similar to the Pearson correlation method (Pearson, 1895). Both methods measure the association between the two measured quantities, that is, land cover classes and the probability of the correct association of the pixels assigned to these classes. In Python, these algorithms were implemented by the code presented in Listing in the GitHub (here, we tested both approaches and selected Kendall method). The correlation between two computed land cover classes is high when observations have a similar values, and low when observations differ substantially.

The procedure of categorising image into several classification clusters in image processing includes data evaluation accuracy, i.e., the comparison of the results with reference data. Since the fieldwork was not envisaged and the reference data limited, the assessment of the results was performed using correlation matrix for evaluating the pixels on the assignment to the correctly classified classes using coefficients for different land cover variables. To this end, the correlation plot was done for each set of 20 classes as an iterative loop for the images.

Complex cases in classification process included similar values of the land cover classes for both plant and water areas which consist of connected lakes, channels, and marshes of the delta. The classification control was done by double correlation of the classes and correctness of the pixels assigned to the different but similar classes according to spectral reflectance values using Kendall correlation matrix. To this end, the probability cases for land cover classes were computed in Python. The Kendall rank coefficient is used as a statistic matrix testing a hypothesis if land cover variables are statistically dependent when comparing pixels associated to these land cover classes. It includes the assessed land cover classes identified on the images for the period of 2013 to 2022 and enables the evaluation of pixels correctly classified in the target classes. The Pearson correlation approach estimates a linear correlation by assigning a value between −1 and 1 which indicates the degree and direction of the correlation between two land cover classes.

The results of the NDVI computation show changes in vegetation from 2013 to 2022 (Figure 7). The comparison of the computed values is summarised in Table 2 based on the SpatRaster from Bands 4 and 5 of the Landsat 8–9 satellite images and shown in histograms of Figure 8. The lowest NDVI values in the study area were obtained from the Landsat-9 image in 2022/11/11 as −0.3072812, while the highest data are 0.5879165 received from the Landsat-8 image in 2013/11/10.

Figure 7.

Figure 8 shows the frequency of the NDVI values by pixels on the image. Although theoretically, the NDVI ranges from −1.0 to +1.0, in this particular area of Mali, the actual range varies from −0.3 to 0.5 for the month of November for all years.

Figure 8.

Here, the areas of bare soil land (loams and clays) and sand have very low NDVI values (less than 0.1), while vegetation types have positive values. Woody canopy, island and marginal palm groves have values over 0.45 to 0.50; grassland − 0.2 to 0.3, rice fields from 0.20 to 0.30, floodplain wetlands with aquatic grassland − 0.40 to 0.45; savannah − 0.30 to 0.35; thickets in uplands − 0.35 to 0.40; vegetation mosaics including shrub and lacustrine woodland − 0.1 to 0.2.

Changes in vegetation associations showing mosaics of areas covered by savannah against floodplain with aquatic plants and other areas were estimated by using a time series of Landsat 8–9 satellite images acquired on different dates between 11.2013 and 11.2022 (Figure 9). The highest SAVI values were 0.4281315 on 2013/11/10, while in 2021/11/16, the highest SAVI was only 0.2790245, that is, the lower level of vegetation health (Table 2). The lowest SAVI values among all these years were obtained from image on 2013/11/10 which shows the most significant range of variation in values (Figure 10), which differs from the NDVI values in relevant years (see Table 2) that shows the lowest NDVI values on 2022/11/11 (NDVI=−0.3072812) but corresponds to the highest values obtained for 2013/11/10 (NDVI=0.5879165).

Figure 9.

Figure 10.

Figure 10 shows the frequency of the SAVI values by pixels on the image representing their associations with land cover types. The higher the SAVI values are, the denser and more healthy is the vegetation coverage. The actual range of values of the Inner Niger Delta, Mali, for the taken images varies from −0.15 to +0.35, although the general theoretical range of values may change between −1.0 and +1.0, similar to the NDVI. Since SAVI is adjusted for soil brightness for areas with low vegetative coverage, it is particularly suitable for the African Sahel region with desert and sandy areas in arid and semi-arid climate. Woodland and palm groves have values over 0.30 to 0.35; grassland − 0.2 to 0.3 which is similar to NDVI, vegetation mosaics including shrub and lacustrine woodland − 0.10 to 0.15, areas covered by irrigated and wild rice from 0.15 to 0.20, thickets − 0.20 to 0.25 and savannah − 0.25 to 0.30, the lush wetland vegetation in flooded areas with aquatic plants ranges between 0.25 and 0.35, i.e., have the highest values of SAVI.

The radiometric and biophysical advantages of EVI (Figure 11) for mapping vegetation canopy in Malian environment is that compared to NDVI which better saturates in high biomass regions such as tropical forests, EVI remains sensitive to canopy variations that is suitable to savannah of Sahel arid and semi-arid regions. Specifically, EVI shows a correspondence between the canopy reflectances of plants as detected on the pixels and its values obtained from the Landsat 8–9 sensor specifically for semi-arid grasslands, thicket and shrubby savannah in uplands, savannah in shallow floodplains of the Inner Niger Delta, and spots of woody savannah.

Figure 11.

The values of EVI (Table 2) were compared to the other two different well-known vegetation indices (NDVI and EVI) and were calculated and used from each satellite image as a time-series format of six scenes. EVI values calculated from Landsat 8–9 satellite images are higher than earlier satellites of the Landsat mission because of their lower resolution, e.g., compared to Landsat-MSS or TM, the extent of minimal and maximal values has a greater span of values. Clearly visible two peaks on the histograms (Figure 12) correspond to the flooded inundated water-covered areas with negative values (EVI from 0 to −0.35 in a general extent). The lower data are typical for floodplains and aquatic plants which contrast against the bright vegetation areas that correspond to the grasslands, savannah and thickets (EVI values from 0.30 to 0.68). Finally, the middle values that are typical for the agricultural areas (rice) have mediocre values (EVI values from 0.20 to 0.30).

Figure 12.

The values of EVI range from −1 to +1. However, for the specific area of the Inner Niger Delta, Mali, on the Landsat 8–9 images taken on November, it varies from −0.25 to +0.40. For healthy wetland vegetation in, it varies between 0.05 and 0.40 which well represents the second peak in the histograms, while the first peak (−0.20 to +0.05) corresponds to flooded areas of delta and wetlands. Since EVI is corrected for aerosol scattering present in the atmosphere and is adjusted for soil and canopy background, it well detects canopy background and identifies noise. Moreover, EVI is sensitive to land cover types with dense vegetation.

The examination of data (Table 2) shows that EVI values were extremal in 2013/11/10 when they ranged between −0.3576084 and +0.6849960, which corresponds to the consequences of the large drought in 2012 in the Sahel region which was hit by a major climate extremes (decrease in precipitation and rise in temperatures) which weakened vulnerable vegetation communities in the consequent years, as reflected in the values of the vegetation indices in 2013. In the consequent years, these effects gradually diminished along with the stabilisation in the ecosystem.

The image classifications aimed at defining land cover types within the Inner Niger Delta of each Landsat 8–9 scene were interpreted for each of the six evaluated years within the 2013–2022 time gap. Land cover types are classified based on the k-means clustering using Landsat 8–9 images of Inner Niger Delta for November and visualised in Figures 13 and 14.

Figure 13.

Figure 14.

The notation for land cover classes are as follows: 1) woodland, 2) floodplain with aquatic grassland, 3) savannah, 4) thicket, 5) wild rice, 6) grassland, 7) irrigated rice fields, 8) savannah in shallow floodplain, 9) bare soil (loams and clays), 10) bare soil (sands), 11) vegetation mosaics, main channel banks, 12) lacustrine woodland, 13) shallow lacustrine woodland, 14) shrubby savannah in uplands, 15) anthropic woody savannah, 16) thicket in uplands, 17) woody savannah in uplands, 18) shrubby savannah in uplands, 19) island palm grove, 20) margins palm grove.

Coloured land cover classes show areas of agricultural and natural zones of vegetation as well regions with water and bare lands. The area of delta floods and vast flat area for the one to three year gaps are visualised in the processed images to show the low-lying floodplain with expanses of wetland grasses and reeds classified against other land cover types. The land cover classes are identified based on the frequency of the values of spectral reflectance of pixels on the images encountered once on each Landsat scene with regard to the identified land cover types. The difference in the northern and southern parts of the delta consists in varying land cover types. For instance, while northern part has sand ridges that are formed during the flood period, southern portion of the delta is a vast alluvial plain.

For each data set, the image scene was divided automatically into 20 classes with similar properties of spectral reflectance. The random data partition is repeated for 5 cycles with 100 samples in each go, and the results of k-means clustering over the 20 classes are reported and presented in the images. The classified land cover classes are sampled and labelled, as shown in the resulting images in Figure 13 for years 2013, 2015 and 2018 and Figure 14 for years 2020, 2021 and 2022, respectively. At each iteration of the k-means clustering, the RSToolbox query one instance label pair of pixels and compares the values of target pixel with that of the centroid pixel representing the given group of the land cover classes from the input Landsat Bands 4, 3 and 2. After all instance-labelled pairs of pixels are queries (target-centroid), a classification model is trained and computed for labelled pixels as data from bands. The obtained statistics is summarised in the Appendices. These values were applied to estimate the effects of hydrological-climate fluctuations on the observed expansion of vegetation associations at target areas covered by the images.

For each time window, the scenes were classified automatically using R libraries terra, raster and RSToolbox with machine-based allocations of pixels to the land cover classes defined according to the classification. The data are evaluated per land cover classes using Kendall correlation matrix to examine changes in classification from 2013 to 2022 (Figure 15). The summary of land cover changes is given in a table with obtained statistics on pixels.

Figure 15.

Using correlation plot by Python libraries (Figure 15), the associations of 20 major land cover types were compared and identified for each scene as a comparison between 2013 and 2022 to find out the correlations between values in previous and later years in spectral reflectance of the pixels on the images and the distribution of the vegetation associations in Inner Niger Delta with the most important land cover classes: agricultural areas – irrigated rice fields (class 7), natural vegetation areas – grassland (class 6), savannah of various types – shrubby, anthropic woody and uplands savannah in shallow floodplain (class 8), thicket in uplands, bare soil – loams and clays (class 9), bare soil – sands (class 10) and specific plants typical for floodplain of the Inner Niger Delta with characteristic aquatic plants (shallow lacustrine woodland) (class 13), island palm groves (class 19) and margins palm grove in the inundated areas and channel banks (class 20) (Figure 15).