ASSESSMENT OF VEGETATION DYNAMICS AND FOREST LOSS USING GOOGLE EARTH ENGINE AND MULTI-TEMPORAL SENTINEL-2 IMAGERY

This study evaluated regional vegetation dynamics and changes between 2015 and 2020 using Google earth engine (GEE) platform and normalized difference vegetation index (NDVI) derived from the multi-petabyte catalogue of sentinel-2 imageries. Using the computational capability of GEE, yearly mean NDVI from 2015 to 2020 were computed using level C-1 product. Subsequently, each of the NDVI images was classified into four land cover classes; water bodies, non-vegetated, grassland/cropland/shrubs, and forest using NDVI threshold values of < 0.01, 0.01-0.20, 0.20-0.30 and > 0.30, respectively. The classified maps allowed for the assessment of yearly variation in vegetation and changes between 2015 and 2020. Result showed that non-vegetated area increased from 18.53% in 2015 to 42.56% in 2020 (~ 25.00% gain), the forest area reduced to 6.78% in 2020 compared to 23.76% measured in 2015 (~ 17.00% loss in forest); whereas water bodies and grassland/cropland/shrubs remained relatively constant (0.21 and ~ 50.00%, respectively) across the years studied. Presently, the forest land was estimated to be about 2, 371.131 km 2 (~ 6.70%) of the total land mass, grassland/cropland/shrubs occupied 17, 770.79 km 2 (~ 50.07%), non-vegetated area was slightly less than half with 15, 274.85 km 2 (~ 43.04%) and water bodies occupied 75.68 km 2 (~ 0.21%).


INTRODUCTION
Remote sensing has contributed immensely to monitoring climate change on a global scale in near real time through analysis of huge volume of data acquired with different earth observation systems (EOS). Continuous acquisition of data over time have resulted in huge amount of data beyond the computational capability of the traditional desktop image processing packages to meet timely analysis of earth processes (Rani et al., 2018). This has led to the development of Google earth engine (GEE) and other online computation mapping platform. GEE is a cloud computing platform for repository, processing, and analyses of satellite imagery and other geospatial dataset. It provides a single platform to access various database of a multipetabyte catalogue of satellite imageries and the computational power required to computes the images (Gardner, 2010;Gorelick et al., 2017). GEE platform hosts diverse dataset from different EOSs such as; landsat, moderate resolution imaging spectroradiometer (MODIS), copernicus sentinel, and shuttle radar topography mission (SRTM).
GEE cloud resource runs on application programming interface (API) of Javascript and python; also, it has gained popularity among the geospatial community as a powerful cloud resource for geospatial analysis (Google, 2021). For example, Jesus et al. (2020) reportedly used GEE for image acquisition and derivation of several indices such as normalized difference vegetation index (NDVI), green normalized difference vegetation index (GNDVI) and visible atmospherically resistant index (VARI). In another study by Khan and Gilani (2021), GEE has assisted in monitoring drought indices such as vegetation condition index (VCI) for detecting changes in vegetation condition, temperature condition index (TCI) for examining temperature changes, soil moisture condition index (SMCI) for accessing degree of dryness and wetness of soil, alongside the precipitation condition index (PCI) computed from satellite derived data.
Vegetation is a major component of the ecosystem; it regulates various biogeochemical processes such as water, carbon cycling, and nitrogen (Obalum et al., 2012;Xu and Trugman, 2021). Vegetation converts solar energy into biomass and forms the base of all food chains and habitat for many organisms (CNVC, 2013). There have been several studies on the use of NDVI for monitoring vegetation (Weier and Herring, 2000;Bhandari et al., 2012;Al-doski, 2013;Gandhi et al., 2015;Jing et al., 2016;Nouri et al., 2017;Schmid et al., 2017;Rani et al., 2018;Gessesse and Melesse, 2019). NDVI provides indicator for many physiological and ecological parameters including surface, vegetation, albedo and photosynthetic activities, land cover chaand drought assessment (Govaerts and Verhulst, 2010). At continental level, multi-temporal data sets of NDVI have been used to classify vegetation cover types for Africa (Goward et al., 1985) and South America (Townshend et al., 1987) using NOAA's global vegetation index (GVI) product at approximately 15.00-20.00 km spatial resolution. Also, Loveland et al. (1991) employed Advanced Very High Resolution Radiometer (AVHRR) vegetation index and unsupervised classification method to define homogeneous land cover regions for the conterminous United States.
The references provided above emphasize the importance of NDVI in vegetation monitoring, though at coarse resolution incapable of detail evaluation of subtle changes in vegetation amount. The present study employed GEE cloud platform to assess yearly vegetation dynamics from 2015 to 2020 using high-resolution sentinel-2 satellite imagery. Sentinel-2 is one of the fleets of dedicated European Union owned satellites designed to provide the wealth of data and imagery that are central to Europe's copernicus environmental programme. Sentiniel-2 has spatial resolution between 10.00 and 60.00 m and revisiting period of 10-day cycle under the same viewing angles with a swath view field of 290.00 km contrast to landsat-8, which has spatial resolution of 15.00 and 100.00 m and a swath width view field 185.00 km swath with a repeat cycle of 16 days (Missions, 2014;SUHET, 2015a;Teluguntla et al., 2018). GEE platform was utilized in this study because of its enhanced access to large satellite dataset forthwith and also allows processing of the dataset in its high performance open cloud platform in quick sequence (Thuận, 2019). This research focused on the dynamic changes of vegetation to non-vegetation area in Kwara State and investigated the pattern in which both climate and anthropology affected these changes.

