Mapping invasive strawberry guava (Psidium cattleianum) in tropical forests of Mauritius with Sentinel-2 and machine learning

ABSTRACT Invasive flora alien species (IAS) pose a serious threat to biodiversity in tropical forests. This study evaluated different classification approaches for remote sensing of the target IAS, strawberry guava (Psidium cattleianum), using open-access satellite imagery to aid research and understanding of IAS impacts in tropical forests. This is likely the first study using freely available open-access imagery to map strawberry guava, and the first mapping of invasive species in Mauritius using remote sensing. Ground-based observations in Black River Gorges National Park (BRGNP), Mauritius, resulted in 4670 ‘ground truth’ samples across three land cover classes representing different levels of strawberry guava canopy cover; 70% of these ground truth data were used for training and model tuning, and 30% were held-out for testing. Classification was performed using Sentinel-2 MSI images and Google Earth Engine. Computation of Jeffries-Matusita distance supported the addition of Grey Level Co-Occurrence Matrix texture measures, and spectral indices (ARI1, and ReNDVI) for classification, since these features increased separability of strawberry guava cover classes substantially over using spectral bands alone. Higher than 80% overall and individual class accuracies were achieved with both Random Forest (RF) and Support Vector Machine (SVM) classification of strawberry guava in BRGNP. RF is recommended for similar future applications for its higher overall accuracy (97.60% ± .20% at 95% confidence), lower variation, and prediction of more homogenous class regions. This study found that strawberry guava canopy cover in BRGNP was highest in central regions and correlated to steeper slopes, with little overall change from 2016 to 2020.


Introduction
Tropical islands have some of the most productive and biologically diverse habitats on the planet and many have thus been termed biological hotspots (Myers et al. 2000), yet they are facing significant threats including habitat loss and invasive species. Control of invasive flora alien species (IAS) is important to maintain global biodiversity, especially in these tropical regions (IUCN -SSC 2017). One such IAS is strawberry guava (Psidium cattleianum). Although native to humid and wet environments, strawberry guava is both drought and frost tolerant, so it can be found in diverse environments from grasslands to tropical montane forests (CABI 2019). With such global invasion potential to sensitive tropical ecosystems, the importance of establishing accurate mapping of this IAS is clear, yet at the time of this study, no published research was found using open-access imagery to map strawberry guava invasion.
Satellite remote sensing represents a valuable tool for studying areas of IAS (Cuneo, Jacobson and Leishman 2009;Hoyos et al. 2010) to direct resources and understand ecosystem interactions. Vegetation absorbs photosynthetically active radiation in the visible region of the spectrum (VIS), from 400 nm to 700 nm, and reflects light outside this region, especially the near infrared (NIR) (McCoy 2004). The steep change from VIS absorption to NIR has been well studied for its significance in discriminating vegetation classes (i.e. Cochrane 2000;Schuster, Förster and Kleinschmit 2012). Yet, approaches for detecting IAS using remote sensing are often species, site, and sensor specific with research often utilizing expensive high-resolution imagery (Huang and Asner 2009). While this method can produce high classification accuracy on a small scale, it is costly and therefore inaccessible to many organizations (Fisher et al. 2018). Advances in openaccess satellite imagery, namely that of Copernicus Sentinel-2 MultiSpectral Imager (MSI) (ESA 2018), provide a possible means for more broadly applicable methods (i.e. Laurin et al. 2016;Puletti, Chianucci and Castaldi 2017). Vegetation classification to the species level often requires high spatial resolution (Fassnacht et al. 2016;i.e. Immitzer, Atzberger and Koukal 2012) and high spectral resolution imagery (Feret and Asner 2013;Barbosa et al. 2016;Ferreira et al. 2019). However, studies also show that the classification accuracy of vegetation to the species level is further dependent on the classification method utilized (Fassnacht et al. 2016;Laurin et al. 2016). Historically, simple classification algorithms like the Maximum Likelihood Classifier (MLC) have dominated land classification research (Yu et al. 2014) due to its ease of use and relatively low computing power required (Huang, Davis and Townshend 2002). With significant increases in computing power, theoretically superior machine learning algorithms, like support vector machine (SVM) and random forest (RF), have become computationally accessible for the average user and thus recommended for higher accuracy image classification (Yu et al. 2014;Maxwell, Warner and Fang 2018). Unlike MLC, these machinelearning algorithms do not make assumptions such as normal (Gaussian) data distribution, but instead algorithm parameters are tuned to the training data (Jensen et al. 2011). Furthermore, empirical studies such as Foody and Mathur (2004) and Huang, Davis and Townshend (2002), as well as a review by Yu et al. (2014), showed that machine learning algorithms were more accurate than the traditionally used MLC.
The aim of this study was to evaluate different approaches for classifying IAS strawberry guava using open-access satellite remote sensing in order to aid research and management of IAS impacts in tropical forests. This included investigating the importance of additional spectral features for separability of IAS cover classes and comparing machine learning classification algorithms for classification accuracy. We focused on the biologically diverse tropical forests of Black River Gorges National Park (BRGNP) on the island of Mauritius, which are heavily invaded by the IAS, strawberry guava (Psidium cattleianum), and where there has been no published research to date using satellite remote sensing to map invasive species. Using Sentinel-2 MSI imagery and in-situ field data of strawberry guava in the canopy of BRGNP we sought to 1) evaluate the importance of additional vegetation indices and texture features (e.g. gray-level co-occurrence matrix, GLCM) in increasing separability of classes of strawberry guava invasion with Jeffries-Matusita distance (JM); 2) assess the accuracy of support vector machine (SVM) and random forest (RF) supervised classifications of imagery; and 3) use classified images to gain an understanding of the distribution of the target IAS in BRGNP.
Image selection, processing, and classification started with ground truth data acquisition and partitioning, and then followed sequentially by image processing. Classification optimization, classification, accuracy assessment, and land cover change detection were completed with Google Earth Engine (GEE) (Gorelick et al. 2017). Code adapted from Google Developers (2020) was used for most data analysis steps except feature selection where RStudio with R programming language was used (R Core Team 2017).

