Woody plant cover, the area of the vertical projection of woody plants (trees, shrubs, and bushes), plays an important role in the structure and function of savanna ecosystems and is needed by the savanna modeling community. Recent problems facing savanna ecosystems such as woody plant encroachment and subsequent habitat fragmentation further underscore the relevance of regional-scale and even larger-scale woody plant cover mapping. The mixture of woody plants and herbaceous vegetation in savanna landscapes lends woody plant cover mapping to fractional representation. This study endeavors to develop a simple and reliable approach for fractional woody plant cover mapping in savanna ecosystems. It was tested in the savanna of central Texas, which features a wide woody plant density gradation. A multiple linear regression model was calibrated between orthophoto-based fractional woody plant cover and metrics derived from time series MODIS products of surface reflectance (MOD09A1) and fraction of photosynthetically active radiation (MOD15A2H). By applying this model, woody plant cover was extrapolated to Texas savanna at MODIS scale (500 m). Validation suggests a mean absolute error of 0.098 and an R-squared value of 0.60. This study demonstrates a potential approach for woody plant cover mapping in other savanna ecosystems of the world. It also highlights the utility of time series MODIS products in savanna woody plant cover estimation.
Woody plant cover, the area of the vertical projection of woody plants (trees, shrubs, and bushes), has huge impacts on landscape connectivity, species diversity, and livestock production of savanna ecosystems (Alofs and Fowler 2010; Ratajczak et al. 2012; Anadón et al. 2014). Large-scale woody plant cover data has been needed in woody plant carbon estimation, potential woody cover modeling, and vegetation dynamics analysis (Sankaran et al. 2005, 2008; González-Roglich and Swenson 2016; Yang et al. 2016). Recently the global phenomenon of woody plant encroachment in savanna ecosystems further underscores the need of woody plant cover map of large scale for the encroachment monitoring and modeling, and for management decision-making (Hudak and Wessman 1998; Van Auken 2009; González 2010; Yang et al. 2016; Stevens et al. 2017).
The increasing availability of remote sensing data provides us an opportunity to map woody plant cover over large savanna areas. The mixture of woody plants and herbaceous vegetation lends savanna woody plant cover mapping to fractional representation that estimates landscape components at subpixel level (DeFries et al. 1995). It is especially true with coarse-resolution remote sensing data. The fractional representation is more accurate and bears greater potential in research and application (Hansen et al. 2003). For instance, it can distinguish different woody plant densities (DeFries et al. 2000; Hansen et al. 2002), which are highly susceptible to resource availability (water, nutrient) and disturbance regime (fire, herbivory) (Sankaran et al. 2005). It can also be used to model savanna structure (Lefsky 2010). In addition, successive fractional woody plant cover maps are capable of monitoring woody plant cover change over time.
A number of global-scale tree cover products covering savanna ecosystems have been developed (Hansen et al. 2013; Sexton et al. 2013; Kobayashi et al. 2016; DiMiceli et al. 2017). However, these products are targeted at trees above a certain height and do not explicitly consider short woody plants that are common in savanna ecosystems (Yang and Crews 2019a). Particularly poor performance of these products has been reported in sparse tree areas (White et al. 2005; Montesano et al. 2009). Therefore, their ability in reflecting savanna woody plant cover has been a concern for researchers.
Recently there have been many efforts mapping woody plant cover in savanna ecosystems. For instance, Urbazaev et al. (2015) calibrated synthetic aperture radar (SAR)-based woody plant cover models with high-resolution airborne lidar data for a small savanna area in southern Africa. However, this method may not be applicable for large savanna areas due to the limited coverage and availability of high-resolution lidar data. Brandt et al. (2016) mapped the fractional woody plant cover of Sahelian drylands of Africa with Moderate Resolution Imaging Spectroradiometer (MODIS) and Satellite pour l’Observation de la Terre (SPOT)–Vegetation (VGT) data. They mainly take advantage of the different phenophases between woody plants and herbaceous vegetation. However, this approach requires deep prior knowledge about the vegetation phenology of the study area. Moreover, the phenology of woody plants or herbaceous vegetation in other savanna areas may not be as uniform as that in the Sahel area. Yang and Crews (2019b) mapped the fractional woody plant cover of Texas savanna with Web-Enabled Landsat Data (WELD) and classification and regression trees (CART). However, this study involves a huge amount of high quality training data, which is usually not easily available in other savanna ecosystems. Hence, this study aims to develop a more applicable and reliable approach for savanna woody plant cover mapping.
A wealth of MODIS products has been developed since 2000, covering a broad range of variables such as net primary production, snow cover, and land surface temperature (Turner et al. 2006; Wang et al. 2009; Zeng et al. 2015). These publicly available products bear great potential in savanna woody plant cover mapping. First of all, while their global coverage ensures spatial comparability, their high temporal frequency (e.g., 8 day) can minimize the effect of cloud contamination. These characteristics make them suitable for large-scale study (Townshend et al. 1994; Schwarz and Zimmermann 2005; Kobayashi et al. 2016). Moreover, the time series MODIS products can cover different phenophases of woody plants and herbaceous vegetation. Specifically, derived seasonal or annual metrics are able to capture the characteristic points of their phenology cycles (Hansen et al. 2002). Furthermore, the visibility of major plant morphological traits (e.g., structure) and human imprint (e.g., cultivation) at MODIS scale underscores the value of MODIS products in such study (Hansen et al. 2003; Vintrou et al. 2014). Therefore, this study will explore the utility of MODIS products in woody plant cover mapping in savanna ecosystems. Specifically, it will develop an appropriate statistical model linking reference woody plant cover to MODIS products and their derivatives.
2. Study area
The proposed method of this study was tested in the savanna of central Texas (Figure 1). This savanna area is underlain by a thick and relatively flat bedrock of hard early Cretaceous limestone (Schmid 1969). Mean annual precipitation ranges from ~530 to ~950 mm in the west–east direction. In recent decades, this savanna area experienced substantial woody plant encroachment, resulting in habitat fragmentation, native herbaceous species loss, and biogeochemical changes (Osborn and Witkowski 1974; Archer 1989; Asner et al. 2003; Fowler and Simmons 2009; Alofs and Fowler 2010, 2013; Creamer et al. 2013). Knowing the abundance and distribution of woody plants across this area is important for management decision-making and the savanna study (Archer 1990; Ansley et al. 2001; Hughes et al. 2006).
Junipers (Juniperus), mesquites (Prosopis), and oaks (Quercus) are the most common woody plants across Texas savanna (Appel 1995; Ansley et al. 2001; Lyons et al. 2009). These evergreen and deciduous woody plants differ in phenology cycles, but are interspersed with each other. In addition, huge variation has been observed in physiognomy, composition, and structure along the precipitation gradient of this area (González 2010). All these characteristics make it a challenge to accurately map woody plant cover in this area by satellite remote sensing.
3. Materials and methods
3.1. Conceptual approach
The concept of this study is to develop an appropriate statistical model of the relationship between woody plant cover that is based on “orthophotos” and metrics that are based on MODIS products (Figure 2). The MODIS products include 8-day surface reflectance and fraction of photosynthetically active radiation (FPAR). Then, by applying this statistical model and the metric layers, it will extrapolate woody plant cover to Texas savanna at 500-m scale. The woody plant cover “result map” will be validated with test data of orthophoto-based woody plant cover.
3.2. Orthophoto-based woody plant cover measurement
Training data and test data were derived from 1-m-resolution digital orthophotos of 2012, from the National Agriculture Imagery Program (NAIP). The orthophotos were taken in growing season and have four spectral bands—red, green, blue, and near-infrared. This high-resolution and leaf-on imagery is capable of distinguishing short woody plants from the background.
In response to the principle of representativeness of training data in land-cover mapping, we purposefully prepared a large amount of training data across the study area (Foody et al. 2006). First, a total of 1200 MODIS pixels (500 m) were randomly sampled from a MODIS surface reflectance layer (MOD09A1) of Texas savanna. The corresponding orthophoto of each MODIS pixel was classified into 10 thematic classes using the “Iso Cluster Unsupervised Classification” tool in ArcGIS. Each thematic class was then attributed as “woody plant” or “the others” (e.g., herbaceous vegetation or bare ground), through visual inspection of the natural-color and color-infrared composites (Figure 3).
Second, to test the accuracy of the above classifications, 100 random pixels (1 m) were sampled across the orthophoto of a MODIS pixel. The reference status of each pixel (woody vs nonwoody) was obtained by checking the natural-color and color-infrared imagery. These reference data were then compared to the above classification result, suggesting an overall accuracy of 87%, an omission rate of 8%, and a commission rate of 5%.
Last, the fractional woody plant cover in each MODIS pixel was calculated as the proportion of “woody plant” pixels. A random portion of 840 (70%) orthophoto-based woody plant cover measurements were used for model calibration. The remaining 360 (30%) measurements were reserved as test data for validation purpose.
3.3. MODIS products and derived metrics
The version-6 MODIS/Terra 8-day surface reflectance (MOD09A1) and FPAR (MOD15A2H) products were used in this study. Both products have a spatial resolution of 500 m. The product MOD09A1 contains surface reflectance estimates for seven spectral bands, namely, band 1 (620–670 nm), band 2 (841–876 nm), band 3 (459–479 nm), band 4 (545–565 nm), band 5 (1230–1250 nm), band 6 (1628–1652 nm), and band 7 (2105–2155 nm). This product is a composite of 8-day observations with selection criteria of high observation coverage, low view angle, and minimum cloud and aerosol loading (Vermote 2015). The Fpar_500m layer of the product MOD15A2H estimates the fraction of incident photosynthetically active radiation (400–700 nm) absorbed by green vegetation (Myneni et al. 2015).
All cloud-free MOD09A1 and MOD15A2H of 2012 were downloaded for the study area from Land Processes Distributed Active Archive Center (LP DAAC). Annual maximum and minimum reflectance layers were derived and named in the convention of B1max–B7max and B1min–B7min. Each FPAR layer is named by the start yearday of its acquisition period (e.g., FPAR001 or FPAR249). Annual maximum and minimum FPAR layers (FPARmax and FPARmin) were also developed. All of these layers were prepared as predictors of fractional woody plant cover.
3.4. Model development
Initially, we explored the relationship between the training data of orthophoto-based woody plant cover and each of the predictors. The predictors fall into two categories, namely, surface reflectance and FPAR. As shown by the two typical scatterplots in Figure 4, woody plant cover tends to decrease with surface reflectance and increase with FPAR. In addition, collinearity was found among the predictors of each category. Based on these preliminary analyses, we chose to build a multiple linear regression model to estimate fractional woody plant cover.
4.1. Multiple linear regression model
To develop a minimal adequate model, we first identified the predictor within each category with most explanatory power of woody plant cover. It was found to be B6min and FPAR249. Second, we developed a bivariate linear regression model between woody plant cover and B6min and FPAR249. Third, we added each of the remaining variables in both categories into the bivariate model in a stepwise procedure, and compared the performance of the model before and after adding the variable. It was found that the incorporation of B7min significantly improved the model performance. Last, a multiple linear regression model was developed with software R between fractional woody plant cover and B6min, B7min, and FPAR249 as below:
In this equation, FWC represents fractional woody plant cover. The standard errors of the above parameters from left to right are 0.305 19, 0.072 73, 0.392 15, and 0.054 01. All of the p significance values are below 0.001. The relative importance of each variable in the regression model was calculated with the function calc.relimp of the R software package relaimpo. Among them, B6min has most power in predicting woody plant cover (42.9%), followed by FPAR249 (31.5%) and B7min (25.6%).
4.2. Extrapolated woody plant cover map
With application of the above multiple linear regression model and the three predictor layers, we extrapolated fractional woody plant cover to Texas savanna at MODIS scale (500 m) (Figure 5a). Small negative predictions were set to zero. Water bodies (the inside white areas) are masked by excluding areas of filled FPAR values (249–255) (Figure 5a). This result woody plant cover map is displayed in comparison with the MODIS tree cover map (MOD44B, 250 m) of 2012 (Figure 5b) (DiMiceli et al. 2017).
As illustrated by the common color ramp from light green to dark green, the extrapolated woody plant cover ranges from 0 to 0.90, whereas the MODIS tree cover barely reaches 0.54. In comparison, this result woody plant cover map exhibits a much higher spatial coherency across the study area. According to Figure 5a, the low woody plant cover related to low precipitation is evident in the western region (Yang et al. 2016). The densest woody plant area mainly lies in the southern region of Texas savanna. Part of the eastern region, however, does not have a high woody plant cover as expected by its high precipitation. This is probably due to its adjacency to the metropolitan areas of Austin and San Antonio and related human activity and urban development. On the contrary, the MODIS tree cover product (Figure 5b) only captures part of the dense woody plants in the southern region and fails to reflect the medium to low-level woody plant cover in the rest of the study area.
4.3. Validation of the result woody plant cover map
We evaluated the accuracy of the extrapolated woody plant cover map with the test data of 360 orthophoto-based woody plant cover measurements. According to Piñeiro et al. (2008), we regressed the orthophoto-based woody plant cover on the extrapolated woody plant cover (Figure 6). The regression line (red) and regression details are overlaid on the scatterplot. It is clear that the regression line is very close to the 1:1 line (black). The regression has an R-squared value of 0.60. The extrapolated woody plant cover has a mean absolute error (MAE) of 0.098. It is easy to notice from Figure 6 the slight underestimation trend of the extrapolated woody plant cover in the high end and the slight overestimation trend in the low end.
This study estimated the fractional woody plant cover of Texas savanna, with application of MODIS products and multiple linear regression. The validation suggests a reasonable pixel-level accuracy (Figure 6). The result map is of great value to management decision-making, especially in face of the ongoing woody plant encroachment in Texas savanna (Basant and Wilcox 2017; Zhou et al. 2018). First, it explicitly depicts the abundance and distribution of woody plants across Texas savanna at MODIS scale (500 m). Second, the MODIS scale is appropriate for management practice. Third, this map is useful for follow-up research such as woody plant carbon storage estimation, savanna structure modeling, and savanna dynamics study (Bucini and Hanan 2007; Sankaran et al. 2008; González-Roglich and Swenson 2016; Krofcheck et al. 2016). However, note that this MODIS-scale map may not be suitable for finer-scale savanna studies such as fire-mediated demographic bottleneck modeling and tree–grass interactions modeling (Hanan and Lehmann 2010; Prior et al. 2010).
The continuous woody plant cover values of this study are more realistic and represent an improvement over the discontinuity of those products based on CART (Hanan et al. 2014; Yang and Crews 2019b). This improvement should be attributed to the application of multiple linear regression (MLR) model. Moreover, this study (the MLR model) has a much lower requirement on the amount of training data than those applying a machine-learning approach (Hansen et al. 2003, 2011). Furthermore, as mentioned in section 4.2, the result map of this study can better reflect woody plant cover in Texas savanna than can the MODIS tree cover product (Figure 5). This improvement draws on two aspects. First, the training data of this study incorporates all trees, shrubs and bushes, while the training data of the MODIS tree cover product is limited to trees above 5 m in height (Hansen et al. 2003). Second, the training data of this study is much more representative of Texas savanna, compared to the global training data of the MODIS tree cover product, which was derived from classified Landsat images.
It is interesting to find that savanna woody plant cover shows a strong negative linear relationship with the annual minimum reflectance of band 6 (1628–1652 nm) and band 7 (2105–2155 nm) of MODIS data. It is probably because of the absorption of radiation at these two shortwave infrared bands by woody plants. The fraction of photosynthetically active radiation in early September (FPAR249, i.e., yeardays 249–256) has a strong explanatory power in woody plant cover in Texas savanna. It is possible that herbaceous vegetation in Texas savanna starts senescing in early September while woody plants remain photosynthetically active. However, FPAR from February to April (e.g., FPAR041, FPAR081, and FPAR105) has very minimal and even no explanatory power in the woody plant cover. It is indicated that FPAR is highly susceptible to the phenology cycles of woody plants and herbaceous vegetation. This new finding highlights the greater value of time series remote sensing data than a static snapshot.
Note that, while this study shows good overall accuracy, it does not perform very well at the high end of woody plant cover. The slight underestimation trend of the extrapolated woody plant cover in the high end is probably related to the limited value range of the predictors. That is, there may be a lower bound value in B6min and B7min in dense woody plant area. The FPAR249 may experience saturation beyond certain woody plant density. Like other similar studies (White et al. 2005; Montesano et al. 2009), this study does not perform particularly well at low densities, as evidenced by the slight overestimation trend (Figure 6). It is expected that these shortcomings will be overcome in future studies by applying more advanced statistical methods.
The R-squared value of 0.60 for this 500-m-scale study is comparable to the 0.77 value of a 1000-m-scale study in the Sahelian drylands of Africa (Brandt et al. 2016). However, this study is much more replicable. First, it uses MODIS products (e.g., FPAR, reflectance), which are publicly available and can be candidate predictors to test in other savanna regions. In particular, the derived annual metrics proved their potential in reflecting the characteristic points of phenophases of woody plants (Hansen et al. 2002). Second, it does not involve design and calculation of complex seasonal metrics. Third, the orthophoto-based training data are much easier to access than field based measurements, especially for remote and fenced areas, and can be more representative of the study area. The high spatial resolution imagery in Google Earth is a possible alternative source of training data when the orthophoto is not available (Kobayashi et al. 2016).
Note also that savanna woody plant cover is anticipated to change due to human activity and climate change (Stevens et al. 2017). For instance, savannas are habitats to one-fifth of the human population and thus are highly susceptible to land clearing, road building, and fire suppression (Lamprey and Reid 2004; Archibald et al. 2013). The elevated atmospheric CO2 concentration is supposed to promote woody plant encroachment (Donohue et al. 2013). Changes in future precipitation regime will also affect savanna woody plant growth (Bucini and Hanan 2007; Yang et al. 2016). All of these expected variations in savanna woody plant cover highlight the need for a large-scale woody plant cover map of global savanna ecosystems for monitoring and modeling purposes.
This study developed a simple and reliable approach to map woody plant cover in Texas savanna at MODIS scale (500 m). It provides us an efficient and affordable way to assess the abundance and distribution of woody plants in Texas savanna in a timely manner. In addition, the result map has great value to follow-up research on woody plant carbon estimation and savanna structure modeling. In addition, this study demonstrates a potential way to map woody plant cover for other savanna ecosystems of the world.
Many thanks are given to the editor and reviewers, whose comments and suggestions greatly improved this paper.