Study Area
This study is a state-wide assessment of vegetation distribution across Kwara State in the North-Central geopolitical zone of Nigeria ( Figure 1). Geographically, Kwara lies within latitude 7.96° N and 10.15° E and longitude 2.72° E and 6.21° E. In terms of climate, the annual precipitation of the State ranges between 0.20 and 172.70 mm, occurring between Apr. and Nov., while annual average temperature ranges between 21 and 38°C (Idrees et al., 2021;Ozor et al., 2021). The State experiences harmattan period between Dec. and Jan., usually associated with cold and dry wind. The State has abundant fertile soil that favours agriculture and several tributaries that drained into Niger River to support irrigated farming. However, there has been report of increased rate of illegal logging for both lumbering and wood coal (charcoal) across the State (Akinyemi, 2018), in addition to urbanization and the impact of climate change which motivated this study.

Data Description
The datasets used in this study were sentinel-2, level C-1 imagery for the NDVI, while the precipitation and temperature data were obtained from the European centre for medium-range weather forecasts (5th generation of European re-analysis). In addition, the vector data representing the administrative boundary of the State was also obtained ( Table 1). The vector data representing the boundary of Kwara State was obtained from the global administrative unit layers (GAUL) spatial database of the FAO. The main purpose of GAUL was to develop, compile and distribute well-founded geographic information (GI) database on administrative units' levels for all the 247 countries (both UN and non UN members). The GI has contributed to the structuring and consistency of spatial dataset of countries administrative units right from the national level to the local level (Team, 2007).
GAUL makes use of three stages; data collection, evaluation and integration for distribution to authorized users (Fabio, 2016). The sentinel-2 level C-1 is one of the products of the sentinel-2 satellite. Sentinel-2 was launched by Europe's copernicus environmental programme for EOS, sponsored by European Union under the administration of European space agency (Navarro et al., 2017;Roteta et al., 2019;Idrees et al., 2021). The grid extent of the sentinel-2 level C-1 product is a cluster of 100.00 × 100.00 km 2 tiles in Universal Transverse Mercator/World Geodetic System 1984 (UTM/ WGS84) projection (SUHET, 2015b). Its product is top of atmosphere reflectance re-sampled into ground sampling distance of 10.00 m (band 2, 3, 4 and 8), 20.00 m (band 5, 6, 7, 8A, 11 and 12) and 60.00 m (band 1, 9 and 10) (SUHET, 2015b). ERA5-land monthly average is a monthly re-analysis of ECMWF climate that spanned over 35-year period and also the fifth generation of ECMWF atmospheric reanalysis of the global climate covering the saidperiod (Copernicus Climate Change Service, 2017). The ERA5-land monthly observes the earth atmospheric, land and oceanic climate variables such as land surface temperature, precipitation, runoff, evaporation, etc. The data has a spatial resolution of 30.00 km grid and temporal resolution of one month (Muñoz-Sabater, 2019). For this study, land surface temperature (LST) and precipitation were the only data extracted from ERA5-land monthly. Precipitation influences land surface components and hydrological processes, whereas LST emphasizes the relationship between land surface features and reduction in annual total precipitation (Jing et al., 2016).

