This study was undertaken in the GGHNP located in the foothills of the Drakensberg Mountains in eastern Free State Province of South Africa. Golden Gate Highlands National Park lies between latitude 28°30’ and 28°45’ south and longitude 28°30’ and 28°37’ east, on the border between South Africa and Lesotho (Figure 1). The park forms a part of the Maloti Drakensberg Catchment Complex that produces about 50% of the total water supply in South Africa. The area falls under summer rainfall region where the rainy season stretches from September to April with mean annual rainfall of 780 mm. The park is underlain by rock formations representing the upper Karoo Sequence which is intruded by dolerite dykes and sills.

The vegetation of the park is predominantly grassland, with a very small percentage of woodland communities along the rivers and in sheltered places that are protected from fire. Mucina and Rutherford (eds. 2006) identified four dominant vegetation units associated with the study area, namely Northern Drakensberg Highland Grassland (Gd 5), Lesotho Highland Basalt Grassland (Gd 8), Eastern Free State Sandy Grassland (Gm 4), and Basotho Montane Shrubland (Gm 5).

The high-resolution imagery provided via Google Earth has been increasingly used in scientific research (Fortin & Edwards 2001; Madin, Madin & Booth 2011; Pringle 2010) to assist in the selection of field sampling locations (Tang et al. 2010). The woodland polygons were delineated using high-resolution images within Google Earth mainly from Maxar technologies responsible for WorldView 1-3, and Geoeye with less than 1-m resolution from 2015 to 2020. Because the study area is a mountainous environment, time-slide in Google Earth was used to access historical image archives to offset shadow effects. The woodland community polygons were converted from Google Earth files to shapefile and topologically corrected in ArcGIS version 10.5.1. Fieldwork was also conducted to ensure that the woodland community polygons are adjusted and that the boundaries are correctly delineated for further analysis. Purposive sampling was undertaken using the delineated woodland polygons to verify mapped and unmapped woodland communities covered by the mountain shadows. The woodland community polygons data were used as training samples for statistical analysis (Estes et al. 2012).

