Study area

This study was conducted in a c. 264ha area in the MKFR (Fig. 1), which covers c. 213 083 ha and encircles the 7150-ha Mount Kenya National Park (MKNP), which begins at 3100 m and extends to the highest point, the Batiaan Peak, at 5199 m (Lange & Bussmann Reference Lange and Bussmann1998). The phonolites from main volcanic events c. 2 million years ago form the bedrock of the study area, but the inorganic body of the soils originates from a later coverage of volcanic ashes and pyroclastic rocks (Lange & Bussmann Reference Lange and Bussmann1998). The mountain lies on the equator (latitude 0°10´S, longitude 37°20´E) and forms one of the most pristine mountain ecosystems globally and a remarkable landscapes due to its peaks with rugged glacier-clad summits and diverse forests. The precipitation pattern consists of long rains from March to May and short rains from October to December (Supplementary Fig. S1, available online). The mountain shows a marked vegetational gradient dictated by altitude and rainfall amount. The lower tree line of the forest belt is due to agricultural and pastoral activities. Ocotea usambarensis, which never constitutes pure stands and prefers humid Nitisols and Acrisols, forms the evergreen submontane forests on the southern, south-eastern and eastern slopes of Mount Kenya between 1500 and 2500 m altitude (Lange & Bussmann Reference Lange and Bussmann1998).

Fig. 1. Maps of the study area in the Mount Kenya Forest Reserve showing the locations of selectively logged Ocotea trees (dots).

Acquisition and pre-processing of satellite data

WorldView-3 data were acquired on 15 September 2019 for detecting canopy gaps in selectively logged sites. WorldView-2 imagery acquired on 30 January 2014 and Google Earth were used for historical comparison. In order to cancel out the haze component caused by additive scattering from the RS data, the dark object subtraction method was applied (Chavez Reference Chavez1988). The WorldView-2 satellite captures panchromatic images (450–800 nm) with a spatial resolution of 0.46 m and multispectral images with eight visible–near-infrared (VNIR) bands (400–1040 nm) at 1.84m resolution, while WorldView-3 acquires panchromatic images with a spatial resolution of 0.30 m, multispectral imagery with eight VNIR bands at 1.2 m and eight shortwave-infrared (SWIR) bands (1195–2365 nm) at 3.7m spatial resolution. Additionally, there are 12 clouds, aerosols, vapours, ice and snow (CAVIS) bands with a spatial resolution of 30 m. Only WorldView-3 VNIR bands were used because SWIR bands covering the study area exhibited extensive cloud cover. The satellite data were pansharpened to obtain new bands with a spectral resolution of the multispectral bands and a spatial resolution of the panchromatic band.

Acquisition of field data

Ground truth points were collected in February 2020 using a handheld Global Positioning System (eTrex^® 20 GPS Receiver; Garmin, Olathe, KS, USA) and a pansharpened WorldView-3 image (pixel = 0.30 m). Seventy vegetated gaps formed after illegal logging of Ocotea trees were located in the field (Fig. 1). In the WorldView-3 imagery, the canopy gaps were either partially/fully illuminated or not illuminated at all (Fig. S2b & c). Since the human-made canopy gaps shared similar reflectance characteristics with natural canopy gaps, GPS coordinates of 301 vegetated gaps and 301 shaded gaps were collected, including the 70 canopy gaps formed after illegal logging of Ocotea trees. A vegetated gap is a forest canopy gap with low vegetation inside it, which is the initial stage of vegetation recovery from forest disturbance. Coordinates of the approximate locations of gap centres were recorded and then overlaid on the pansharpened WorldView-3 image. Using a geographic information system (GIS; ArcGIS® v. 10.3; ESRI, Redlands, CA, USA), points were set on the vegetated and shaded canopy gap pixels on the WorldView-3 imagery, and by following the edges of the pixels, the points were made into polygons. Locations of closed forest canopy could be identified using the WorldView-3 imagery; thus, a total of 301 polygons were extracted. Canopy gaps formed after the logging of Ocotea trees had their dimensions collected in the field, such as their dripline measurements, maximum length and compass orientation, together with the maximum width perpendicular to the length. Using the GIS, a map of canopy gaps was generated from ground survey data. The accuracy of canopy gap delineation using RS was assessed by comparing them with the dimensions collected from the field. The ground reference data were then randomly split into train and test datasets of 70% and 30%, respectively.

Feature extraction and selection