Study Location
The island of Mauritius (1,865 km 2 ), located in the southwest Indian Ocean, is part of the small island developing state the Republic of Mauritius (Figure 1(a,b)). Since the establishment of the first colony and fortified settlement in Mauritius by the Dutch in 1638, there has been rapid loss of native forests (Norder et al. 2017) and introduction of invasive alien species (Florens et al. 2017;Vaughan and Wiehe 1941). A 2010 review by Baider et al. determined 691 species of angiosperms were native to Mauritius, including 273 endemic plants to the island, but further noted <5% of intact habitat was left. Consequently, of the island's endemic species, 30 were extinct and 81.7% were threatened (Baider et al. 2010). The review concluded that IAS posed the greatest threat to native species, which has been further substantiated by studies quantifying density, biomass, and basal area of invasive versus native species over time (Virah-Sawmy et al. 2009;Florens et al. 2017).
The island has a tropical maritime climate with a warm, wet summer from November to April and a drier, cooler winter from June to October with ~70% annual rainfall occurring in summer when tropical cyclones can affect the island (MEPU 2010;MMS 2020a). This study focused on the Black River Gorges National Park (BRGNP) located in southwest Mauritius (Figure 1(c)).
BRGNP was formally established in 1994 and is the largest protected area on the island at approximately 66 km 2 (NPCS 2017). Elevation ranges from about 93 m in the northeast valley to 730 m on the peaks of the Central Plateau. Land cover is predominantly forest with upland rainforest (>400 m elevation) on the eastern side and drier lowland forest (<400 m elevation) on the western side of the park (Nigel, Rughooputh and Boojhawon 2014;NPCS 2017). The park contains the largest area of native and endemic species on the island (Vaughan and Wiehe 1937;NPCS 2017) but has also been highly invaded by strawberry guava (Vaughan and Wiehe 1941;Virah-Sawmy et al. 2009;Florens et al. 2016).

Strawberry guava in Mauritius
Strawberry guava (Psidium cattleianum) is an evergreen tree native to Brazil that has been designated an invasive species in 23 tropical and subtropical regions, including Mauritius (CABI 2019). Both in Mauritius and globally, strawberry guava tends to grow in dense, monospecific stands that strangle out native vegetation and pose serious threats to native plants (Macdonald et al. 1991;Patel 2012;Tng et al. 2016;Florens et al. 2017). Strawberry guava can establish in shady understories and reach 8 m tall, and once mature overtake the canopy (CABI 2019; Barbosa et al. 2016). Alternatively, it can establish in the full sun as a shrub growing 2-4 m tall (CABI 2019). The fruits are eaten by humans (Patel 2012) and it is likely for this reason that the plant was introduced to Mauritius around the 1830s (Lorence and Sussman 1986;Kueffer and Mauremootoo 2004). Strawberry guava spreads prolifically from seeds and can also re-sprout from stumps, making removal and eradication complicated (Vaughan and Wiehe 1937;Fowler et al. 2000;NPCS 2017). Recent mapping of invasive species and estimates of spread in Mauritius have been conducted by field surveys and expert knowledge. At the time of this study, no published research was found on remote sensing of IAS in Mauritius, or on using open-access satellite imagery to map strawberry guava. Given the fact that strawberry guava tends to grow in dense monospecific stands with widespread range, it is an excellent candidate for using satellite mapping to track its distribution on the island (Patel 2012).

Field surveys
Field surveys were conducted to gather training and testing data for classification. Many forests of BRGNP have multiple strata, but since satellite imagery only represents layers exposed to the sun, strawberry guava cover was estimated only for the sun-exposed canopy as suggested in McCoy (2004). For simplicity, only strawberry guava cover was estimated in each plot, and thus 'cover' refers to only strawberry guava canopy cover and the remaining area was considered non-invaded.
Given BRGNP has steep cliffs and dense vegetation, trails were used as survey transects to ensure accessibility. Transects were selected to sample a range of environmental conditions, invasion levels, and ecosystems in BRGNP to maximize representation of the spectral profile for each class (Hoyos et al. 2010;Campbell and Wynne 2011). To determine which transect locations were the most appropriate, the most recent cloud-free Sentinel-2 image of BRGNP (Level 1C, acquired 18 July 2017) was examined for metrics known to indicate vegetation spectral differences. Trail shapefiles were extracted from OpenStreetMaps (OSM; Geofabrik 2018).
Since the distribution of strawberry guava was not known beforehand, plots were opportunistically sampled along transects to gather a sufficient sample of all levels of strawberry guava invasion. Plots were chosen to be as homogenous as possible, avoiding overlaps in classes (Campbell and Wynne 2011). It is recommended that plots contain multiple pixels (Campbell and Wynne 2011), so areas of >20 m × 40 m were used to encompass multiple pixels from both 10 m and 20 m Sentinel-2 bands. Plots were recorded as polygons using a handheld GPS (Trimble Juno 3B, Trimble, USA) with nodes recorded to within 5 m accuracy.
A visual estimate of strawberry guava crown cover in the canopy layer was recorded for each plot using an approach similar to Barbosa et al. (2016). To minimize bias, two field workers independently estimated the proportion of guava in the canopy layer compared to the total area of the plot and then agreed on a final guava canopy cover proportion. From the agreed guava canopy cover proportion, plots were then assigned a cover scale value based on the ordinal classification of the Braun-Blanquet cover scale (Braun-Blanquet 1932; Figure 2). The scale values ranged from 1 (least cover, covering less than 20% of the area) to 5 (most cover, covering more than 75% of the area) (Braun-Blanquet 1932). In this way, if the visual proportion of guava cover in the canopy was estimated by observers to be 10%, then the plot would be assigned a cover scale value of 1. This method sacrifices detail to gain information on a larger quantity and area of plots and this approach has been successfully used in previous vegetation remote sensing applications (Zak and Cabido 2002;Damgaard 2014;Barbosa et al. 2016).

