Spatio-temporal monitoring of vegetation phenology in the dry sub-humid region of Nigeria using time series of AVHRR NDVI and TAMSAT datasets

2.1 Study area

This study focused on Kwara state in the middle belt region of Nigeria where significant environmental changes are visible. The state geographically shares boundary with Niger state to the North, Oyo, Osun and Ekiti to the south and Kogi state to the east and lies between latitudes 8°–10°N and longitude 3°–6°E. Kwara state is bounded by River Niger from the east to the North and it is often referred to as savannah ecotone because it forms the transition zone between forest and savannah region of Nigeria [28, 29]. The state covers an area of about 36,825 square kilometer with population of 1548412 (in 1991), 2365353 (in 2006) [30] and population density ranging from 20-120 people per km² as shown in Fig. 1.

Figure 1

Population density of the study area for 2010 and 2014 (Data is 100 × 100 m) Source: www.worldpopulation.org.UK [31]

The annual rainfall of Kwara state is between 1000 mm to 1500 mm with double rainfall peaks in June and September [28, 32]. The region is characterized by two season’s namely dry and wet season, the raining season begins in March and ends in October while the dry season begins in November and ends in February. The temperature is between 25°C–30°C in wet season and 33°CO-34°C in dry season. The soil is mainly classified as ferruginous tropical soils and support crop production [33]. The vegetation is mainly forest and savanna and it constitutes about 47.78% and 35.04% respectively [29]. The common tree species found in this region includes Vitellaria paradoxa (shea butter), Acacia spp, Parkia spp, Afzelia Africana, Terminalia spp. Kwara state enjoys favorable climatic condition which supports the growth of tall grasses such as Andropogon gayanus, Ageratum conyzoides (goat weed grass), Pennisetum purpureum (elephant grass) and Heteropogon contortus (Spear grass) [34]. Due to the increasing rate of anthropogenic activities in this region, assessment of change in vegetation phenology is very important. Figure 2 shows the mean monthly NDVI pattern of the study area from 1983-2011 as derived from the

Figure 2

Location of the study area showing selected areas across different land cover types in Kwara state, Nigeria (a) open broadleaved deciduous forest/woodland, (b) Cropland, (c) shrubland and (d) grassland. Source. www.GADM.org [35]

2.2 Datasets

