Using remote sensing to assess plant health and drought response in game reserves and adjacent farmland overtime in the Eastern Cape, South Africa

The aim of the study described in this article was to investigate the vegetation health and drought response of naturally occurring Albany thicket and neighbouring farmland vegetation, that appears in an area of the Eastern Cape, South Africa. Google Earth Engine was used to manipulate Landsat 5, 7 and 8 datasets to produce a 30-year temporal dataset, after which the normalized difference vegetation index (NDVI) and the normalized difference water index (NDWI) were then applied to create a time series analysis. The Mann-Kendall and Spearman correlation statistical tests were used on the time series to observe trends and correlations between the NDVI and the NDWI datasets. The Spearman correlation test results showed that there were high correlations between the NDVI and the NDWI datasets (all above 0.805). Furthermore, the Man-Kendall test showed that all the datasets had positively increasing trends, while the NDVI datasets all had monotonic positive trends. Large differences in the NDVI and the NDWI were seen for the different vegetation types during times of drought, and farmland was the most severely affected with an average of 19% decrease in the NDVI and an average of 71% decrease in the NDWI.


Introduction
Due to exponential human population growth, natural land areas have been cleared or transformed from their original state for agricultural or urban purposes. Moreover, rapid urban expansion has been identified as placing extreme pressure on natural environments and resources (Oyedotun, 2018).
Urban areas expand into natural vegetation and, owing to the increase in population, more agricultural land is needed to supply food for the residents (Pandian, Rajagopal, Sakthivel, & Amrutha, 2014).
Understanding and monitoring ecosystems can be impeded by evolutionary changes, such as past land cover changes or stress from droughts (Lausch, et al., 2016). Knowing past evolutionary changes is essential for understanding current ecosystems' biodiversity (Lausch, et al., 2018). The use of remotely sensed data is an ideal method of gathering data on past events utilizing tools such as a time series analysis.
The study focused on natural and agricultural land classes and explored how these two classes have changed over time. Conservation efforts have seen agricultural land reverted to a natural area by becoming a game reserve. This study also looked at the possibilities of using remote sensing for monitoring variances in vegetation health and drought tolerances between game reserves and adjacent farmland. The objectives of the study included assessing the quality of the vegetation of two-game reserves that bordered each other. Moreover, the study observed whether any health and drought tolerance changes could be seen over the thirty years since the reserves were started and the land-use activities drastically changed. The results of this observation were compared to the farmland that neighboured the two reserves.
Landsat 5, 7 and 8 aerial imagery datasets were merged to create the temporal analysis. Results were obtained from the Google Earth Engine platform by applying the normalized differential vegetation index (NDVI) and the normalized differential water index (NDWI). The availability of remotely sensed data has allowed for many vegetation related studies to be completed (Guo, Li, He, Xu, & Jin, 2018). Time series analysis has regularly been used to acquire results from extensive temporal studies, and many methods have been developed to interpret the results (Neeti & Eastman, 2011). In this study, the Mann-Kendall (MK) trend test and Spearman correlation were used to assess the time series results.
The stress of vegetation can be caused by soil nutrient loss, soil toxicity, lack of solar radiation, drought (Luus & Kelly, 2008), soil conditions and other factors affecting it. Vegetation indices can calculate vegetation stress, Numata, et al. (2007) discovered that the vegetation index decreases during times of vegetation stress and increases during less stressful times. Other studies have also confirmed that vegetation indices are successful in estimating productivity and phenology over time (Beck, Jonsson, Hogdas, & Karlsen, 2007) (White, Nemani, Thornton, & Running, 2002). Drought stress can be measured through the reading of net primary production (NPP) according to a study done by Asner, et al. (2004). This is done by modelling variables such as chlorophyll absorption, temperature, and the moisture content of a plant. Similar results were also obtained by Jensen (2005).
The most common index used for analysing vegetation is the NDVI, which is based on differences between red and near-infrared (NIR) reflected by vegetation (Tucker, 1979). The NDVI can estimate the state and changes of vegetation due to the index being related to photosynthetically active radiation (Ivanova, Kovalev, Yakubailik, & Soukhovolsky, 2019). A high NDVI shows that foliage has increased photosynthetic capacity, while lower values show a reduction in photosynthesis, which may result in leaf wilting (plant stress) (Gu, et al., 2008).
The NDWI is sensitive to changes in water content within leaves. This is done by the index using the NIR and shortwave infrared (SWIR) bands, which allows the detection of changes in water content and spongy mesophyll (Gao, 1996). SWIR is responsible for the detection of vegetation water content and spongy mesophyll structure, while the NIR value is determined by the leaf structure and dry matter content. By combining the two bands, various inclusions from the leaf structure are removed, allowing for a more accurate result (Ceccato, Flasse, Tarantola, Jacquemond, & Gregoire, 2001). Hazaymeh and Hassan, (2016) have found that data quality becomes uncertain over very thinly vegetated areas (Ghulam, Li, & Qin, 2008).