In minimizing the effects of data saturation when mapping canopy gaps in dense forests, methods such as vegetation indices (VIs), image transform algorithms, texture measures and spectral mixture analysis have been utilized previously (Malahlela et al. Reference Malahlela, Cho and Mutanga2014, Dalagnol et al. Reference Dalagnol, Phillips, Gloor, Galvao, Wagner and Locks2019). Table 1 presents 55 features (23 means, 23 standard deviations (SDs), 8 ratios and 1 brightness feature) extracted from pansharpened WorldView-3 imagery, and these features are subsequently referred to as variables in the analyses. The VIs were chosen from those sensitive to greenness and plant senescence (Malahlela et al. Reference Malahlela, Cho and Mutanga2014). To accurately extract shaded gaps, the shadow detection index (SDI) of Shahi et al. (Reference Shahi, Shafri and Taherzadeh2014) was used in the modelling. Shade is associated with canopy gaps as nearby trees appear on the edges of gaps (Dalagnol et al. Reference Dalagnol, Phillips, Gloor, Galvao, Wagner and Locks2019). In reducing the redundancy and intercorrelation among the list of potential features, a subset of the best-performing features was extracted from the initial 55 features before classification; for these, approximate optimal thresholds were determined from the reference data using histograms and then adjusted to determine thresholds that resulted in the highest matching accuracy compared to the reference data.

Table 1. List of features used for detecting canopy gaps in the Mount Kenya Forest Reserve: 23 means (of the 15 vegetation indices and 8 visible–near-infrared bands), 23 standard deviations (of the 15 vegetation indices and 8 visible–near-infrared bands), 8 ratios (of the 8 visible–near-infrared bands) and 1 brightness feature (average of the means of bands 1–8).

Spectral separability

Spectral separability measures the distance between two signatures, and the separability between any combinations of variables can be used in the classification. Only the subset from the best-performing features was used in the analysis. The digitized vegetated gaps, shaded gaps and forest canopy samples (Fig. S3) were assigned spectral information using the mean pixel values within their polygons and then used to generate respective signature files. The transformed divergence (TD) index and Jeffries–Matusita (J–M) distance separability measures were used. Divergence (D) is calculated from the mean and variance–covariance matrices of the data representing feature classes (Kavzoglu & Mather Reference Kavzoglu and Mather2000):

(1) \begin{align}{D_{ij}} & = {1 \over 2}tr\left[ {\left( {{\sum\nolimits _i} - {\sum\nolimits _j}} \right)\left( {\sum\nolimits _j^{ - 1} - \sum\nolimits _i^{ - 1}} \right)} \right]\\ &\quad\; + {1 \over 2}tr\left[ {\left( {\sum\nolimits _i^{ - 1} + \sum\nolimits _j^{ - 1}} \right)\left( {{\mu _i} - {\mu _j}} \right){{\left( {{\mu _i} - {\mu _j}} \right)}^T}} \right]\end{align}

The TD is introduced to reduce the impact of well-separated classes that may raise the average divergence value and make the divergence measure misleading (Kavzoglu & Mather Reference Kavzoglu and Mather2000):

(2) $$T{D_{ij}} = c\left[ {1 - {e^{{{ - {D_{ij}}} \over 8}}}} \right],$$

where tr[·] is the trace of a matrix, which is the sum total of the diagonal elements of the matrix, and Σ _i and Σ _j are the variance–covariance matrices of classes i and j; μ _i and μ _j are the corresponding mean vectors; c is a constant value defining the range of TD values.

The J–M distance between distributions of two classes ω _i and ω _j has been defined as follows (Richards & Jia Reference Richards and Jia1999):

(3) $$J\!\!-\!\!{M_{ij}} = \;2\left( {1 - {e^{ - {B_{ij}}}}} \right),$$

where B _ij is the Bhattacharyya distance (Kailath Reference Kailath1967) computed as:

(4) \begin{align}{B_{ij}} &= {1 \over 8}{({\mu _i} - {\mu _j})^T}{\left( {{{{\sum _i} + {\sum _j}} \over 2}} \right)^{ - 1}}({\mu _i} - {\mu _j})\\ &\quad + {1 \over 2}ln\left( {{1 \over 2}{{|{\sum _i} + {\sum _j}|} \over {\sqrt {|{\sum _{_i}}||{\sum _j}|} }}} \right),\end{align}

where μ _i and μ _j are the mean reflectances of species i and j, Σ _i and Σ _j correspond to their covariance matrices, with |Σ _i | and |Σ _j | being the determinants of Σ _i and Σ _j , respectively, $ln$ is the natural logarithm function and T is the transposition function. The J–M distance improves the Bhattacharya distance by normalizing it to between 0 and 2.