Data Preparation
Processing image in GEE involves a number of steps, including data search, filtering, sorting, subset to study area, image processing and analysis. All these involve some level of coding in the GEE API. Figure  2 presents the data processing and analysis steps.

Data search
GEE is a multi-petabyte analysis-ready data catalogue (https://developers.google.com/earth-engine/datasets) which is built for effective performance computation service of geospatial datasets (Chastain et al., 2019;Hay-Chung et al., 2021). To analyze geospatial data in GEE, the first task is to access data repository website (https://code.earthengine.google.com/) in the browser by calling out the dataset using Javascript language function (ee.FeatureCollection and ee.ImageCollection). This process is achieved through internet accessible API coded in Javascript and python on a web-based interactive development environment that facilitate quick feature model and result visualization (Gorelick et al., 2017). Using these resources, the administrative unit of Kwara State, sentinel-2 and climate data (precipitation and land surface temperature) were retrieved.

Data filtering and sorting
Filtering is the process of retrieving all the required dataset for the study. From the GUAL 2015 firstlevel administrative units of the GEE repository, the shape file of boundary was filtered using the table schema (ADM1_NAME). Next, sentinel-2 level C-1 products from 2015 to 2020 were filtered from the GEE repository using both the acquisition date and the boundary file of the study area. For quality assurance and control, the resulting filtered data were further sorted using the image properties (Gorelick et al., 2017). For example, sentinel-2 data with cloud cover (cloudy pixel percentage) greater than 1.00% were excluded from the filtered images using in-built functions of the GEE. Thereafter, the mean of each image band for each year was generated and clipped to area of interest (Schmid et al., 2017).

Image processing Generation of NDVI image
Using the mean annual image bands, the mean NDVI image of the respective year were generated using near-infrared (NIR) and red bands using the expression in Equation 1 (Gessesse and Melesse, 2019). In this study, the red and NIR bands of sentinels -2A and -2B with wavelengths of 664.60 nm and 664.90 nm, and 832.80 nm and 832.90 nm for the Red and NIR bands, respectively were used (Bhandari et al., 2012). (1); where NIR and RED are the near infrared (band 8) and red (band 4) of the sentinel-2 imagery. The resulting NDVI images were re-projected to the universal transverse mercator (UTM) coordinate system, zone 31 north, world geodetic system 1984 (WGS84) datum (Brown et al., 2018;Idrees et al., 2021).  By standard, NDVI values range from -1.00 to +1.00 (Al-doski, 2013;Nouri et al., 2017;Schmid et al., 2017;Rani et al., 2018;Idrees et al., 2021) where surface features with high moisture content such as water, snow and cloud reflect more in the visible band than in the near-infrared band, resulting in negative NDVI value (-1.00 to 0.01). Characteristically, hard surfaces including natural and manmade like bare soil, rock, road, and building are usually identified with NDVI value of around zero (< 0.10), while vegetation exhibits stronger nearinfrared reflectance resulting in varying NDVI values from around 0.30 to + 1.00, depending on the health or degree of vegetation greenness (Aldoski, 2013;Gandhi et al., 2015). The NDVI values obtained across the six images ranged from -0.32 to -0.72. The values steadily dropped between 2015 and 2020 (Table 2), which signals decrease in either vegetation amount or greenness.

Land Cover Classification and Change Detection
To evaluate vegetation distribution for each year, the NDVI images were classified into four land cover classes (water bodies, non-vegetated area, grassland/ cropland/shrubs and forest) using NDVI thresholds (Table 3) determined based on literature sources (Bhandari et al., 2012;Al-doski, 2013;Gandhi et al., 2015). Changes in biomass (vegetation) between two consecutive years were assessed using NDVI differencing (DNDVI) (Al-doski, 2013).
NDVI difference (Equation 2) is a pixel-based operation computed by subtracting NDVI image of earlier year from the NDVI image of previous year. The change detection analysis is a direct approach to measuring changes between a pair of images that represent an initial state and final state of the environment in the imagery. The input images may be single-band images of any data type. The difference is calculated by subtracting the initial state image from the final state image (Equation 2), and the classes are defined by change thresholds. A positive change identifies pixels that became brighter (i.e., the final state brightness was greater than the initial state brightness), while a negative change identifies pixels that became dimmer (i.e., the final state brightness was less than the initial state brightness) (Al-doski, 2013).  where NDVIf is the image of the latter (final) year and NDVIi is the image of previous (initial) year. In addition to assessing biomass change between consecutive years, the amount of vegetation loss between 2015 and 2020 was estimated using postclassification change detection technique (Al-doski, 2013). Pixel-based change detection does not allow quantifying class-to-class change; so, a spatial context change detection technique using classified NDVI images was implemented. For each initial state class, the analysis identifies the classes into which those pixels changed to in the final state image, making it easy to identify not only where changes occurred but also the class into which the pixels changed to (Gandhi et al., 2015). Finally, the precipitation and temperature data extracted from ERA5-land was analyzed in order to understand the extent climatic factors impacted vegetation amount and distribution across the years.

Land Cover Classification
The land cover classification maps produced from the mean annual NDVI using threshold values allowed assessing vegetation amount, variation and distribution for each year from 2015 to 2020. Figure  3 presents the vegetation amount and distribution across the years under study. In the figure, it can be observed that the south eastern part of the state and some parts of the western side were generally more forested than any other part. This is further explained from the NDVI value (> 0.40), which indicates a measure of greenness that remained fairly consistent across the years considered, except 2018 and 2020; where significant reduction in vegetation amount and vigour was recorded. Also, considerable change from grassland to non-vegetated (NDVI < 0.20) is noticed across the northern fringe but intensified in the northeast. The northwest and the central area exhibited vegetation characteristics of grassland, cropland and shrubs, occasionally intruded with light forest. The result has shown a pattern that indicates alternating year of increase followed by decrease in vegetation amount (compare Figures 3a  and 3b, 3c and 3d, and 3e and 3f). Table 4 provides the quantitative evaluation of the maps in Figure 3. For example, the class percentage for water bodies varied between 0.21 and 0.23%, which shows that there is no significant change in the percentage coverage of water bodies. Similarly, across the years, the grassland/cropland/shrubs maintained averagely 55.50% area coverage (minimum of 50.45% in 2020 and maximum of 61.38% in 2017). In contrast, there was significant increase and decrease in the non-vegetated and forest land, respectively (Figure 4). The non-vegetated area increased from 18.53% in 2015 to 42.56% in 2020; conversely, forest area reduced to 6.78% in 2020 from 23.76% measured in 2015. Figure 4 shows the comparative class percentage plot.    Figure 5 is the change detection map where positive change represents pixels that became brighter as result of greater brightness in the final image than the initial image, while negative change identifies pixels that became dimmer as a result of less brightness in the final image than the initial state of the image. Since the brighter (higher) the NDVI value the greener the vegetation, positive change indicates gain in vegetation while negative change signifies vegetation loss. In Figure 5, white colour shows area where change did not occur between two epochs considered, purple indicates negative change, and the green colour depicts areas with positive change. Table 5 and Figure 6 quantitatively presented the estimated percentage area for positive, negative and no-change. The area covered by the no-change pixels varies between 42.39 and 59.97% (average 52.64%). Similarly, the area occupied by the negative change vary between 2.36 and 56.73% (range 54.37%) while the positive change has a range value of 44.83% (minimum of 2.79% and maximum of 47.62%). Unlike the no-change pixel, the negative and positive changes are inversely proportional, alternating in reverse order in successive change map. Irregularity in the change status reflects the observed fluctuation in vegetation amount and distribution (Figure 3). It can be inferred from these analyses that the study area experienced vegetation loss by approximately 26.05%.

Time Series Change Detection
The post-classification change detection map between 2015 and 2020 provides map summary of changes that have occurred (Figure 7). In the image, the land cover class represented in white colour indicate the totality of non-vegetated area, other land cover types converted to non-vegetated inclusive. Also, forest and grassland/shrubs are depicted in dark and light green, respectively; while the blue colour indicates open water. In Table 6, statistics of the total area converted to each land cover type is presented. The total land area occupied by forest is estimated to be about 2, 371.13 km 2 , approximately 6.70% of the total landmass. About half (50.07%, ca. 17,770.79 km 2 ) of the State's landmass predominantly comprised grassland/shrubs, while the non-vegetated area is slightly less than half (43.04%, ca. 15,274.85 km 2 ). The water bodies occupied 0.21% of the land mass (ca. 75.68 km 2 ).

DISCUSSION
Determining the vegetation amount, distribution and loss is challenging due to complex interaction of anthropogenic activities, natural (climatic) factors and their spatial heterogeneity (Adubasim et al., 2018). The time series NDVI change detection analysis allows identifying and quantifying change between images of the study area at different times. Both the qualitative and quantitative analysis has shown that there was no significant change in the percentage area coverage for the water bodies and grassland ( Figure 3; Table 4). For instance, the percentage area occupied by water bodies varied between 0.21% and 0.23% while the land/cropland/shrubs, land cover type have a mean of 55.50% (min. of 50.45% and max. of 61.38%) of the total area with a shortfall of about 7.00% between 2015 and 2020. In contrast, there was significant increase and decrease in the non-vegetated and forest, respectively (Figure 4). The non-vegetated area increased from 18.53% in 2015 to 42.56% in 2020 (ca. 25.00% gain), whereas the forest area reduced to 6.78% in 2020 compared to 23.76% in 2015 (ca. 17.00% loss in forest). Investigation into the degree of change using pixel-based and post classification techniques yielded quantifying time series and initial-final state changes (Figures 5 and 8  while the water bodies and grassland/shrubs remained relatively constant. Comparison of the class statistics of the land cover map of 2020 (Table 4) and the change map (Table 6) validates the accuracy of this work. For example, the water bodies produced 0.21% for both maps while the non-vegetated land cover (42.56 and 43.04%), grassland/shrubs (50.45 and 50.07%) and forest (6.78 and 6.68%), respectively. In terms of vegetation dynamics, an interesting pattern emerges that indicates alternating year of increase followed by decrease in vegetation amount (compare Figures 3a and 3b, 3c and 3d, 3e and 3f). This provoked further investigation into possible influence of climatic factors, specifically precipitation and LST. Figure 8 presents the plot of LST and rainfall, both of which are inversely proportional.
High LST was experienced during dry or low rainfall time and reduced with increase in rainfall. Averagely, the State recorded min. and max. mean annual temperature of 20 and 29°C, with precipitation of 243.00 and 1273.00 mm, respectively, typical of tropical savannah climate (Gleixner et al., 2020). Overall, average of seven months of rainfall was recorded, with peak intensity in Aug. and Sep. The amount of rainfall varied in space, time and intensity, while some years had early onset of rainfall that lasted till Nov. (Figures 8b, 8d, and 8e). Figures 8a, 8c and 8f, respectively representing 2015, 2017, and 2020, had late start and early cessation of rainfall with fluctuations in intensity. Obviously, the regional climatic factors (precipitation and SLT) produced similar pattern of variation seen in the vegetation analysis. It is, therefore, valid to infer that regional climate influenced distribution and amount of vegetation (Figure 3), as also observed by others (Schultz and Halpert, 1993;Trenberth and Shea, 2005). Even with the obvious variations in vegetation amount, which often impacts groundwater recharge (Ouyang et al., 2019), reduction in the surface area of water bodies during the investigation was minimal.

CONCLUSION
This study investigates the use of Google earth engine (GEE) using sentinel-2 level C-1 to examine vegetation dynamics and loss between 2015 and 2020. GEE optimized the process of data collection, processing and analysis without the usual tasking computational process involved in the traditional desktop remote sensing image processing. The study has indicated that both anthropogenic and natural factors contributed to variation in vegetation distribution and amount with resultant 25.00% increase in non-vegetated area and 17.00% loss in forest area. Activities such as urbanization, increase in lumbering activities, agricultural practice, bush burning arising from increased temperature and longer months of dry season exposed forested area to non-vegetated land. The NDVI allows evaluating multi temporal measures of time series gain and loss in vegetation. Stability of the aerial coverage of grassland/shrubs to approximately 50.00% of the landmass across the years revealed high exchange between urban and forest area rather than other environmental factors. Further study is necessary to investigate in details, the interaction between vegetation loss and climate change so as to predict the trend of future vegetation dynamics.