Land cover classes
For classification, field data was separated into three strawberry guava land cover classes to represent 'none to very low' (Braun-Blanquet value 1 and 2), 'mixed' (Braun-Blanquet value 3 and 4), and 'high' plant canopy cover (Braun-Blanquet value 5; Figure 2). These categories are based on Hoyos et al. (2010). Additionally, Congalton (1991) estimated at least 50 samples per class should be used for accuracy assessment and the larger the land area and spectral variation in each class, the more samples are desired. Accordingly, sites totalled approximately .024 km 2 for 'none to very low', .086 km 2 for 'mixed' and .143 km 2 for 'high' guava. After resampling at 10 m pixels this resulted in 2207 samples for 'none to very low', 924 for 'mixed', and 1539 for 'high' guava cover totalling 4670 ground truth samples.

Data partitioning
Ground truth data were partitioned in GEE using stratified random sampling for each class to separate into 70% training data and 30% held-out for testing as is common (i.e. Laurin et al. 2016). For parameterization of SVM and RF classifiers, five-fold cross-validation was performed using the training data. Accordingly, the training data was split further into five groups of 80% training and 20% validation, as is common (Foody and Mathur 2004). All training data was thus used in model tuning and parameter values could then be selected to minimize potential overfitting. Performance was evaluated as the mean of overall accuracies of the five-fold cross-validations.