Pairwise feature comparison

A correlation matrix was computed from the reference samples (Table 2). Similarity scores were calculated between each pair of variables to determine whether or not two variables were co-referent. Pairwise comparisons are in form of a matrix: C = [c _kp ] _{n × n} , where c _kp is the pairwise comparison rating for kth and pth criteria. The matrix C is reciprocal; that is, c _pk = $c_{kp}^{ - 1}$ , and all of its diagonal elements are unity; that is, c _kp = 1, for k = p (Malczewski Reference Malczewski, Huang and Cova2016).

Table 2. Pairwise correlations between the best-performing variables extracted from the WorldView-3 visible–near-infrared bands, computed from the reference samples. (See Table 1 for acronym definitions.)

RF and SVM models

The results in this study were obtained by training the RF and SVM models in R software (R Core Team, Vienna, Austria). In RF, to differentiate between predefined categories, decision trees recursively partition the source set into subsets with bagged samples by univariate splits at internal nodes (Breiman Reference Breiman2001). Before running the model, the number of decision trees (ntree) and the number of predictor variables (mtry) randomly selected at each node are defined. The RF model aggregates predictions from all decision trees, then the majority vote of all trees assigns a final class for unknown features (Breiman Reference Breiman2001). Using the grid search method, the mean decrease in accuracy (MDA) was used to extract a subset of the best-performing variables (Breiman Reference Breiman2001). The MDA shows how much accuracy the model losses by excluding each variable; therefore, the higher the MDA value, the more important the variable in the model.

The SVM assigns a class from one of the two possible labels when test data are introduced after the training phase (Vapnik Reference Vapnik2000). The SVM separates the original data while maximizing the margin between classes and minimizing the misclassification error (Vapnik Reference Vapnik2000). An advantage of ML classifiers such as SVM is that they are suited for extreme case binary classification. For any two distinguishable classes with k samples represented by (x ₁, y ₁), …, (x _k , y _k ), where x ∈ R ⁿ is an n-dimensional space and y is a class label with values of +1 or −1, SVM will look for an optimal hyperplane defined by w = (w ₁, …, w _n ) and b, such that (Huang et al. Reference Huang, Song, Kim, Townshend, Davis and Masek2008):

(5) $${y_i} = \left[ {w\cdot\;{x_i} + b} \right]1i\; = \;1,\; \ldots \;,k$$

The hyperplane can be located by minimizing the norm of w or the following function under the above inequality constraint (Huang et al. Reference Huang, Song, Kim, Townshend, Davis and Masek2008):

(6) $$F\left( w \right) = 0.5\left( {w\; \cdot \;w} \right)$$

Using kernel functions, SVM applies non-linear decision boundaries and introduces a cost parameter C and gamma parameter γ to quantify the penalty of misclassification errors and to give the curvature weight of the decision boundary, respectively. The robust radial basis function was selected as it has fewer parameter values to predefine. A parameter search must be done to select the best C and γ for a certain classification problem. Therefore, the γ parameter needs to be predefined (Huang et al. Reference Huang, Song, Kim, Townshend, Davis and Masek2008):

(7) $$K\left( {_i,\;{x_j}} \right) = {e^{ - \gamma {{\left( {{x_i} - \;{x_j}} \right)}^2}}}$$

The cost parameter C also needs to be predetermined for the canopy gap-mapping problem. A cross-validation quantitative analysis of pairs of values for the C and γ parameters was carried out. The combination of parameters with the lowest error was chosen to train the algorithm. Both RF and SVM classifiers were trained using 70% of the ground reference data, and for robust classification results the ten-fold cross-validation method was repeated ten times.

Measures of model performance

The performance of the RF and SVM models was evaluated using 30% of the ground truth data. Confusion matrices with overall accuracy, kappa coefficient and producer’s and user’s accuracies were computed and averaged over ten repetitions. Overall accuracy is computed by summing the number of pixels correctly classified divided by the total number of pixels, while producer’s accuracy is the percentage of particular classes on the ground that are indicated as such on the classified map (Mutanga et al. Reference Mutanga, Adam and Cho2012). The user’s accuracy shows the probability that a pixel indicated as a specific feature is classified as such on the classification map (Mutanga et al. Reference Mutanga, Adam and Cho2012). The kappa coefficient is the difference between the observed accuracy and the agreement that would have been expected by chance. The McNemar test was used to compare and indicate the statistical significance of any differences between the two classifiers (Foody & Mathur Reference Foody and Mathur2004).