Study area
Sibuya Game Reserve, a portion of the Kariega Game Reserve, as well as the immediate surrounding farmland formed part of the study area, which is on the Southern Coast of the Eastern Cape, South Africa (See Figure 1). Sibuya Game Reserve was started in 2003 when three farmers merged their land, creating a reserve of 2200 hectares. The section of Kariega Game Reserve that was included in the study surrounds the upper reaches of the Kariega River and borders Sibuya Game Reserve. The area, which is around 1500 hectares, was formed in 1989; natural fauna was introduced to the area in 1991. In both reserves, the fields where farmland was cleared are still visible and are now open grassy plains.
The Kariega and Kasouga rivers flow through the study area, which has led to the topography of the area consisting of valleys and flat river flood plains. The region predominately receives rainfall during autumn and winter months from passing cold fronts, but large electrical storms infrequently bring rain in the summer months. The Cape Supergroup forms the underlying geology, specifically the Witteberg and Bokkeveld Groups. Lithologies of these groups include sandstone, quartzitic sandstone, siltstone, claystone, mudstone and shale (Johnson & Thamm, 2009). The clay-rich soil that is predominant across the area is a common soil type in the Eastern Cape. It was eroded from the shales and mudstones of the Beaufort Group, which is part of the Karoo Supergroup (Matike, Ekosse, & Ngole, 2011). The indigenous flora consists primarily of Albany thicket, specifically the Kowie thicket (Hoare, et al., 2006); however, small amounts of fynbos are found in nutrient-poor sandy loams (Hodson, Vijver, & Peijnenberg, 2011). Thicket species thrive in clay-rich soil (Hoare, et al., 2006) and are concentrated in the valleys. However, they do also occur in areas of moderate topography where clearing for farmland has not happened.
Current farming in the study area is either cattle or cropland. Both agricultural activities have led to the land being cleared and the soil degraded, either by overgrazing or by ploughing. In the valleys that run through farmland, no clearing occurred due to the uneven topography. The natural vegetation has been left alone, and some small animals are still present. In the study, these areas were looked at separately and called natural areas. The different sections of the study area (Sibuya Game Reserve, Kariega Game Reserve, the farmland, and the natural areas) are visible in Figure 1.