Bi-weekly NDVI 3g datasets (15 days composite) from Global Inventory Modeling and Mapping Studies with spatial resolution of 8 km was used in this study (https://ecocast.arc.nasa.gov/data/pub/gimms) for the annual phenological assessment from 1983-2011. The decomposition and reconstruction of the NDVI datasets had already been done in order to correct solar zenith angle, orbital drifts or artifacts such as volcanic aerosols which might cause systematic trends in the datasets using a technique known as Empirical Mode Decomposition (EMD), this makes it different from the older version of NDVI product [27, 36, 37, 38]. The longer span of the new NDVI3g datasets (more than 30 years) makes it suitable for monitoring seasonal vegetation trends in the dry sub-humid region of Nigeria where technical resources and adequate information about change in vegetation phenology is limited. The second data used in this study is the Tropical Applications of Meteorology SATellite (TAMSAT) which is basically for Africa (https://met.reading.ac.uk/tamsat/data). The monthly rainfall datasets had been calibrated using historic rain gauge records, and it has a spatial resolution of 4 km [39]. GlobCover 2000 and 2009 datasets from the European Space Agency (http://www.esa-landcover-cci.org) was used in this study in order to have in-depth understanding of vegetation change dynamics across different land cover types. The global land cover map for 2009 has a spatial resolution of 300 m and it is derived from MERIS sensor on ENVISAT while the 2000 global land cover has a spatial resolution of 1km and is derived from SPOT-Vegetation [40].

2.3 Methods of image pre-processing

All the datasets used in this study were reprojected to the Universal Transverse Mercator (UTM) and subsets of the satellite imagery were created using the administrative boundary of Kwara State. Prior to the Seasonal Trend Analysis (STA), the bi-weekly NDVI composite data were aggregated to a monthly composite using Maximum Value Composite (MVC) method which further reduces noise in the time series [41]. The rainfall data was resampled to 8 km spatial resolution in order to match with the spatial resolution of using bilinear interpolation resampling method. The global land cover map of 2009 was resampled to 1km resolution so as to merge with the same spatial resolution of the 2000 global land cover. The land cover class for both maps was aggregated to six major classes, an overlay analysis was applied to identify the change matrices for the time period of 2000 to 2009. Assessment (STA) was performed using three approaches as described below.

2.3.1 Analysis of seasonal trend

Change in NDVI trend across the dry sub-humid region of Nigeria was analyzed using STA approach as described by Eastman et al., 2009. This method was selected based on its ability to detect seasonal signals in any time series datasets whilst rejecting noise and inter-annual variability [42, 43]. STA uses two stages of time series analysis to detect NDVI trend. This stage is then preceded by visualization of the produced image (amplitude and phase image) at the final stage [42]. STA performs harmonic regression of annual (yearly) images in the time series as described in equation 1. In the first stage of analysis, each annual imagery in the time series (1983-2011) are subjected to harmonic regression in order to approximate the seasonal curve thereby producing five greening parameters, amplitude 0 (annual mean image), amplitude 1 and phase 1 (annual cycle), amplitude 2 and phase 2 (the semi-annual cycle) [42, 44, 45]. The overall greenness (mean annual NDVI for each year) is described as Amplitude 0, the magnitude of peak of annual greenness as amplitude 1, while the timing of annual peak of greenness is referred to as phase 1 [44]. An increase or decrease in phase 1 indicates shift in the timing of greening event to an earlier or later time of the year respectively. The generalized equation for the harmonic regression is presented as follows:

where t is referred to as time, y is the series value, T is the length of the time series, n is a harmonic integer multiplier, e is error term, α₀ is the mean of the time series while a_n and b_n are regression parameters. The harmonic regression curve can also be expressed as:

after the term of equation (1) had been rearranged and error term had been omitted or ignored. Hence, α_n and φ_n are referred to as amplitude and phases [42].

In the second stage of the analysis, annual trends in the greenness parameters are calculated [42]. This stage involves the calculation of annual trends in the greenness parameters (amplitude 0, amplitude 1 and phase 1) using Theil Sen trend detection estimator. Theil-Sen (TS) trend estimator is a robust non-parametric technique which calculates all pairwise combination of images in time, it takes the median of all slopes to generate trend images for each of the greenness parameters [46, 47, 48]. The TS trend estimator is insensitive to outliers and it has a breakdown bound of 29% [42, 43, 49, 50]. This simply implies that the values of the TS slope tolerate outliers up to 29% of the observations used in the time series without degradation of its accuracy [44, 51]. STA was also applied to the rainfall (TAMSAT) data for better comparison with NDVI.

2.3.2 Significance test of STA using Contextual Mann-Kendall (CMK)

Trends in vegetation dynamics for each pixels (NDVI time series) was statistically examined in order to identify areas undergoing similar trends [44]. In order to evaluate the significance of trends during the Theil-Sen slope calculation (see section 2.3.1), a non-parametric statistic test known as Contextual Mann-Kendall (CMK) was applied [52]. CMK test is resistant to outliers and uses similar approach like that of Mann-Kendall [52]. To perform CMK test, three stepwise approach are involved. The first step account for the elimination of serial correlation in the NDVI and rainfall time series which might have influence on the trend test using a method known as prewhitening [44, 53]. Prewhitening removes serial correlation from the residuals of both NDVI and rainfall while maintaining the trends in the dataset [55, 56]. The trend in the prewhitened series is similar to the original data while having no serial correlation [53]. In the 2^nd stage, trend consisting of 3 × 3 neighborhood around each pixel are evaluated, while in the 3^rdstage, the calculated 3 × 3 neighborhood result which is produced in the 2^nd stage was used for calculating spatial autocorrelation [44]. The result of the CMK trend test are in the form of image i.e statistical p-value and z score each describing statistical significant trend (positive and negative trend) in the greenness parameters (amplitude 0, amplitude 1 and phase 1). Generally, a positive z score indicates upwards trend where as a negative z score represents a downwards trend [52].

2.3.3 Interpretation of phenological curve

An interpretational approach as described by Eastman et al, 2009 was used to describe variation in the start and end of the series based on the slope and the intercept of the five harmonic parameters (amplitude 0, amplitude 1, phase 1, amplitude 2 and phase 2) across different land cove types for selected locations. The selected locations were sampled on a single pixel base (64 km²), seasonal curve for the beginning and the end of the series were produced by using temporal profile [Eastman, 2012]. The differences between the two curves indicate resultant changes which can be used to describe trends in the observed phenological pattern of vegetation which are useful for environmental studies and ecosystem monitoring [27, 44]. The term greenup and green-down are used inter-changeably to represent the starts of season (SOS) and end of season (EOS) respectively.

3.1 Spatial analysis of vegetation trend pattern

Time series of (NDVI) from 1983-2011 was used to provide information on vegetation dynamics in the dry sub-humid region of Nigeria. The results of the spatial pattern of statistical significant trends in overall greenness (amplitude 0), peak of annual greenness (amplitude 1) and timing of annual peak greenness (phasel) over the study area as shown in Figures 3a-c. Both positive and negative trends were observed from the results of the greenness parameters. Positive pixels in the spatial image (significant image) indicates an increasing trends while negative pixels indicate decreasing trendsin vegetation dynamics over the observed period. Significant positive trends in overall greenness were observed in Fig. 3a while negative and nonsignificant trend pattern areobserved for fewer pixels (p< 0.05).

Figure 3

Spatial map of Kwara State showing the result of the Contextual Mann-Kendall statistical significance (p-value) calculated from GIMMS NDVI3g data for 1983–2011: (a) overall greenness, (b) peak of annual greenness and (c) timing of annual peak greenness

Figure 3b shows significant pattern in peak of annual greenness.Positive trends are observed particularly in the northern part while negative trend were observed in the western and southeastern part. Pixels with nonsignificant trend in the peak of annual greenness are distributed across the study area.Figure 3c shows the significant trend pattern in timing of annual peak of greenness, more than 90% of the pixels used in this analysis showed a negative trend, this indicates a shift in the pattern of timing of annual peak of greenness to a later time, while the fewer pixels which shows a positive trend indicates earlier green-up.

Although large clusters of positive NDVI trends in overall greenness (amplitude 0) were observed in the study area during the examined period (1983–2011) as shown in Fig. 4, areas with negative NDVI trends for example the black circle (Fig. 4a) can be attributed to land degradation due to anthropogenic activities.

Figure 4

Spatial pattern of linear trend of mean annual NDVI for the period of 1983-2011 (NDVI/month): (a) Annual mean NDVI trends (amplitudeO), (b) p-value

The results of the mean rainfall time series analysis (amplitude 0) show significant positive trends with variations in the slope direction (Fig. 5). This variations however has impact on vegetation phenology.

Figure 5

Spatial pattern of linear trends of mean annual rainfall calculated from monthly TAMSAT for the period of 1983-2011 (mm/month): (a) Annual mean rainfall trends (amplitudeO), (b) p-value

The result of the spatial regression between rainfall data and NDVI over the study region is presented in Figure 6. The result shows a high relationship (coefficient or determination, R²) between NDVI and rainfall which indicates that rainfall is the main driving force for vegetation growth in sub-Sahara and tropical regions like Nigeria [5, 56, 57]. Other studies also showed that incoming shortwave radiation and air temperature are not constraints of NDVI trends compared to rainfall distribution [54].

Figure 6

Spatial regression of NDVI and rainfall over the study region for the period 1983-2011, the values indicate pixel-wise R² of NDVI vs rainfall

Figure 7 shows the result of timing of annual rainfall (phase 1). Shift to a later time was observe from the result. This implies that the early months of the year received less rainfall than the later months of the year. An increasing or decreasing phase angle depicts a shift to earlier or later time of the year.

Figure 7

Spatial pattern of Timing of annual rainfall (phase1) calculated from monthly TAMSAT for the period 1983-2011 (mm/month)

3.2 Visual interpretation of phenological curves

Four locations were selected to visualize the monthly mean NDVI profile curves for the beginning and end of the time series (1983-2011) across different land cover types in the study region (Fig. 8). The seasonal curves for the examined years were compared for (a) open broad leave deciduous forest/woodland (b) cropland (c) mosaic grassland and (d) closed to open shrubland as shown in figure 8. Seasonal variation in vegetation cover was observed from the results of the NDVI curves (phenology) for the selected locations. The differences between the two curves (1983 and 2011) however indicates changes in trend pattern.

An increase in overall greenness (NDVI) for all months (except August, as there is a decline in 2011 in this month) was observed for the land cover type (Open broadleaved deciduous forest/woodland) in location a with an increase in vegetation productivity in 2011 while decreasing trend was observed in the other greenness parameter (peak of annual greenness and timing of annual peak of greenness). The maximum difference between the two curves occurred in the January months, October, November and December while minimum difference occurs in March, July, August and September (Fig. 8a). This location is within the Borgu game reserve where agricultural activities is not allowed, although some locals still do some illegal cultivation in the buffer zone of the game reserve.

Location b shows significant increase in overall greenness and peak of annual greenness, while a significant decrease was observed in the timing of annual peak of greenness for this region. The land cover type in this location is associated with cropland where intensive agriculture is being practiced almost throughout the year. The minimum difference between the two curves occurred in June and September while the maximum difference occurred in November and December (Fig. 8b).

Location c & d exhibits significant increase in overall greenness and significant decrease in other greenness parameters for the study period. Two peak of greenness was observed from the seasonal curve of both locations (Fig. 8c and d). The land cover types in these location are shrub and grassland respectively, an increase in NDVI productivity was observed during the growing season, it reaches its peak in the month of June and gradually decline in August i.e a short pause in the raining season lasting for two to three weeks which is known as August break. Minimum difference between the two curves (1983 and 2011) occurred at the start of the growing season (start of season) in February while maximum difference occurred in July, August, October, November and December months for location c (Fig. 8c). For location d however, minimum difference occurs in March, April, July and September while Maximum difference occurred in January, June, and October (Fig. 8d).

Figure 8

Seasonal curves of NDVI for four selected locations in kwara state (1983 in green and 2011 in red). The graph shows the year from January through December on the x-axis, and NDVI × 10000 on the y-axis.

Vegetation greenness in the dry sub-humid region of Nigeria generally increased over the period of study (1983–2011), this can be seen from the result of amplitude 0 (significance image of overall NDVI greenness) in Figure 3a. Variation in the phenological trend pattern was confirmed by the result of the NDVI seasonal curve across different land cover types in the selected locations.

Change in land use land cover was revealed by the result of the GlobCover between 2000 and 2009.

Table 1 summarizes the distribution and size of land cover classes in the year 2000 and 2009. It indicates a decline in sizes of forest and grassland while an increment in the size of cropland and shrubland was observed. This clearly indicates land use transformation is dominantly occurring with the influence of anthropogenic activities.

An assessment on the land use dynamics using change matrices for the year 2000–2009 (Fig. 9) also confirmed significant shifts in vegetation dynamics. Accordingly, there is a substantial transition from the natural ecosystem (Forest/Shrubland) to human dominated activities, dominantly cropland.

Figure 9

Change dynamics for the period 2000–2009 (FO-Forest, CL-cropland, GL-Grassland, SL-shrubland)