Transect selection
Satellite images from Sentinel-2 were acquired from USGS EarthExplorer (https:// earthexplorer.usgs.gov) and the analysis for transect selection was conducted in ESA SNAP 6.0 S2TBX (Sentinel-2 Toolbox). The image was subset to bands B2 -B8A and B11 -B12, and 20 m resolution bands were resampled to 10 m (Table 1). Metrics analysed using the inbuilt functions in SNAP were leaf area index, normalized difference vegetation index (NDVI), and unsupervised classification using Kmeans with 100 iterations (Baldeck and Asner 2013) and an arbitrary 10 clusters. Since vegetation is known to change with elevation in BRGNP (Vaughan and Wiehe 1937) a 30 m resolution digital elevation map (DEM) from the Shuttle Radar Topography Mission (SRTM NASA V3) was also utilized. Prioritized trails were those passing through and including a range of index values, K-means classes, and elevation ( Figure 1).

Pre-Processing
Image retrieval and analysis was conducted in GEE since it is open-source and has inbuilt image databases (Gorelick et al. 2017). Sentinel-2 Level 1C imagery was used since it is pre-processed to orthorectified top-of-atmosphere reflectance (ESA 2017).

Image selection and masking
To overcome the issue of cloud-cover, Sentinel-2 images from before and after the survey dates were aggregated into a single cloud-free composite of median values. Multiple time window sizes inclusive of field survey dates were examined. A phenology-based cloud mask adapted from Simonetti et al. (2015) by Ceccherini, Verhegghen and Dario (2017) was applied to the image collection to create a composite cloud-free image. The mask algorithm was developed to produce high accuracy results even in tropical regions where cloud-free observations may be limited (Simonetti et al. 2015). The algorithm classifies clouds and cloud shadows using fuzzy rules based on cloud and shadow spectral signatures in the blue, green, red, NIR, and SWIR bands, as well as NDSI (Normalized Difference Snow Index) and NDVI (Simonetti et al. 2015). In each image of the collection, pixels classified as clouds or shadows are masked out, or set to 'no data', and then the median value of each pixel across the collection is used to produce a single composite image (Ceccherini, Verhegghen and Dario 2017). User-defined model parameters include specifying a maximum cloud cover to filter the image collection by, as defined in Sentinel-2 image metadata. For this study, maximum cloud cover was set to 100% to include all images in the date range, as recommended in the model (Ceccherini, Verhegghen and Dario 2017).
Smaller time windows produced increasingly patchy composites with some cloud remnants due to lack of cloud-free areas. The final 2018 composite used, ±4 weeks of field survey dates, was chosen for its proximity to survey dates combined with favourable visual inspection (minimum patchiness and adequate cloud masking) and was created from 15 scenes taken from 30 March 2018-18 June 2018 (Appendix: Table A1). The 2019 and 2020 images used the same eight-week timeframe as the 2018 image (30 March−18 June). Sentinel-2B only began imaging Mauritius in July 2017, so frequency of images for 2017 and 2016 composites were much lower. The image collection period was therefore expanded for the 2017 and 2016 images based on when cloud-minimal imagery was available to produce satisfactory composites (Appendix: Table A1). The 2017 composite used 12 scenes taken between 30 March 2017-18 July 2017, and the 2016 composite used 9 scenes taken between 14 April 2016-13 July 2016. Only the 2018 composite image was used in feature selection algorithms, classification algorithm parameter tuning, and classifier training.
Since only classes of vegetation were desired, Mare Longue Reservoir in the northeast of BRGNP was masked. The version 1.7 land cover maps of Hansen et al. (2013), available in the GEE data catalogue, were used to set all pixels in the 'water' class to 'no data'. Water remnants remained at the lake edge, so a threshold approach was applied to mask all remaining pixels with reflectance <4.5% in band 6. Visual inspection verified sufficient masking of water pixels.

Subsetting and resampling
The images were then subset to the 10 bands of interest, B2 -B8A and B11 -B12, and the 20 m bands were resampled to 10 m using nearest neighbour resampling as recommended by Richards and Jia (2006) so that all bands had the same resolution.

Vegetation indices features
We computed four vegetation spectral indices: NDVI, ReNDVI, SRI, and ARI1 as in Laurin et al. (2016) and Puletti, Chianucci and Castaldi (2017). NDVI, one of the most widely used vegetation indices (Khatami, Mountrakis and Stehman 2016), computes the normalized difference between near infrared and red bands to indicate the density of green leaves (Sellers 1985): where ρ 842 is reflectance in the near-infrared band, Sentinel-2 B8, and ρ 665 is reflectance in the red band, B4. ReNDVI minimizes the differences in leaf surface reflectance in order to be more sensitive to changes in chlorophyll content (Gitelson and Merzlyak 1994;Sims and Gamon 2002). ReNDVI was calculated as: where ρ 740 is maximum reflectance in the red-edge, B6, and ρ 705 is minimum reflectance in the red-edge, B5. Additionally, SRI calculates a non-normalized estimation of chlorophyll content (Sellers 1985): ARI1 was also calculated to estimate anthocyanin content (Gitelson, Merzlyak and Chivkunova 2007): where ρ 560 is reflectance in the green band, B3.

GLCM texture features
Texture features were computed using the GEE GLCM function. The operation requires defining a kernel and size of neighbourhood, also known as a sliding window, over which the functions are applied. The larger the window, the larger the area each class must occupy in the scene to discriminate class patterns (Conners, Trivedi and Harlow 1984).
Since class differentiation in this study depended on finer detail, a smaller window size of 3 × 3 was used (Conners, Trivedi and Harlow 1984;Laurin et al. 2016). The GEE GLCM function computes 18 texture measures for each input band, described in Haralick, Shanmugam and Dinstein (1973) and Conners, Trivedi and Harlow (1984), creating an output layer of each texture feature, thus totaling 180 texture features for each image. For simplicity, feature bands were subset to include only texture measures supported in relevant literature to increase spectral separability of vegetation classes. These consisted of sum average, variance, contrast, dissimilarity, entropy, angular second moment (ASM), and correlation (Laurin et al. 2016;Hall-Beyer 2017). It should be noted, Hall-Beyer (2017) used the directional GLCM mean, while this study used the alternative, non-directional GLCM sum average, which is computed in GEE GLCM features.

Feature selection
Each composite image with indices and texture features totaled 84 bands, or features. To increase efficiency in both computing power and time, and well as decrease possible noise from extraneous features, spectral separability analysis was carried out on the 2018 composite with training data to select a subset of the most useful features (Swain and King 1973;Mausel, Kramber and Lee 1990;Dutra and Huber 1999).
The widely used Jeffries-Matusita distance (JM) (Swain and King 1973;Dutra and Huber 1999;Richards and Jia 2006) was computed to determine the statistical spectral separability of classes using iterative feature subsets. This method is independent of the classifier and assumes the greater the spectral separability between classes, the lower the probability of classification error (Swain and King 1973). Therefore, the feature subset of n features with the highest JM is selected. Feature selection was performed using the exported training samples and composite image in RStudio Version 1.1.423 with R programming language (R Core Team 2017), 'raster' package (Hijmans 2017), and 'varSel' package (Dalponte and Ørka 2016). Since this study evaluated multi-class distance, average JM between class pairs in each subset was used, as in Dutra and Huber (1999). To determine which features to include in each subset, an automatic Sequential Forward Floating Selection (SFFS) method was used as this produces reliable selection with fast computation (Pudil, Novovičová and Kittler 1994;Kavzoglu and Mather 2000;Dalponte et al. 2013).
The varSel JM SFFS algorithm limits feature inputs to 65, so a semi-automatic method was first used to reduce the number of features from 84 to 63. To do this, separability was evaluated using all 10 spectral bands, all 4 vegetation indices, and one texture measure stack at a time (e.g. entropy of each of the 10 spectral bands). The 7 highest ranking bands for the texture measure were then selected, based on bands added earlier to the subset ranking more important. This was re-done for each texture measure. A JM was then computed using the 10 spectral bands, 4 vegetation indices, and each of the 7 texture measures in their 7 highest ranking bands (49 total texture features). A final subset was chosen from this result.

SVM classification and parameterization
To optimize the SVM classifier, the influence of parameter values of cost (C, allowed misclassifications) and gamma (γ, distance of influence of training points) on model performance was explored by comparing five-fold cross-validation performance for varying C and γ values. A C and γ parameter value was then selected which optimized model performance and minimized potential overfitting.
The SVM classifier maps training data, or the input vector, in a multidimensional feature space using a kernel function, and then uses optimization algorithms to maximize the distance between classes (Huang, Davis and Townshend 2002). The optimization algorithms include introducing slack variables to relax constraints and allow misclassified pixels according to a user-defined cost (C) (Cortes and Vapnik 1995). There are multiple methods to convert SVM from binary to multi-class, but in this study the GEE default method of voting was used as in Gualtieri and Cromp (1998) and Hoyos et al. (2010). For kernel mapping, radial basis function (RBF) was implemented (Foody and Mathur 2004;Hoyos et al. 2010;Nemmour and Chibani 2010). Huang, Davis and Townshend (2002) showed that RBF provided sufficient accuracy and the parameter value (γ), which determines the distance of the influence of a training point, was less sensitive to change than that of alternatives.
In this study, the γ parameter was determined by evaluating performance resulting from repeating five-fold cross-validation, with γ varying from .005 to 1.0 as in Foody and Mathur (2004). The SVM cost parameter (C) was similarly determined by varying C from 1 to 10, common values found in the literature (Foody and Mathur 2004;Nemmour and Chibani 2010;Laurin et al. 2016).
Variation of C > 3 had little effect on performance, as expected (Cherkassky and Ma 2004), with a variation of only 1.6% (Appendix: Fig. B1(a)) so a value of 3 was selected. Performance was more sensitive to changes in γ (Appendix: Fig. B1(b)) with a decrease from 77.4% accuracy at .005 to 50.1% accuracy at 1.0, indicating overfitting to training data. A γ of .020 was selected.

RF classification and parameterization
Optimization of the RF classifier was conducted similar to SVM. The parameters n tree (number of trees) and m try (number of features) were varied and the influence on model performance was assessed by evaluating the five-fold cross-validation performance across parameter values.
Using RF (Breiman 2001) a decision tree was created for each bootstrap to determine what decisions, or splits, should occur to assign each pixel to a class. At each node of the tree a random subset of features was selected from which to determine the split (Breiman 2001). The resulting classifier was applied to the selected images with the class with the highest number of predictions, or votes, assigned to pixels (Breiman 2001).
RF has the advantage of having only 2 input parameters, the number of trees (n tree ) and number of features selected for each split (m try ), which are easily estimated without overfitting (Liaw and Wiener 2002). For n tree , Breiman (2001;2003) notes that accuracy and stabilization increase as n tree increases and suggests starting at 1000 and setting the value generously high. However, a review of RF application to remote sensing by Belgiu and Drăgut (2016) noted most studies use 500 trees. Thus, performance resulting from the 5-fold cross-validation was assessed for values of n tree from 100 to 1000. To select m try , Breiman (2003) suggests using the default, the square root of the number of features, but also advises testing half and double the default. Accordingly, performance for m try was assessed for values between 2 to 8.
Varying n tree had little effect on performance, with total variation in performance from examined values of only .6% accuracy (Appendix: Fig. B2(a)). However, higher n tree resulted in significantly longer computing time so the commonly used value of 500 was selected. The number of features selected per split also had little effect on performance with a variation of only .4%, so the default and generally accepted setting of the square root of the number of features was selected (m try = 4) (Appendix: Fig. B2(b)).

Accuracy assessment
Classification accuracy of the 2018 image was assessed by calculating overall accuracy, a confusion matrix, and kappa coefficient (K) for each classification method using the heldout ground truth testing dataset. Training and testing data was then resampled 10 times, using a 70% training and 30% testing split, with classification accuracy and K calculated each time to obtain overall accuracy mean and confidence interval.
The confusion matrix is a common measure of land classification accuracy and is considered appropriate when presented alongside a quantitative metric with confidence limits (Foody 2002;Congalton and Green 2009). From the matrix, one can calculate the producer's accuracies, which show how often reference data was classified correctly on the map, and the user's accuracies, which show how often the classified area is present on the ground. The kappa coefficient, first proposed by Cohen (1960), measures the difference between the observed probability of agreement between reference and classified data, and the probability of agreement that could be expected by chance (Campbell and Wynne 2011). K ranges from 0 to 1 with values closer to 1 indicating more perfect agreement.

Land cover change detection
Change in land cover classes was assessed by post-classification comparison (Singh 1989;Warner, Almutairi and Lee 2011). Since classes were scaled sequentially from 1 (low/no strawberry guava) to 3 (high strawberry guava), image differencing could show changes in pairs of images (Warner, Almutairi and Lee 2011). To do this, the classified earlier year image was subtracted from the classified later year image for each image pair of each classifier. With this method, areas that stayed the same class had values of 0, areas that decreased in strawberry guava class had negative values (−1 or −2), and areas the increased in strawberry guava class had positive values (1 or 2). For example, if a region was class 3 (high strawberry guava) in 2017 but class 1 (low/no strawberry guava) in 2018, it would have a change of −2 for shifting down 2 strawberry guava cover classes. Land cover areas were assessed calculating the total area of pixels in each class for each year. It is important to stress, the results of this method are not linear or continuous, but rather highlight changes between the discrete classes.
Code for the GEE processing and classification is available at https://github.com/ laurmk/PCattleianumInMauritius. Version 1.1, used in this paper, is archived at 10.5281/ zenodo.5812245.

Field data observations
Conservation management areas (CMAs), fenced areas heavily managed for invasive species control, were encountered along trails including the Macchabee and Mare Longue (Figure 3(a,b)). CMAs were plotted as 1 on the Braun-Blanquet scale, or the low/no strawberry guava class (class 1) for classification, after visual confirmation of guava absence. Areas with large dense monospecific strawberry guava stands (Figure 3 (c,d)) were plotted as 5 on the Braun-Blanquet scale, or the high strawberry guava class (class 3) for classification. In these areas, it was noted that plant stem density was so high the forest was near impenetrable. Where strawberry guava was growing in shade or part-shade, it tended to grow taller and less dense with reduced branching (Figure 3 (e)). These areas were often considered mixed classes, or 3 to 4 on the Braun-Blanquet scale, when the canopy of taller trees was more open. Many areas of strawberry guava, particularly at higher elevations in the central and eastern central region, were fruiting (Figure 3(f)).

Feature selection
Using SFFS and JM, the 7 highest ranking bands (rank of 1 being highest importance) for each texture measure were selected for an overall JM analysis (Figure 4). B3 and B4 consistently ranked higher, while B8 ranked in the top 7 for only 2 measures, making it the least important for computing useful texture measures.
As expected, the increase in separability gained from adding an additional feature to the subset decreased as the subset got larger, indicating model saturation was being approached. The 10 spectral bands, 4 vegetation indices, and 49 texture features produced a JM of approximately 1.41 using the largest subset of 62 of 63 features ( Figure 5). To maintain computational efficiency and minimize dimensionality, the subset of 20 features was selected which had a moderate JM separability of approximately 1.35 (Appendix: Fig. C1). For comparison, the JM was .74 for the 9-feature subset using the 10 spectral bands, and 1.11 for the 13-feature subset using the 10 spectral bands and 4 vegetation indices.
The selected subset included 20 of the top 21 most frequently selected features, with sum average in band B12 excluded ( Figure 6). The subset consisted of 2 vegetation indices, 2 contrast measures, 5 variance measures, 6 sum average measures, 3 dissimilarity measures, 1 ASM measure, and 1 entropy measure (Appendix: Fig. D1). This subset did not include any of the 10 spectral bands, as these were only included in subsets of > 24 features.

Land cover classification and change detection
All classified images using both classifiers showed similar spatial distribution of canopy cover (Figure 7). Areas of high strawberry guava were concentrated in the central and north-western areas of the park, while areas of low/no guava were most concentrated in the southwestern area of the park. SVM did not classify many areas as mixed, but RF concentrated substantial areas of mixed in the eastern side of the park. SVM tended to classify more small, isolated areas, creating a more speckled image. Conversely, classes in RF images tended to be more homogenous areas.

Strawberry guava and slope
A visual association between high plant cover and steep slopes was observed, so the distribution of slopes in pixels classified as high strawberry guava was compared to areas classified as low/no strawberry guava and mixed in the 2018 random forest classified image. A Mann-Whitney U test showed that there was a significant difference at 95% confidence between slopes in areas classified as high compared to low/no strawberry guava (p < 2.2 e −16 ), and areas classified as high strawberry guava compared to mixed (p <   2.2 e −16 ). Mean slope in areas classified as high strawberry guava was 18.60° with a standard deviation of ±11.06°, compared to low/no strawberry guava mean slope of 12.14° with a standard deviation of ±8.18° and mixed mean slope of 12.50° with a standard deviation of ±8.93° (Figure 8).

Land cover change
Total area within the BRGNP shapefile boundaries was 66.97 km 2 , of which 65.92 km 2 was classified in the 2016 image, 65.91 km 2 in 2017 and 2018, and 65.94 km 2 in 2019 due to fewer water pixels masked. In 2020, 65.76 km 2 were classified due to an additional small masked area in the composite from persistent cloud cover in the image collection.
The low/no strawberry guava class showed small variation in predicted total area between classifiers and between most years with a 5-year mean of 25.62 km 2 and standard deviation (s) of 2.92 km 2 for RF and 22.68 km 2 and 1.49 km 2 , respectively, for SVM (Figure 9). High strawberry guava showed more yearly variation in predicted total area with RF, increasing 3.24 km 2 from 2016 to 2018, and then decreasing 7.42 km 2 from 2018 to 2019 which is mostly maintained into 2020. There was larger variation in high strawberry guava between classifiers with a 5-year mean of 26.72 km 2 (s = 3.71 km 2 ) for RF and 37.46 km 2 (s = 1.76 km 2 ) for SVM. The mixed class accounted for the smallest area in all images, with relatively large variation between classifiers. SVM predicted small yearly variation of total mixed area with a mean of 5.75 km 2 and standard deviation of .40 km 2 . Conversely, RF predicted a 5-year mean of 13.54 km 2 and standard deviation of 2.63 km 2 . The deviation in the RF mixed mean comes from a 38% decline from 2016 to 2018, and then a 54% increase from 2018 to 2020. Both classifiers showed a relatively large decrease in high and increase in low/no strawberry guava from 2018 to 2019.
When examined spatially, class changes from year to year varied more readily ( Figure  10). From 2016 to 2017, both RF and SVM showed general decreases in predicted strawberry guava canopy cover along the eastern side and south-central region of BRGNP. However, in 2017 to 2018 the eastern side generally increased in strawberry guava cover class, especially in SVM predictions. The southwest region showed general  increases in SVM and RF predicted canopy cover from 2016 to 2017, but general decreases from 2017 to 2018. Both classifiers, but especially RF, showed the most widespread decreases in strawberry guava classes from 2018 to 2019. This was followed by widespread increases in 2019 to 2020, although RF showed much of this only being +1 class. RF did also show decreases in 2019 to 2020, especially in central areas. From 2016 to 2020, the central region showed the least amount of change, with large areas unchanged using both RF and SVM, apart from 2019 to 2020 when using RF.

Accuracy assessment
Overall accuracy was 92.52% for SVM with kappa of .88, and 97.81% for RF with kappa of .97 when assessed using the held-out ground truth test data (1323 points). After repeated random sampling of the ground truth data 10 times into 70% training and 30% testing data, mean overall accuracy for RF remained higher than SVM at 97.60% with a 95% confidence interval of ±.20%, versus 91.66% ± .34% for SVM ( Figure 11). Mean kappa for RF was also higher than SVM at .961 with a 95% confidence interval of ±.003, versus .862 ± .007 for SVM. The null hypothesis that the distributions of overall accuracies for SVM and RF were equal was rejected at 95% confidence using the Mann-Whitney two-tailed U test (p = .00018).
The confusion matrix for the initial data partition, used in image production and analysis, showed where misclassifications occurred from both the point of view of the map maker (producer) and map user ( Table 2). The producer's accuracies, which show how often Figure 8. Distribution of slope of pixels classified as low/no guava, mixed, and high guava in the 2018 random forest classified image, with mean (×), and box width proportional to number of pixels in the class. Slope calculated from a 30 m resolution SRTM NASA V3 digital elevation map. Mean slope of high guava pixels was 18.60° with a standard deviation of ±11.06°, compared to low/no guava mean of 12.14° with a standard deviation of ±8.18° and mixed mean of 12.50° with a standard deviation of ±8.93°.
reference data was classified correctly on the map, were highest in both RF and SVM for low/ no strawberry guava, followed by high strawberry guava, and lowest for mixed. The user's accuracies, which show how often the classified area is present on the ground, show more variation between classifiers. In SVM, the most reliable classes were almost equally low/no and mixed, with the least reliable being high strawberry guava. In RF, high strawberry guava was the most reliable by a slight margin, followed by mixed, and slightly less for low/no guava. All class and overall accuracies were >80%. Both SVM and RF had K > .88, indicating almost perfect agreement between observed and expected accuracy (Landis and Koch 1977).

Discussion
The overall aim of this study was to evaluate different approaches for classifying IAS strawberry guava using open-access remote sensing to aid research and management of IAS in tropical forests by studying strawberry guava invasion in BRGNP, Mauritius. Results of this study have shown that additional texture features and vegetation indices increased separability substantially of strawberry guava cover classes over using only spectral bands. In addition, RF produced higher accuracy and less variable results than SVM for the studied application, although both were able to classify images with >70% overall and individual class accuracies. Finally, classified imagery was able to show spatial patterns in the distribution of the target IAS, strawberry guava, in BRGP. In the following sections we discuss how our results have addressed our original objectives related to 1) evaluating the importance of additional vegetation indices and texture features, 2) assessing the accuracy of support vector machine (SVM) and random forest (RF) supervised classifications of imagery; and 3) using classified images to gain an understanding of the distribution of the target IAS in BRGNP.