Methodology
In the study, a sizeable temporal range of thirty years (1990-2020) was used to obtain results and identify patterns of change. Landsat 5, 7 and 8 satellite imagery were used for the study. Google Earth Engine (GEE) software was used for acquiring, editing and pre-processing of Landsat data for this study; GEE has been developed for educational, non-profit use by Google. Processes to datasets are done using JavaScript, which enables a wide variety of functions to be applied to the data (Gorelick, et al., 2017). The three datasets had to be harmonised, as Landsat 8 wavelengths of bands differ slightly from the wavelengths of Landsat 5 and 7 datasets (Roy, Kovalskyy, Zhang, Vermote, & Yan, 2016), the harmonization of the bands correct for slight differences in spectral characteristics between Landsat 8 and 7/5. Braaten (2019) has developed a script for harmonising the datasets; the code is open access and free to use. Further prepossessing includes masking out clouds and their shadows from datasets, as they can significantly affect the accuracy of the results. Landsat datasets consist of a band called "pixel_qa", which identifies clouds and their shadows in pixels and sets the pixel to a value of zero (Braaten, 2019). Zhu, et al., (2015) have developed an open-source cloud masking function, which uses "cloudShadowBitMask", "cloudsBitMask" and "pixel_qa" to mask clouds and their shadows out of images that have been used for cloud masking. The harmonisation, cloud masking and indices were written as functions, but before they could be applied, the datasets had to be clipped to the area of interest and the date range specified, and the datasets merged. Due to the merging of the three datasets, there was a large amount of overlap and some months had up to four indices values A function was created to calculate the average indices value for every year in the data range.
After renaming of bands, harmonization, cloud masking, applying spectral indices (NDVI/NDWI), setting the date range and cropping to the area of interest, and reducing the dimensionality of data to yearly values, the raw data was worked into two result formats: time-series graphs and spectral maps.
Time-series results have been used for detecting growing periods, greening and browning trends, vegetation activity changes and vegetation degradation over time (Guo, Li, He, Xu, & Jin, 2018). Forkel, et al. (2013) conducted a study comparing many of the available types of analysis on the NDVI time series. In this study, the Mann-Kendall Trend Test and the Spearman correlation were the focus of the analysis.
The Mann-Kendall test is non-parametric and detects monotonic trends in data series (de Jong, de Bruin, de Wit, Schaepman, & Dent, 2011) (Mbatha & Xulu, 2018). Non-parametric tests may be used with missing data, as the validity does not depend on datasets being normally distributed (de Jong, de Bruin, de Wit, Schaepman, & Dent, 2011) and is not affected by outliers, making the test ideal for environmental, climatic and hydrological data (Pohlert, 2020). Time series datasets derived from vegetation indices often do not meet parametric standards like normality and homoscedasticity (Guo, Li, He, Xu, & Jin, 2018).
The Mann-Kendall statistic (S') value is initially set to 0 at the start of testing and represents no trend. If data from one period in the time series are greater than the initial data, S' is increased by a value of 1. If the data are less, then S' is decreased by 1 (Khambhammettu, 2005). This procedure is done for the entire time series, and the final value of S' will indicate an increasing trend (when S' is positive), a decreasing trend (when S' is negative) or no trend (when S' is zero). The probability of S', and thus the significance of the trend, is represented by the variance of S' (Var(S') (Guo, Li, He, Xu, & Jin, 2018). The p-value indicates if a monotonic trend is present in the dataset. If p is <0.05, there is a monotonic trend, and if p is >0.05, there is a non-monotonic trend (Karmeshu, 2012) (Guo, Li, He, Xu, & Jin, 2018). However, if the p-value is smaller than the alpha value, the null hypothesis is rejected, indicating that there is a trend present in the data, and the results are statistically significant (Karmeshu, 2012).
For the correlation test, the Spearman correlation coefficient was used. This measures the degree of correlation between variables (Rezaee, Aliabadi, Dorestani, & Rezaee, 2020). A correlation between variables is the measurement of the strength of the relationship between variables and ascertains whether they are positively or negatively related (Obilor & Amadi, 2018). The Spearman correlation coefficient is based on the ranks of the observations and is thus non-parametric; it defines the relationship between two or more variables using a monotonic function (Rezaee, Aliabadi, Dorestani, & Rezaee, 2020). The results of the Spearman correlation are numerical and range from -1 to +1. The value -1 indicates a perfect negative relationship, +1 a perfect positive relationship and 0 indicates no correlation between the variables (Rezaee, Aliabadi, Dorestani, & Rezaee, 2020).

NDVI analysis and results
The NDVI was developed to indicate plant health conditions and vegetation density. The results from an NDVI analysis are rated between -1 to 1. The values -1 to 0 indicate that the vegetation is dead or very unhealthy and near death, 0 to 0.33 shows unhealthy/stressed vegetation, 0.33 to 0.66 represents moderately healthy vegetation and 0.66 to 1 signifies very healthy and unstressed vegetation (Ivanova, Kovalev, Yakubailik, & Soukhovolsky, 2019).
In the study, the time-series graph of NDVI values was most valuable for visualising trends in the data over the study period. The yearly NDVI averages were used to plot a time series graph for the three areas; this is presented in Figure 2. The graph showed that the reserve, the natural area, and the farmland all followed a similar pattern, which could be expected, as these areas all experience the same climatic conditions; however, the graph also showed clearly how they differ. The natural region's yearly line was above the other two areas for the whole period; the farmland yearly line was below the other two. The reserve's yearly line was sandwiched between the natural area and the farmland's line, which could have been due to the reserve having both old farmland and wild fauna/flora. A spectral map of each region was produced for the study; the images were an average of NDVI reflectance between 1990 and 2020. These images, which are displayed in Figure 3, showed the spatial patterns of NDVI across the areas. Typically, NDVI is measured between -1 and 1. However, the scale was adjusted for these images to 0 to 1, as the images showed no data below zero. In the images, water was detected as being very unhealthy due to it absorbing NIR and red. Therefore

NDWI analysis and results
The results of the NDWI are measured from -1 to 1. The value -1 represents vegetation with little to no water content, 0 represents a moderate amount of water content, while 1 represents vegetation with high water content. The NDWI formula was applied to the datasets in the form of a function, producing an NDWI layer. This layer and the associated data were used for the analysis.
Just as the NDVI time series were shown in Figure 2 above, the yearly averages for the three areas from 1990 to 2020 are presented in Figure 4. The three lines in the graph each had similar patterns, with highs and lows occurring at the same time due to the areas being in proximity to each other and experiencing the same climatic conditions. The reserve's line was sandwiched between the natural area line (above the reserve line) and the farmland line (below the reserve line), which was due to the reserve having both old farmland and natural bush. In the year 2015, a drop in the graph's line was visible; this was caused by the start of an extended drought in the area. Throughout the time series, other notable decreases could be observed, which also correlated to drought conditions or low rainfall years. Figure 4: Yearly Averages of NDWI over the three areas from 1990 to 2020.
The spectral maps for the NDWI followed the same pattern as those for the NDVI, in that each image depicted the yearly averages for the three areas from 1990 to 2020. The NDWI spectral images are displayed in Figure 5. The NDWI is measured from -1 to 1, but as the values in the images all occurred between -0.5 and 0.5, this was set as the scale to be used. A colour ramp of three colours (red, yellow, and blue) was used to represent the scale: -0.5 was symbolised by red, 0 was represented by yellow and 0.5 was represented by blue. Very few red areas were displayed on the spectral maps; only small portions of the reserve and the farmland showed red. For the most part, the yellow colour was widespread throughout the farmland and reserve areas. However, in the valleys in the reserve and the natural areas, blue was far more prominent, as this was where natural bush occurred.

Mann-Kendall trend test
The Mann-Kendall test was done on both the NDVI and NDWI datasets at a 95% significance level, and a continuity correction was performed. The Mann-Kendall results for the NDVI data are seen in Table 1. The Mann-Kendall statistic was positive for the three areas and had a high value, which showed an overall increasing trend. The p-value was lower than the alpha value for the three regions, which further indicated that the null hypothesis was rejected. Each area had a p-value less than <0.05, which showed that a monotonic trend could be identified. The results from the Mann-Kendall NDWI test results that are presented in Table 2 differed significantly from the NDVI results. The areas all had a positively increasing NDWI trend from the Mann-Kendall statistic, but the natural region had a much stronger trend than the reserve and the farmland. Each region rejected the null hypotheses by having a p-value lower than the alpha value.
However, only the natural land had a monotonic trend as the reserve and the farmland both obtained a p-value higher than 0.05.

Spearman correlation test
The results from the correlation are measured between -1 and +1, with -1 specifying a perfect negative relationship, +1 a perfect positive relationship and 0 no correlation between variables. The Spearman correlation test was run on a combination of the NDVI and NDWI datasets, the results of which are displayed in a matrix table in Table 3. The NDVI and NDWI for each area showed a strong positive correlation to each other, with all being above 0.8, while correlations dropped when comparing different regions to others. However, a positive correlation was still seen across the table. As stated above, a strong positive correlation was seen when comparing NDVI and NDWI between the same areas. The three scatter plots presented in Figure 6 showed the NDVI and NDWI correlation between each area. The densely packed data points forming a linear pattern indicated a strong correlation. Multiple outliers could still be seen; the natural area had the most severe outliers, which was likely responsible for the lowest correlation value of the three areas.

Trend and correlations tests
The Mann-Kendall test was used for the analysis of trends in both the NDVI and NDWI timeseries datasets, results from the test were displayed above in Table 1 for NDVI and Table 2 for NDWI.
The NDVI S' results all showed a strong, increasing trend, with the natural area having the highest (2424), and the reserve area having the weakest (1859) of the three areas. Var(S') was very high for the NDVI areas; 103251-103272 was the range, which indicated that the S' results were highly probable, and the data was significant. The NDWI S' results were weaker than the NDVI results; the natural land had the highest S' value (1130) The p-value is a valuable output as it indicates whether a dataset is significant, whether a trend is present and whether this trend is monotonic or not (Guo, Li, He, Xu, & Jin, 2018). For both the NDVI and NDWI, the null hypothesis was rejected, the p values for NDVI were minimal (<0.0001), and the p values for the NDWI ranged from 0.002 -0.2849. All p-values were positive, which showed a positive trend across all the datasets. Rejecting the null hypothesis indicates that the datasets are statistically significant and that a trend is present in the datasets. To determine if the trend is monotonic or non-monotonic, the p-value must be smaller than 0.05. From the NDVI results, all three areas showed a monotonic trend, and the NDWI results showed that only the natural area had a monotonic trend. In contrast, the reserve and the farmland had p values greater than 0.05 and therefore showed a non-monotonic trend. A positive monotonic trend seen in the results of a test reveals an entirely non-decreasing trend (Royden & Fitzpatrick, 2010). A non-monotonic trend shows that the trend within seasons increases and decreases. This trend occurred in the case of the NDWI, which showed that the two areas (the reserve and the farmland) responded poorly to drought conditions. In contrast, the natural area was able to maintain a monotonic trend throughout the time series. The S' results from the NDWI supported the p-value results and determined the farmland and the reserve as having a weak trend.
The Spearman correlation test was the ideal correlation coefficient to use, as it can work with nonparametric data and uses a monatomic function to measure data (Rezaee, Aliabadi, Dorestani, & Rezaee, 2020). The results were displayed in Table 3 and the scatter charts in Figure 5. The correlation between the reserve's NDVI and NDWI was the highest at 0.817, which indicated a robust correlation A study conducted by Dubovyk et al (2015) on the use of medium resolution satellites for monitoring vegetation dynamics across southern Africa used the Mann-Kendall trend test and the Spearman correlation coefficient for detecting trends and correlations. Results from the study showed that the two statistical tests worked well together. Statistically significant positive trends were obtained for the southern area of South Africa, which was the area of the study described in this article. The results of the study conducted by Dubovyk et al (2015) support those obtained from this study.