Different environmental data sets were used to determine how they influence the spatial distribution of the woodland communities across the landscape. The environmental variables were downloaded from different databases sources (Table 1). Soil chemical and physical properties were one of the environmental variables used. Soil chemical variables included pH, nitrogen (N), soil organic carbon, cation exchange capacity (CEC), whilst soil physical properties included silt, coarse fragments, bulk density, sand, and clay content (Hengl et al. 2017). Shuttle Radar Topography Mission (SRTM) Digital Elevation Models (DEM) was also acquired (https://www2.jpl.nasa.gov/srtm/). Slope and aspect were computed from DEM using ArcGIS 10.5.1.

Topographical variables such as slope and aspect were computed using ArcGIS version 10.5.1. Topographic Position Index (TPI), Topographic Wetness Index (TWI), Topographic Ruggedness Index (TRI) and Topographic Roughness were computed using QGIS version 3.0 (Guisan et al. 1999; Mokarram, Roshan & Negahban 2015; Muddarisna et al. 2020; Radułaa, Szymura & Szymura 2018; Riley, DeGloria & Elliot 1999). The TPI computes the differences of elevation at a specific pixel (central pixel) and the average elevation around it within a defined radius (Gallant & Wilson 2000; Weiss 2001). The elevation is an input to the computations of the TPI. The TWI is one of the hydrologically based topographic index. Topographic Wetness Index describes the tendency of a cell to accumulate water (Gruber & Peckham 2009), and can be used as an indicator of soil moisture content. The input variable to calculate TWI is the slope as computed from the elevation data. The TRI was computed according to Riley et al. (1999) providing summary of change in elevation over 3 × 3 pixel window. It computes the terrain heterogeneity measurement for every location, and each pixel contains the difference in elevation from a centre pixel and the eight pixels surrounding it. The extraction of values from various environmental variables were extracted using a tool ‘Extracting values to points’ embedded in ArcGIS version 10.5.1. The fire frequency was done using Moderate Resolution Imaging Spectroradiometer (MODIS) burn area index, and the VCI for February 2020 was computed from 2001 to 2020 MODIS normalized difference vegetation index (NDVI) (Table 1). Vegetation condition index is a good indicator of how vegetation is affected by drought (Kogan 1990).

To prepare data for analysis, delineated woodland polygons were used to extract data from a series of environmental data. The first step was to create a regular point shapefile layer with a spacing of 250 m using QGIS version 3.0. The point data were then overlaid on the delineated woodland polygons to create a binary layer (250 m spacing) indicating the presence (ones) and absence (zeros) of woodland polygons. The 250 m spacing was selected because most of the environmental layers had 250 m spatial resolution. The binary layer was used to extract all the environmental variables associated with the presence and the absence data for statistical analysis (Table 1). There were about 1500 presence points and approximately the same amount of absence points for the park. The binary layer with associated environmental variables was converted into a text file for further statistical analysis. For statistical analysis, the binary field in the latter layer is the dependent variable and all the environmental variables are independent.

Descriptive analysis (minimum, maximum, mean, standard deviation and coefficient of variation (standard deviation/mean) × 100) was done to determine the spread and variation of the data sets across the study area, associated with the presence and absence locations. Cross-correlation was also done to determine a level to which explanatory variables are significantly related based on the 95% confidence level.

A generalized linear model (GLM) assuming a binomial distribution (LRM) (Bolker et al. 2009; Hosmer, Lemeshow & Sturdivant 2013), was used to determine optimal environmental variables influencing the spatial distribution of the woodland communities. To minimise multi-collinearity in highly correlated independent variables and over-fitting in training and test datasets, the stepwise logistic model was used: P=11+e−y[Eqn 1]

Where P is the probability of occurrence for the forest patches as explained by several environmental variables, and y = β¹x¹+ β²x²+ β³x³… + c, where β = slope, x = variables and c = intercept: y=β1x1+β2x2+β3x3…+c[Eqn 2]

To determine the optimal number of environmental variables influencing the spatial distribution of woodland communities, several significance tests were done at 95% confidence level (p < 0.05). ‘Optimal environmental variables’ in this article refer to significant variables for explaining the spatial distribution of the woodland communities or a model with the best fit. The optimal variables, and their statistical estimates and confidence levels were also recorded (Table 2). To validate the LRMs used, the data was split into 70% calibration and 30% validation and tested accordingly. The area under-the-curve (AUC) was used as an accuracy indicator for the model (Bradley 1997; Hand 2009). Area under the curve values range from 0 where the model cannot separate the classes (i.e. woodland vs. no woodland) and 1 indicates a good model (two classes clearly distinguished). All the models were run in R statistical programming software using a Rattle Package: AUC=∫01Pr[TP](v)dv[Eqn 3]

Where TP stands for true positive and FP for false positive, and pr[Tp] is a function of v= pr[Fp]

The Spearman rank order correlation coefficient, (ρ) was also used to measure the strength and direction of association that exists between two variables (Eqn 3): ρ=1−6Σdi2n(n−2−1)[Eqn 4]

Where d_ι = difference in paired ranks and n = number of cases.

Figure 2 provides a summary of methods and procedures followed in data collection and analysis.

Thirty-two woodland polygons were identified in GGHNP and delineated using the Google Earth programme. The woodland polygons were converted to a shapefile and edited to ensure topology in ArcGIS version 10.5.1. Field verification was done for grouping woodland polygons into communities. Using GGHNP broad vegetation map and classification by Manfred (1990), woodland polygons were grouped into four different categories (Figure 3), namely Afromontane forests (Olinia, Podocarpus, Kiggelaria), Euclea woodland, Leucosidea woodland, and Protea woodland. The largest extent of woodland community, Leucosidea woodland was associated with the drainage lines and forest margins, although encroaching the footslopes and midslopes in some areas. Afromontane forests (Olinia, Podocarpus, and Kiggelaria) was confined to sheltered gorges and deep valleys. The Euclea woodland community was found at base of cliffs and shelter of large boulders, and also in small clumps in rocky open grassland. The Protea woodland community occurred in open grassland and associated with well drained soils in the footslopes and midslopes.

The descriptive statistics showed relative variability amongst the environmental parameters with appreciably different means (Table 2). For example, variation was relatively higher for TRI (68.78%), topographic roughness (68.03), aspect (60.04%), coarse fragments (37.46%) and TWI (31.33), whereas pH, bulk density, sand and clay content had relatively less variation (2.39%, 3.23%, 7.56% and 8.46%, respectively). However, the topographical position index showed a relatively higher negative value (−818.18%). According to studies on topographic position and landform analysis, negative TPI values signify locations that are lower than their surroundings (gorges, deep valleys), whereas positive TPI values represent locations that are higher than the average of their surroundings (De Reu et al. 2013; Weiss 2001).

The Spearman Rank Correlation showed no significant correlation between most of the environmental parameters (Figure 4). However, significant and positive correlations were found between a few variables such as Topographic Roughness versus TRI (r = 0.985), CEC versus Soil Organic Carbon (r = 0.823); Silt versus CEC (r = 0.783); Elevation versus Soil Organic Carbon (r = 0.757), Clay Content versus Soil Organic Carbon (r = 0.64) etc., whilst negative and significant correlations were found between variables such as Sand versus Silt (r = −0.879); Sand versus Clay Content (r = −0.825); Soil Organic Carbon versus Sand (r = −0.832); CEC versus Sand (r = −0.799).

In determining the optimal number of environmental variables influencing the spatial distribution of woodland communities, topography, edaphic factors (soil chemical and physical properties), fire frequency and remote sensing-based VCI were found to be significant at 95% confidence level (Table 3). The environmental variables such as topographic roughness, TWI, soil organic carbon, CEC and coarse fragments were significant and correlated positively with the woodland assemblages (Table 3). Elevation, soil nitrogen, soil pH, clay content, and fire were also significant but negatively correlated with the woodland assemblages.

The AUC (= 0.81) of the receiver operating characteristic (ROC) was used in evaluating the GLM model for predicting the accuracy of the distributions of the woodland assemblages (Figure 5). The ROC curve was generated by plotting the True Positive Rate (TPR) against the False Positive Rate (FPR). The TPR shows the proportion of positive samples that are correctly predicted, whilst FPR shows the proportion of positive samples that are incorrectly predicted (Wang & Zheng 2013). The AUC indicated that the optimal environmental variables selected as determinants of the spatial distribution of woodland assemblages were based on the good or parsimonious model.