Vegetation indices and texture features
Spectral separability, as evaluated by SFFS using JM, greatly increased when vegetation indices were added to spectral bands and increased considerably again when texture features were added ( Figure 5). This indicates that the probability of greater accuracy also increased with the addition of vegetation indices and texture features, and the overall accuracies of >90% achieved would not have been possible without these additional Figure 11. Boxplots of overall accuracy and kappa coefficient for RF and SVM classification of the 2018 composite with 10 random samples partitioning ground truth data to 70% training and 30% testing. Table 2. Guava canopy cover classification accuracy confusion matrices using SVM and RF classifiers trained on the 2018 composite with ground truth training data, and tested using held-out ground truth test data. Diagonals in italics show producer's accuracy.

Classified (%)
Reference ( features. This is in line with the similar studies of Laurin et al. (2016) and Puletti, Chianucci and Castaldi (2017) who noted vegetation class separability increased with texture features and/or vegetation indices. Using SFFS with JM, ARI1 and ReNDVI were determined to be among the 20 most important bands, with NDVI and SRI among the 6 least important according to how often they were chosen for subsets ( Figure 6). Not only were habitat ranges diverse, but lower classes had more diverse vegetation than high classes, so selection of indices that are sensitive to more specific spectral differences within vegetation types could be expected, like ARI1 and ReNDVI, rather than more generalized SRI and NDVI. The importance of ARI1, which estimates anthocyanin levels, is supported by the study of Luximon-Ramma, Bahorun and Crosier (2003) which detected anthocyanin in wild fruits on Mauritius. Since strawberry guava was fruiting at the time of the study, it is possible that fruitspecific spectral properties assisted class separability. ReNDVI was more important than NDVI likely due to its more sensitive measurement of chlorophyll concentration and thus decrease in noise from background pigments (Gitelson and Merzlyak 1994).
As expected, texture features played an important role in spectral separability of classes with sum average features accounting for 6 out of the 20 selected features (calculated for B2, B3, B4, B7, B8a, and B11) followed by variance with 5 (for B3, B4, B5, B8, and B12), dissimilarity with 3 (for B4, B8a, and B12), contrast with 2 (for B4 and B8a), and entropy and ASM each with 1 (each for B3). The selection of these features to be included in the subset was supported by Hall-Beyer (2017) and Laurin et al. (2016) who examined the importance of GLCM texture features for classification.
Using principle component analysis, Hall-Beyer (2017) determined which texture measures were most important in describing class differences where there were visual edge features within class patches. Specifically, contrast or dissimilarity should be chosen to distinguish classes with edge components (Hall-Beyer 2017), both of which were included in the 20 selected features in this work. This finding supports our field observation of high strawberry guava cover characterized by dense, monospecific stands of mostly uniform height ( Figure 3). Likewise, it was visually observed that as strawberry guava canopy cover decreased, species diversity in the canopy generally increased, likely also leading to increase within-patch edge features. This is also supported by previous studies indicating decreases in BRGNP forest diversity is associated with increases of invasive species density, including specifically strawberry guava (Lorence and Sussman 1986;Virah-Sawmy et al. 2009;Florens et al. 2017).
The study by Hall-Beyer (2017) also suggested choosing an interior feature to discriminate between interior textures of class patches, specifically recommending mean as first choice and correlation as second. The redundancy of mean and correlation noted in Hall-Beyer (2017) supports the selection of 6 sum average and 0 correlation features in this study. Although Hall-Beyer (2017) noted ASM and variance showed little importance in describing unique characteristics, the selection of these two features does agree with Laurin et al. (2016) where ASM and variance increased JM separability of forest types in wet evergreen forests. However, this likely explains why only 1 ASM measure was included and further why this feature ranked relatively low in overall selection ( Figure 6). Still, Hall-Beyer (2017) did note entropy could be chosen when doing a detailed analysis, and entropy and ASM were always contrasted, likely explaining why both entropy and ASM were selected for B3.
None of the 10 spectral bands were selected for the 20-feature subset, possibly because the selected vegetation indices and texture features sufficiently highlighted class differences while minimizing background noise. All 10 spectral bands were represented in either a texture feature or vegetation index calculation. Furthermore, the additional red-edge bands, B5 to B7, were supported in importance for vegetation monitoring by the selection of ReNDVI, B5 variance, and B7 sum average features in the subset.

SVM and RF comparison
After resampling training and testing data 10 times, strawberry guava canopy cover in BRGNP was classified by RF with a mean overall accuracy of 97.60% (±.20% at 95% confidence), and by SVM with 91.66% mean overall accuracy (±.34% at 95% confidence). This supports the suitability of satellite remote sensing methods to adequately map strawberry guava in BRGNP. Both SVM and RF produced high accuracy classifications; however, RF produced slightly higher accuracy and kappa with lower standard error and greater homogeneity of classes. Accordingly, the null hypothesis that the distributions of overall accuracies for SVM and RF were equal was rejected at 95% confidence. Overall accuracy of RF classification was also less sensitive to parameter input values, with default values based on previous research producing satisfactory results. Therefore, our results further support the use of RF over SVM for strawberry guava canopy cover classification in BRGNP using Sentinel-2.
Both classifiers successfully classified the 2018 image with overall and individual class accuracies >70%. The mixed class included the widest range of strawberry guava canopy cover and thus was the least homogenous class. Consequently, misclassifications occurred more often in mixed classes which was like the studies of Hoyos et al. (2010) and Immitzer, Vuolo and Atzberger (2016).
A notable difference between RF and SVM was in the mixed class where RF had 14% greater producer's accuracy (Table 2). In the mixed class, SVM showed high error of omission (pixels in reference data misclassified in predicted data) to high strawberry guava, leading to the lower producer's accuracy. This is observable in Figure 8 where areas classified as mixed in RF were more likely to be high strawberry guava in SVM. This indicates SVM underestimated the mixed class and overestimated high strawberry guava. This could also explain the discrepancy between SVM and RF predicted total class areas from 2016 to 2020 (Figure 9) where RF predicted bigger yearly changes in total mixed and high strawberry guava area than SVM where total mixed and high strawberry guava remained more stable.
The mixed class included the largest variation of the species of interest, 25% to 75% strawberry guava canopy cover, and was therefore the least uniform class. Furthermore, mixed was the rarest class with only 20% of total pixels used for training and testing. Since SVM only uses support vectors, or boundary conditions, in classification it has been observed to be more sensitive to training data quality and rare classes than RF (Lawrence and Moran 2015;Maxwell, Warner and Fang 2018), possibly leading to the lower producer's accuracy of mixed. Alternatively, because of the ensemble method of RF, it is particularly robust in the presence of noise and outliers (Rodriguez-Galiano et al. 2012;Adam et al. 2014). This likely played a role in the ability of RF to produce more homogenous class regions with less confusion between classes. The greater homogeneity expressed in RF classified images makes it easier to interpret the general distribution of classes. This is desirable for park management applications as it is easier to ascertain where large areas of high strawberry guava canopy cover may be found, and therefore easier to allocate resources.

Strawberry guava distribution in BRGNP
Although areas in BRGNP are consistently being cleared of strawberry guava, the results of manual removal over the entire 66 km 2 area in 5 years is not easily observed from total class areas ( Figure 9). However, when examined spatially, small changes can be observed between 2016 and 2018 and relatively larger changes between 2018 and 2020 including larger areas with decreases in strawberry guava canopy cover ( Figure 10). In 2016 through 2020, areas consistently predicted to be low/no strawberry guava did include areas currently under restoration for invasive species such as the southwest region and northwest of Mare Longue (Ragen et al. 2017).
General distribution of classes was similar from 2016 to 2020 (Figure 7). Classified images of both SVM and RF predicted the highest concentration of strawberry guava with the least amount of change in the central region. This is in line with field observations of dense guava along the southern end of the Parakeet Trail and central to western side of the Paille-en-Queue Trail, where it was common to see people picking fruit. Strawberry guava, including wild-harvested, is a common ingredient in food and drink on the island, possibly making control in these areas more difficult due to cultural significance and seed dispersal by humans.
When observed in relation to slope, high strawberry guava canopy cover was associated with areas of steep slopes more than mixed or low/no strawberry guava areas (Figure 8). Previous studies in BRGNP forests showed strawberry guava occurred readily in the understory of upland forests (Vaughan and Wiehe 1941;Florens et al. 2016Florens et al. , 2017, which includes the eastern edge of the park north of the Savanne Trail. However, much of this upland area was predicted to be low/no strawberry guava canopy cover ( Figure 12). This implies strawberry guava invasion in these areas has been restrained to the understory or removed (Florens et al. 2016(Florens et al. , 2017. The low level of predicted strawberry guava canopy cover in the upland forests on the eastern edge, and higher levels on steep slopes is likely a combination of management practices and influence of forest disturbance. Since control of strawberry guava is limited mostly to physical removal, it is substantially more difficult and dangerous to remove on the steep slopes of the central region of the park so it reasonable that efforts have been concentrated in flatter upland areas. Furthermore, it has been supported that strawberry guava tends to invade more readily in disturbed areas (Virah-Sawmy et al. 2009;Tng et al. 2016), including specifically steep slopes in BRGNP that are more vulnerable to disturbance from cyclones (Virah-Sawmy et al. 2009). Since there have been no cyclones strong enough to cause substantial forest disturbance from 2016 to 2020 (MMS (Mauritius Meteorological Services) 2020b), it would be expected that these steep sloped areas remained relatively stable in strawberry guava cover during this time. Change detection in this study predicted steep-sloped regions to be mostly unchanged from 2016 to 2018 (Figure 10), although decreases in predicted high strawberry guava canopy cover in some steep-sloped central areas could be observed in 2019 and 2020.

Project limitations and future research
This study was limited by the lack of a single cloud-free image for classification. Although mosaiced images are commonly used for classification, misclassifications could have been introduced because of uneven data, especially in composites from fewer images. Visual inspection showed some relationships between areas where the images were mosaiced to fill cloudy areas and class boundaries in classified images. This observation implied a misclassification, possibly from misrepresentation of reflectance values or texture between pixels on either side of the cloud-masked area. Composite images from larger collections or those with fewer cloudy pixels are less likely to succumb to this error since the median pixel value would be drawn from a larger sample. The 2016 image, which was a composite of only 9 images, likely had the largest misclassifications induced by cloudmasking. However, the low standard error of the 2018 image overall indicated these errors were low. Since no ground truth testing data were available for images outside field survey dates, assessment of 2016, 2017, 2019, and 2020 classification accuracy was based on field observations, expert knowledge, and previous studies. It should be noted that error is compounded in the post-classification change comparison (Singh 1989;Warner, Almutairi and Lee 2011), so error of land cover change images was greater than that of the 2018 image.
Gathering more detailed canopy percentage and auxiliary data (i.e. other dominant species and/or stem counts) was not feasible but could have been useful for a more detailed classification and evaluation of errors. For instance, future studies could gather more detail to delineate other dominant species and test classification of more species classes. Additionally, more detailed canopy density measurements could facilitate testing RF regression, rather than RF classification used in this study. Regression models have shown promise for predicting canopy density (Barbosa et al. 2016), but the level of canopy density detail in this study was not sufficient to apply such models.
Although this study supported classification of strawberry guava canopy cover, plants occurring in the understory were not examined resulting in an underestimation for BRGNP. Strawberry guava also occurs in the understory of BRGNP (Vaughan and Wiehe 1941;Florens et al. 2016Florens et al. , 2017, which poses a harder to detect threat to biodiversity.
Lastly, future research using open-access satellite imagery to map strawberry guava in invaded forests of other global regions would be valuable to determine if these results could be more widely applied. Due to the high impact of strawberry guava invasion to native forests discussed earlier, and the global extent of invasion, validated maps of strawberry guava in additional regions could be extremely beneficial for global management.

Conclusions
This study has demonstrated that classification of IAS in tropical forests using moderate resolution open-access imagery from Sentinel-2 MSI with derived texture features and vegetation indices can achieve high classification accuracy using SVM and RF machine learning algorithms. Using the 30% of ground truth data held-out from training and parameter tuning (1323 samples), overall accuracy was 92.52% for SVM and 97.81% for RF. When resampling of training and testing was done, mean overall accuracies were similarly 91.66% (±.34%) for SVM and 97.60% (±.20% at 95% confidence) for RF. Given the choice between SVM and RF for strawberry guava canopy classification in BRGNP, this study recommends RF for its higher accuracy, lower variation, and more homogenous classified regions. This study also strongly supports the addition of carefully selected texture features and vegetation indices for increasing class separability of vegetation species.
Using Sentinel-2 imagery, this research was able to predict, with >80% classification accuracy, that areas of high strawberry guava canopy cover in BRGNP were heavily concentrated in the central regions with lower strawberry guava canopy cover along the eastern border and in the southwest. Furthermore, a significant positive association was identified between areas of high strawberry guava canopy cover and steep slopes. Although total high strawberry guava canopy cover did not seem to be decreasing substantially from 2016 to 2020 despite removal efforts, this is likely due to the short time span of imagery, large area analysed, and possible reinvasion indicated by changes in spatial distribution. This study indicates that as the Sentinel-2 legacy gets longer, it may be feasible to observe significant changes in strawberry guava canopy cover in BRGNP. Furthermore, the inbuilt image catalogues and open-access nature of GEE make continued observations of BRGNP practical and efficient using the methods described here.