NDVI and NDWI for monitoring plant health and drought response
In the study, significant seasonality trends were seen in the NDVI and NDWI time-series datasets, the results of which were displayed in Figures 2 and 4 above as time-series graphs. The two figures showed the NDVI and NDWI reflectance values for all three areas. A general examination revealed that the areas had similar patterns, increases and decreases occurring at similar times. The similar NDVI and NDWI patterns were due to the strong correlation seen between the NDVI and NDWI in the Spearman Correlation results. In a study done by Karamihalaki et al. (2016), the NDVI and NDWI were used to monitor drought effects on a forest. From their results, a strong correlation between the two indices was observed, and the NDWI showed a strong correlation between precipitation values.
In the study described in this article, the time series showed that the natural area had the highest NDVI and NDWI values, while the farmland had the lowest NDVI and NDWI values of the three regions. The reserve has both old farmland and natural vegetation, which averaged out the NDVI and NDWI values. Thus, the reserve showed trend patterns that were similar to those seen in the natural area and farmland. From the NDVI time-series graph, it was clear that the natural area had the healthiest vegetation of the three areas, with the highest value (0.743 NDVI) being recorded for the year 2015, while the farmland obtained 0.685 NDVI and reserve 0.727 NDVI. The same results were seen for NDWI: the natural area obtained 0.276, the reserve 0.259 and the farmland 0.185.
The shape and pattern of the time series graphs indicated a history of drought and plant stress within the area from 1990 to 1992. Very low NDVI and NDWI values were recorded, which was the case again in 1996, 1998-2000, 2004, 2008-2010 and 2016 to the present. However, between 2019 and 2020, an increase was seen in the NDVI and the NDWI when drought conditions still affected the area. There was no drought in 1996 and 2004, but these were relatively dry years. The average decrease in NDVI during drought conditions (excluding 1996 and 2004) was 16% for the natural area and 19% for both the reserve and farmland areas. The average NDWI decrease during the drought mentioned above was 53% for the reserve area, 44% for the natural area and 71% for the farmland.
During times of drought, the natural area was the least affected of the three areas. The NDWI decrease (and thus the moisture content) was likely the lowest due to the thick vegetation canopy protecting the lower portion of the trees and the soil below. In addition, natural flora had evolved to survive in the environment. These factors were likely responsible for the low stress/health decrease during drought. The farmland had a very high reduction in the NDWI. This could have been due to many factors, such as the vegetation not being able to tolerate drought conditions, over-grazing causing soil exposure thus increasing evaporation, the foliage of the crop being susceptible to transpiration or the harvesting of plant biomass negatively affecting water contents. With such a high NDWI decrease value, the low NDVI decrease value was unexpected. However, possible farming techniques might have been in place to keep plants healthy during drought periods. The reserve's NDVI and NDWI decrease was slightly higher than in natural areas. This could be partially explained by the reserve having natural vegetation and old farmland. However, the reserve area had a significant amount of old farmland, and it is plausible that due to the decrease in over-grazing, which gave the flora a long time to adapt to the climate, and due to the halting of harvesting, the old farmland showed recovery and resilience to drought conditions. The spectral maps presented in Figures 2 and 4 above provided support to the results from the time series. It could be seen on the NDWI images that the farmland area had a poor distribution of flora that could retain water, thus during drought conditions, the farmland had a decline in water content, which was also seen in the NDWI time series results. Natural flora has a high-water retaining ability, and this could be seen in the NDWI images for the reserve and the natural area. The NDVI spectral images also represented the natural vegetation as healthy green vegetation. When comparing old farmland areas on the reserve to current farmland, the old farmland had slightly healthier vegetation.
However, the NDVI colours were similar; in the NDWI images, the fields also showed signs of higher water content.

Concluding remarks
The results from the Mann-Kendall trend analysis showed that the NDVI results carried more significance than the NDWI results, while high correlation trends were seen between NDVI and NDWI regarding the three land classes. The NDVI and NDWI results showed how severely the land change and subsequent overgrazing, poor land management, bush clearing, and other farming techniques had led to the land degradation seen across old and current farmlands within the study area. Natural bush in the area had evolved to deal with droughts, to a point, and other stress-inducing climatic conditions, which is evident from this class having the highest NDVI & NDWI results. The results showed this class to be more resistant to droughts and less affected during a time of decreased rainfall, with NDVI decreasing by 16% and NDWI decreasing 44% during past droughts. As the vegetation species varied between the farmland and the natural area, it was expected to see a difference in NDVI and NDWI, during times of drought both areas NDVI decreased by 19%.
However, the NDWI during times of drought for the farmland (average decrease of 71%) and for the natural area, which showed a decline of 44%, revealed the massive difference between the vegetation types. The results from assessing the quality of the vegetation of two-game reserves that bordered each other was proof that the land use activities greatly affected the natural vegetation which occurred.
Furthermore, the study observed that health and drought tolerance changes could be seen over the thirty years since the reserves were started and the land-use activities drastically changed.