Mapping the fractional coverage of the invasive shrub Ulex europaeus with multi-temporal Sentinel-2 imagery utilizing UAV orthoimages and a new spatial optimization approach

Mapping the occurrence patterns of invasive plant species and understanding their invasion dynamics is a crucial requirement for preventing further spread to so far unaffected regions. An established approach to map invasive species across large areas is based on the combination of satellite or aerial remote sensing data with ground truth data from fieldwork. Unmanned aerial vehicles (UAV, also referred to as unmanned aerial systems (UAS)) may represent an interesting and low-cost alternative to labor-intensive fieldwork. Despite the increasing use of UAVs in the field of remote sensing in the last years, operational methods to combine UAV and satellite data are still sparse. Here, we present a new methodological framework to estimate the fractional coverage (FC%) of the invasive shrub species Ulex europaeus (common gorse) on Chiloé Island (south-central Chile), based on ultrahigh-resolution UAV images and a medium resolution intra-annual time-series of Sentinel-2. Our framework is based on three steps: 1) Land cover classification of the UAV orthoimages, 2) reduce the spatial shift between UAV-based land cover classification maps and Sentinel-2 imagery and 3) identify optimal satellite acquisition dates for estimating the actual distribution of Ulex europaeus. In Step 2 we translate the challenging co-registration task between two datasets with very different spatial resolutions into an (machine learning) optimization problem where the UAV-based land cover classification maps obtained in Step 1 are systematically shifted against the satellite images. Based on several Random Forest (RF) models, an optimal fit between varying land cover fractions and the spectral information of Sentinel-2 is identified to correct the spatial offset between both datasets. Considering the spatial shifts of the UAV orthoimages and using optimally timed Sentinel-2 acquisitions led to a significant improvement for the estimation of the current distribution of Ulex europaeus. Furthermore, we found that the Sentinel-2 acquisition from November (flowering time of Ulex europaeus) was particularly important in distinguishing Ulex europaeus from other plant species. Our mapping results could support local efforts in controlling Ulex europaeus. Furthermore, the proposed workflow should be transferable to other use cases where individual target species that are visually detectable in UAV imagery are considered. These findings confirm and underline the great potential of UAV-based groundtruth data for detecting invasive species.


Introduction
Invasive species pose a severe threat to ecosystems across the world (Huang & Asner, 2009). Particularly in areas of high biodiversity, the introduction of competitive foreign species can endanger local flora and fauna and affect important ecosystem services (Mittermeier et al., 2011). Central Chile has been declared a global biodiversity hotspot with a high fraction of endemic species and is listed in six out of the nine common templates for global conservation priority regions (Brooks et al., 2006). This status is contrasted with an ongoing history of dramatic land-use changes, i.e. the replacement of native forests with plantations of exotic tree species (Echeverria et al., 2008) and the presence of a large number of alien and invasive species (Fuentes et al., 2013). In consequence, the region has been reported to be one of the most threatened biodiverse ecosystems in the world (Dinerstein et al., 1995).
One of the most prominent invasive species in south-central Chile is Ulex europaeus (U. europaeus), also known as common gorse (Norambuena et al., 2000). U. europaeus is a nitrogen-fixing, perennial, evergreen shrub, native to the Iberian Peninsula. It has a noticeable yellow flower during spring time and possesses conspicuous spines (Clements et al., 2001). The ability of U. europaeus to quickly invade disturbed and open areas including pastures, juvenile forests and rural roads, forces local people to constantly remove U. europaeus from these areas to maintain forest plantations, infrastructure and agricultural activities (Norambuena & Escobar, 2007). As the removal of U. europaeus after establishment is highly challenging, we hypothesize that a good knowledge of the current distribution as well as a sound understanding of the invasion dynamics of U. europaeus is a crucial requirement for preventing further spread to so far unaffected regions. With this hypothesis we follow the recent study of Altamirano et al. (2016) who investigated how landscape heterogeneity affects the distribution of U. europaeus in the Los Ríos Region in Chile, based on a change detection of two Landsat scenes from 1986 and 2003. The study reported an overall decrease of U. europaeus cover in the study area, mostly due to expanding Pinus radiata plantations. U. europaeus was found to be stable on sites with more fires and the spread potential was positively related to the distance from seed sources and landscape heterogeneity. Although this Landsat-based spatial analysis provided initial information on the invasion dynamics and distribution of U. europaeus in Chile, the exact pattern of the current distribution is unknown. Hence, the aim of this study is to develop an efficient and transferable approach to map the current distribution of U. europaeus so that it can be used by different institutions, e.g. local or national nature conservation authorities. For that purpose, we combine unmanned aerial vehicle (UAV) data and freely available Sentinel-2 data. Whereas, older occurrences of U. europaeus can cover larger areas with dense patches, newly established patches might not fully cover a complete Sentinel-2 pixel. Hence, we focus on the fractional coverage (FC%) of U. europaeus to detect also the most recent spreads.
In our approach we attempt to make optimal use of phenological information provided by Sentinel-2 and the UAV data. We consider the yellow flowering of U. europaeus as a key trait to distinguish it from other plant types (Altamirano et al., 2016;Müllerova et al., 2017). The application of intra-annual time series of optical satellites to improve vegetation classifications has been used in several fields including for example the mapping and classification of NATURA 2000 habitats (Fenske et al., 2020;Schmidt et al., 2014;Stenzel et al., 2017;Franke et al., 2012) or inventorying agricultural landscapes Van Niel & Vicar, 2004;Alcantara et al., 2012;Murakami et al., 2001), forest types (Key et al., 2001) and land cover (Carrao et al., 2008). Given the increasingly high temporal resolution of recent sensor systems, the identification of suitable acquisition dates to optimize the efficiency of mapping tasks (i.e. high classification accuracies applying a minimal number of satellite scenes) becomes an important task. Several earlier studies compared the performances of individual as well as combinations of satellite scenes acquired over the course of a year (Schmidt et al., 2014;Carrao et al., 2008).
In traditional remote sensing studies (especially high-resolution and multi-temporal studies), the reference data were typically collected during labor-intensive and elaborate field campaigns or were derived from other very-high resolution remote sensing data. UAV systems are an interesting and low-cost alternative data source (Tsai & Lin, 2017) when combined with conventional digital cameras (Torres-Sanchez et al., 2014). These systems are now available as off-the-shelf products that can be easily configured and operated after a short training phase. Although the spectral resolution of UAV-based RGB-orthoimages is rather low, the spatial resolution can capture more detail in the spatial ranges (Adão et al., 2017). The ultra-high spatial resolution (<GSD of 10 cm) of data collected with UAV systems and their temporal flexibility is superior to traditional aerial sensor platforms (Turner et al., 2014). The potential of UAV-based monitoring of invasive species was already proven by several studies Martin et al., 2018;Alvarez-Taboada et al., 2017;Müllerová et al., 2017;Michez et al., 2016, Müllerová et al., 2013. Kattenborn et al. (2019) have additionally shown that UAV-based reference data have a high potential to supplement or even replace the traditional sampling of presence or absence data in the field.
One challenge for coupling UAV with satellite data is their coalignment. Preprocessing the UAV images to accurate geo-referenced orthoimages remains challenging (Zhuo et al., 2017). Different software solutions are available for processing the single UAV images into orthoimages in a semi-automatic way (Turner et al., 2014). For this process, accurate information of the orientation (φ = roll, θ = pitch, ψ = yaw) and position (x, y and z-coordinates) of the UAV during the flight are crucial (Whitehead & Hugenholtz, 2014). Concerning the georeferencing procedure, Benassi et al. (2017) distinguish between direct and indirect georeferencing approaches. The direct georeferencing approach uses the information of the onboard IMU (inertial measuring unit) and GNSS module (Global Navigation Satellite System) of the UAV, whereas the indirect georeferencing approach additionally includes precise ground control points (GCP) with sub-decimeter accuracy. The GCP can for example be measured with a precise Real-Time Kinematic (RTK) technique or Post-Processed Kinetmatic (PPK) technique (Jurjević et al., 2020). Several studies have proven that indirect georeferencing leads to higher accuracies in the processed orthoimages than direct georeferencing (Jurjević et al., 2020;Padró et al., 2019). However, collecting GCP is time-consuming, expensive and can be difficult depending on the terrain (Padró et al., 2019;Grayson et al., 2018;Turner et al., 2014) and the accessibility of the study area (Zhuo et al., 2017). In comparison, direct georeferencing is a fast and low-cost way to process UAV raw data but may cause larger horizontal and vertical errors (Jurjević et al., 2020;Turner et al., 2014). Conventional on-board GNSS/IMU modules can show horizontal and vertical errors of several meters (2 -10 m) (Aasen et al., 2018;Turner et al., 2014). The usage of high-end and expensive GNSS/IMU (RTK) modules for UAVs, which offer sub-decimeter information, can improve the accuracy of the direct georeferencing (Padró et al., 2019;Grayson et al., 2018;Benassi et al., 2017). However, even solutions using RTK modules on UAVs can show spatial inaccuracies in a decimeter range due to camera triggering uncertainties or offsets between camera and GNSS RTK position (Ekaso et al., 2020). Woo et al. (2018) propose a method to derive UAV-based GCP without a GCP field survey. Based on a statistical point estimation method, GCP for clearly identifiable urban objects were generated. However, this method seems not to be applicable for UAV flights over rural areas with fewer or even no clear infrastructure elements that could serve as GCP. Another way of increasing the accuracy of directly georeferenced UAV orthoimages is the co-registration to other high-resolution aerial or satellite imagery (Zhuo et al., 2017). Unfortunately, high resolution aerial or satellite images are costly as well and may not be available everywhere. Linking direct georeferenced UAV orthoimages with freely available satellite data as a reference is quite challenging as they provide only medium resolution images. The large differences in spatial grain (pixel size) of UAV data (<10 cm) and medium resolution satellite data (>10 m) complicates the co-alignment of images. A method to automatically correct this misalignment, even for rural areas, would hence be highly valuable.
Here, we address this co-registration issue between moderate resolution satellite imagery and direct georeferenced UAV orthoimages by translating it into an optimization problem where varying sets of FC% of land cover classes were regressed against the spectral information of the satellite images. The variation in the FC% of land cover classes is caused by systematically shifting the position of the underlying UAV-based land cover maps against the satellite image whose position is fixed. It is assumed that the position of the UAV where the overall best fit (lowest model error) between all FC% and the spectral information occurs, indicates the best co-alignment between satellite and UAV data. It is assumed that an optimal co-alignment of the UAV and satellite data will improve the mapping of the targeted invasive species U. europaeus.
The main hypothesis of this study is that the actual distribution of the invasive shrub species U. europaeus can be mapped by utilizing low-cost UAV RGB-orthoimages and freely available Sentinel-2 imagery without additional ground truth data (as reference for mapping and georeferencing). We will demonstrate that i) UAV data can be used as alternative to traditional field information, ii) a spatial optimization approach can improve the spatial accuracy of direct georeferenced UAV orthoimages and iii) using optimally timed satellite acquisitions can improve the differentiation of U. europaeus from other plant or shrub types.
In summary, we present a new methodological framework to estimate the FC% of U. europaeus on Chiloé Island (south-central Chile), based on ultra-high-resolution UAV images and a medium resolution intra-annual time-series of Sentinel-2. Our framework is based on three steps: 1) Land cover classification of the UAV orthoimages, 2) reduce the spatial shift between UAV-based land cover classification maps and Sentinel-2 imagery and 3) identify optimal satellite acquisition dates for estimating the actual distribution of U. europaeus by combining all available Sentinel-2 images.

Study area
Our study focuses on Chiloé Island in south-central Chile (Los Lagos Region) (Fig. 1). The area has a strong oceanic influence with corresponding high humidity and rates of precipitation and low temperature amplitudes (annual average rainfall 2,090 mm, mean annual temperature 12 • C). Chiloé Island covers an area of approximately 9,180 km 2 and its main vegetation-stocked land cover classes include native forests (66.9%) and agricultural land (27.4%) (Barrena et al., 2014). The native forests consist of old-growth north-Patagonian rainforests (Alvaro et al., 2004). In the uplands (maximal elevation 800 m) of Chiloé Island these forests are accompanied by larger areas of Magellanic moorlands (Ruthsatz & Villagran, 1991). The current landscape of Chiloé Island is divided into continuous forest areas (mostly in the southern and western parts of the island) and areas showing mosaics of forest patches and agricultural fields and pastures. The latter was shaped by a history of fire-clearances in some areas of the forest and selective logging in other areas (Barrena et al., 2014).

Sentinel-2
Chiloé Island can be fully covered by four Sentinel-2 tiles (18GWU, 18GXU, 18GWT, 18GXT). For 2016, all images with a cloud cover of <20% were downloaded as Level-1C (top of atmosphere reflectance)  products from the Copernicus Open Access Hub (Table 1). To obtain Level-2A (bottom of atmosphere reflectance) products and to detect clouds and shadows, the Sen2Cor stand-alone algorithm (Li et al., 2018) was used. In a next step, clouds and shadows were detected with the classification approach integrated in Sen2Cor, masked out and resulting gaps were linearly interpolated over the entire time-series. We processed the data at 20 m resolution to retain the consistent original spectral information of the bands. Furthermore, this enables the comparability to other medium resolution satellite sensors, like Landsat-8 (30 m resolution). In this study, the 60 m bands (B1, B9, B10) and the 10 m broad NIR band (B8) were removed (Table 2). Following the study of Korhonen et al. (2017) the broad NIR band (B8) was excluded, as its main advantage is the high spatial resolution, which was removed by resampling to 20 m resolution. In addition, the 20 m narrow NIR band (B8a) has a much smaller bandwidth, which is more sensitive to plant chlorophyll and seems to be better suited for the calculation of vegetation indices (Wittke et al., 2019). For each date the four tiles (100 km*100 km) were merged into one orthoimage. Additionally, we used the NDVI (Tucker, 1979) to increase the spectral information of the S2 data stack. Including the NDVI, a total of ten spectral features were available for each date.

UAV orthoimages
Ten ultra-high resolution UAV orthoimages were recorded during a two-week field campaign from 20th of November to 2nd of December 2016 on Chiloé Island. In the period between October and November U. europaeus has an intense yellow flower but at the time of the field campaign, the flower had already started to wither. Nevertheless, even with withered flowers it was possible to reliably distinguish U. europaeus from all other shrub species on the island. The flight areas were preselected on Google Earth using WorldView images in which larger patches of U. europaeus were clearly visible in the last years due to their yellow flower. The flight areas were distributed over the agricultural area of Chiloé Island, mainly in the north and east, as U. europaeus has successfully spread there. Each flight covered an area of approximately 0.5 km 2 .
The UAV system used consisted of an octocopter (OktoXL -HiSystems GmbH), a GNSS navigation module with compass, an inertial measurement unit (IMU) and a gimbal (gyroscope) with an attached Canon EOS 100D camera with a 18-55 mm lens. The total weight of the system was about 5 kg. We used two 4,400 mAh six-cell LiPo batteries, resulting in a flight time of approximately 10 to 15 min. Depending on the local terrain, the flight altitude was between 90 m and 160 m (Table 3). To obtain gapless orthoimages we used side and forward overlapping of 70%. Based on these flight settings, a total of 500 to 1000 images were captured per flight. The ground sampling distance (GSD) of the UAV images was about 2-3 cm per pixel but was aggregated to consistent 10 cm resolution during the direct georeferencing process (using orientation and positioning information of the on-board GNSS/ IMU modules) to account for smaller distortions. This task was done with the AGISoft Photoscan Professional software (AgiSoft PhotoScan, 2016), since several studies (Fraser & Congalton, 2018;Turner et al., 2014) have shown its high potential for generating orthoimages from single UAV images. The range of the number of tie-points for the image alignment task was between 115,000 and 450,000. The total camera error of the UAV orthoimages was about 1 m on average with an average reprojection error of 2.87 pixel.

Spectral and 2D textural features
For improving the land cover classification of the UAV orthoimages (see section III), we included several spectral and 2D textural features in this study. We included RGB-based vegetation indices (VI vis ) as no NIR information was available. Some studies have already shown the potential of UAV-based VI vis (Lussem et al., 2019;Viljanen et al., 2018;Lussem et al., 2018), but mainly for above-ground biomass estimations. We used the Visible Atmospherically Resistant Index (VARI) (Gitelson et al., 2002), Green Leaf Index (GLI) (Louhaichi et al., 2001) and Normalized Green Red Difference Index (NGRDI) (Tucker, 1979) as tested in Lussem et al. (2018) (Table 4).
2D textural features can reduce noise effects caused by very small objects, detectable due to the ultra-high resolution of UAV orthoimages (Kwak & Park, 2019). The gray-level co-occurrence matrix (GLCM) approach (Haralick et al., 1973) is a very common way to reduce this problem and was successfully applied for invasive species monitoring , Michez et al., 2016 and crop classification

Louhaichi et al., 2001
Normalized Green Red Difference Index NGRDI Tucker, 1979 (Kwak & Park, 2019). We used the glcm package (Zvoleff, 2020) in R (version 1.6.4) to calculate the mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment and correlation for all directions and for the red band of the orthoimages. All formulas for the eight texture features can be found in the original work of Haralick et al. (1973). The red band was chosen because it shows the highest contrast for vegetation and in order to avoid the typical auto-correlation amongst predictors when more than one band is used (Feng et al., 2015). Kattenborn et al. (2019) used ten different kernel sizes (from 0.25 to 4 m) to account for the different leaf forms and branch structures of the different invasive species. Feng et al. (2015) found that an optimal kernel size for classifying urban vegetation is around 2 m. We only used 2 kernel sizes (0.5 and 1 m), to prevent an aggregation of smaller individuals of U. europaeus. Finally, a data stack with 22 features (3 visible bands, 3 VI vis , 16 texture features (kernel size 0.5 m and 1.0 m)) was available for each UAV flight.

Methods
For mapping the actual distribution of U. europaeus, we implemented a three-step workflow (Fig. 2). Steps 1 (supervised classification of UAV images) and 2 (spatial optimization between UAV and Sentinel-2 data) were accomplished for each UAV orthoimage individually (UAVlevel).
Step 3 then applied the fractional coverage (FC%) of U. europaeus of all corrected UAV-based land cover maps (S2-level). The primary goal of step 1 was to semi-automatically derive reference data (FC%) from the UAV orthoimages by using a Random Forest-based land cover classification with up to eight classes for the UAV orthoimages based on manually delineated reference polygons.
In the spatial optimization approach (step 2), the land cover maps were systematically shifted against the fixed Sentinel-2 (S2) images. For each shift, the changing FC% of each land cover class were derived for the S2-pixels and were included in several Random Forest (RF) (Breiman, 1999) models (one per class) to estimate the FC% per S2-pixel per land cover class. The shift position with the lowest average model error per UAV orthoimage indicates the real position of the UAV orthoimage to the S2 data. All ten U. europaeus maps were then corrected according to these detected shifts. The reason for considering several classes in this step is that correcting or mitigating spatial shifts based on only one class would be strongly biased to the occurrences of this class and the error of the utilized model trying to estimate these occurrences. We assume that if more land cover classes are included, this bias would be averaged over all classes (models) and the detected shifts will be more representative for all land cover classes and hence the image.
In step 3, FC% for U. europaeus per S2 pixel were derived based on the corrected UAV-based land cover maps. The FC% of U. europaeus were separated into training and validation data, according to the UAV orthoimages. Seven UAV-maps were used for predicting the distribution of U. europaeus on Chiloé Island. The remaining three maps were used to assess the accuracy. Finally, to test which S2 acquisition had the highest impact on the predictions, we then tested all possible acquisition combinations.

Land cover classification (Step 1)
Based on the ultra-high resolution UAV orthoimages, reference data for the land cover classification were derived (UAV-level). For this task we manually delineated polygons of the main land cover classes (Table 5) of each UAV orthoimage. The reference data were crosschecked with high resolution images from Google Earth to ensure a correct determination of the classes. The differentiation of U. europaeus from other shrub types was improved due to its yellow flowering during the field campaign. Additionally, our knowledge of the study area and photographs acquired during the field campaign were used to improve our assessment.
Within these polygons 1,000 sample points per class were selected (stratified random selection, with a minimum distance of 0.5 m). The sample points were separated (1:1) into a training and a validation set. A Random Forest (RF) model was trained for each UAV orthoimage including the UAV data stack (RGB bands, VI vis , 2D texture parameters)   (Ramírez et al., 2014;Villagrán, 1988) as predictors. We used the Random Forest implementation of the caret package (Kuhn, 2020) in R (R Core Team, 2020) within the RStudio IDE (RStudio . To optimize the model performance, the mtry parameter (the number of randomly selected predictors) was determined empirically in a grid search function. The number of trees (ntree parameter) was set to the default value of 500, because this parameter has less influence on the model performance and the default value of 500 was successfully used before (e.g., Fernández-Delgado et al., 2014). This setting was used for all Random Forest models in this study. To improve the overall accuracy of the classification, we tested several RF models with different feature group (RGB bands, VI vis , 2D texture parameters) combinations to check which individual feature group and which group combination performs best.

Spatial optimization (step 2)
The purpose of the spatial optimization was to find the best coalignment between the UAV-based orthoimages and Sentinel-2 images. Since each UAV orthoimage had an individual spatial offset to the Sentinel-2 images, all UAV-based orthoimages were considered individually and independently of each other (UAV-level). During the spatial optimization, the position of the UAV-based land cover maps was systematically shifted against the position of the Sentinel-2 images, similar to a moving window approach (Fig. 3). Starting from the basic position (no shift), the maps move in 2 m intervals for a maximum of 20 m in each direction. A total of 441 different shift positions were examined. Three tasks were completed in each shift interval: 1) determining the FC% of each land cover class for the underlying Sentinel-2 pixel (20 m). 2) Train a unique Random Forest (RF) model for each land cover class to estimate the corresponding FC%. The 60 spectral features of the Sentinel-2 pixels were used as predictor variables, the FC% values as response variable. Sentinel-2 pixels which cover the edges of the orthoimages were excluded from the analysis, because these areas normally show distortions due to the lack of overlapping images in the image alignment process. On average about 500 (min 171; max 1012) Sentinel-2 pixels were available per UAV orthoimage for training the RF model. A 5-fold cross-validation (80% of the pixels for training and 20% for validation) was used for assessing the RF model accuracy. The estimated FC% were averaged per pixel.
3) The accuracy of the predictions was assessed with the normalized root mean squared error (nRMSE), because this metric considers the value range of the target variable. This makes it more reliable when comparing models with different scales (Bennett et al., 2013). The nRMSE values of all models for all land cover classes were then averaged per shift position (nRMSE avg ). According to our assumption, the lowest nRMSE avg of the 441 different shift positions indicates the real position of the UAV-based orthoimages and the Sentinel-2 imagery.

Mapping the actual distribution of U. Europaeus (step 3)
For this next step, we used all ten corrected U. europaeus maps to derive and validate the final FC% of U. europaeus for all covered Sentinel-2 pixels (S2-level), except for the pixels lying on the borders of the orthoimages. We separated the FC% data into a training and a validation set. The FC% of seven U. europaeus maps (2989 Sentinel-2 pixel in total) were used for predicting the FC% for Chiloé Island with a Random Forest (RF) model and the FC% of the remaining three maps (1001 Sentinel-2 pixel in total) for validating the predictions. We chose three maps with low, medium and high FC% of U. europaeus as the validation maps. This ensures, that the validation data cover the entire FC% range (0 to 100%).
To check whether the spatial optimization approach improves the model accuracy, we analyzed three different scenarios ( Table 6).
The Wilcoxon-Signed-Rank-Sum-Test (Wilcoxon, 1945) was then used to check for significant differences between the scenarios. This nonparametric hypothesis test allows the comparison of dependent and paired samples, to assess whether their means are from the same or from different populations. Despite both models relying on different subsets with slightly different cover fractions, it can be assumed that the cover fractions have a strong correlation, since they were only a few meters apart or even overlapped.

Identify optimal Sentinel-2 acquisition dates (step 3)
Finally, to further improve the model accuracy, we searched for the  Table 6 Tested scenarios.

ID Scenarios Explanation
A) no shift no spatial optimization was applied B) best shift shift position with the lowest nRMSE C) worst shift shift position with the highest nRMSE optimal Sentinel-2 acquisitions to detect the actual distribution (FC%) of U. europaeus. For this task, we tested all six available S2 acquisitions (including all spectral features) in all possible combinations. This approach was successfully used in Schmidt et al. (2014) to search for the most important RapidEye acquisitions to classify semi-natural grassland.
According to the binomial coefficient (Equation I) a total of 63 unique date combinations were possible without considering the order of the acquisitions. For each model run the nRMSE was calculated to assess the accuracy. The influence of a single acquisition depends on the total number of used images (Schmidt et al. 2014). The more acquisitions were used in a model, the lower the possible influence of a single acquisition. Therefore, we grouped all examined combinations according to the number of acquisitions used to find the influence of each acquisition for these groups. For each acquisition, all combinations of a specific group (e.g. all combinations with 3 acquisitions) were further separated into a sub-group A) where the acquisition was included and a sub-group B) where the acquisition was excluded. For both sub-groups the average nRMSE was calculated and then subtracted (nRMSE diff ) (Equation II). The direction of the nRMSE diff indicates whether the acquisition has a negative or positive impact. (1)

Land cover maps (UAV-level)
The overall accuracies (OA%) of the land cover maps for the different tested feature groups (RGB bands [RGB], VI vis [VI], 2D texture parameters [GLCM]) combinations and UAV flights vary between 0.55% and 0.94%. To find the best feature combination, we averaged the OA% of each UAV flight for the different feature groups. When comparing the feature groups individually, GLCM shows the highest average accuracy (0.84 OA% and a standard error (SE) of 0.04), followed by RGB (0.77 OA %; 0.05 SE) and VI (0.65 OA%; 0.06 SE) (Fig. 4). The feature combination with the highest average accuracy (0.91 OA%; 0.02 SE) includes all three feature groups (RGB-VI-GLCM). Therefore, we used this feature group combination as reference data for all following analyses. A complete overview of the feature importance of the Random Forest model including all feature groups is in the appendix section ( Figure A-1).

Spatial optimization (UAV-level)
The accuracy of the Random Forest models estimating the FC% for the land cover classes were averaged per shift position for each UAV orthoimage (nRMSE avg and R2 avg ). Thus, 441 different nRMSE avg and R2 avg values were obtained for each UAV orthoimage. The range of the nRMSE avg of all UAV orthoimages was between 0.092 and 0.219. An even larger difference could be detected for the R2 avg where the range was between 0.154 and 0.799 (see appendix section). When comparing the individual UAV orthoimages, differences in the accuracy ranges can be identified. The largest range was found for UAV orthoimage #10 (Δ0.077 nRMSE avg ; Δ0.477 R2 avg ) and the lowest for #5 (Δ0.04 nRMSE avg ; Δ0.28 R2 avg ).
For a better interpretation of the results, we visualized the nRMSE avg values of all 441 shift positions of all UAV flights as heatmaps. In Fig. 6, exemplary heatmaps for the UAV orthoimages #7 and #8 are shown. A complete overview of all heatmaps is in the appendix section ( Figure A-2 and A-3).
Most heatmaps, clearly show a single hotspot (reddish colors) reaching lowest errors (nRMSE avg ), while areas further away from this hotspot achieve constantly lower accuracies. In Table 8, the model errors of the three scenarios (Ano shift, Bshift with lowest error, Cshift with highest error) and their detected x-and y-shifts are shown. The largest improvements (scenario B compared to scenario A) were detected for UAV orthoimage #9 (-0.02 nRMSE avg ; +0.072 R 2 avg ) and the lowest for #2 (-0.003 nRMSE avg ; +0.028 R 2 avg ). The largest deteriorations of the model error (scenario C compared to scenario A) were detected for UAV orthoimage #10 (+0.077 nRMSE avg ; − 0.477 R 2 avg ) and the lowest for #5 (+0.04 nRMSE avg ; − 0.28R 2 avg ). The largest absolute shift (absolute sum of x-and y-shift) of scenario B was detected for UAV orthoimage #9 (|18| m) and the lowest for #7 (|4| m). On average, the detected shift of scenario B was about |2.8| m in the x-direction and about |7.8| m in the y-direction. The average shift of scenario (C) was | 19.6| m in x-direction and |18.6| m in the y-direction. All UAV-based U. europaeus maps were then corrected according to the detected shift positions of scenario B.
We provide an additional plausibility check for these results in the Supplementary material.

Detecting the actual distribution of U. Europaeus
Based on the seven pre-selected and corrected U. europaeus maps, we trained a Random Forest (RF) model to estimate FC% of U. europaeus. We tested all possible S2 acquisition date combinations as input variables for the RF model. The accuracy of the estimations was assessed with three hold-out maps.
The nRMSE values of all acquisition date combinations ranged between 0.138 (0.478 R 2 ) and 0.108 (0.68 R 2 ). Fig. 7 shows the nRMSE and the R 2 values of all date combinations grouped by the number of acquisitions used in the RF model. On average, the model accuracy was better when more images were used (Table 9). The average nRMSE of all mono-temporal models was 0.126 with an average R 2 of 0.562. In comparison, all models including two acquisition dates reached an average error of 0.117 nRMSE and an average R 2 of 0.631. When using all six S2 acquisition dates, an error of 0.108 nRMSE and a R 2 of 0.68 was reached. However, the gain of accuracy leveled off with each added image. Including a random second acquisition in the RF model, the improvement was − 0.09 nRMSE (+0.069 R 2 ) on average, while the inclusion of a random sixth acquisition improved the model by only − 0.002 nRMSE (+0.008 R 2 ). When comparing the unique acquisition date combinations, the combination including all six acquisitions achieved the best overall accuracy (nRMSE = 0.108; R 2 = 0.68). Although the accuracy was only slightly higher compared to some other combinations.
Based on the best acquisition date combination, we examined whether there are significant differences in model accuracy between the three scenarios (Fig. 8). Scenario A (nRMSE = 0.127; R 2 = 0.618) and C (nRMSE = 0.171; R 2 = 0.312) show much lower accuracies than scenario B (nRMSE = 0.108; R 2 = 0.68). In scenario B the overprediction of U. europaeus with low FC% as well as the underprediction for high FC% from scenarios A and C is clearly reduced. A Wilcoxon-Signed-Rank-Sum-Test indicates that the FC% predictions of the three scenarios differ significantly (p-value < 0.05).
In Fig. 9, the FC% estimations of U. europaeus for Chiloé Island based on scenario B are shown. For a better visualization, the FC% were aggregated to a pixel size of 500*500 m where each pixel shows the maximal value of all pixels contained in the area. The map shows two main hotspots of U. europaeus occurrences on Chiloé Island: A) the area around the city of Ancud including the northern shoreline and along the main route "Ruta 5" heading southwards and B) the densely populated area around the cities of Castro, Chonchi and Dalcahue. The islands between the Gulf of Ancud and the Gulf of Corcovado are also infested by U. europaeus. Furthermore, it can be shown that U. europaeus occurrences are lower in the southern and the western parts of the island.
We can see that also smaller fractions of U. europaeus were detected properly based on the S2 data (Fig. 10). The results for UAV orthoimage #3 in particular indicate that smaller patches along the roads and at the edges of pastureland could be correctly identified by the RF model. Larger patches of U. europaeus located on UAV orthoimages #1 and #9 were also clearly differentiated from the surrounding vegetation. We also observed some false occurrence predictions with low FC% values within dense shrub areas, where no U. europaeus was observable in the UAV orthoimages. Almost no FC% were estimated within pasturelands and arable lands.

Acquisition date importance
In Fig. 11, the nRMSE and the R 2 values for all combinations where the single acquisitions were included or excluded are shown. The results suggest that only the acquisition from 2016 to 11-20 (spring) has a high   positive influence on the model predictions. On average, the model performance was more accurate when the November acquisition was included (improvement of − 0.008 in nRMSE; and + 0.044 R 2 ). The other five S2 acquisitions show only slightly positive influences on the model accuracy. The S2 acquisitions from 2016 to 03-05 and 2016-07-03 improve the model accuracy by approximately − 0.004 nRMSE (both) and + 0.025 and + 0.027 R 2 , respectively. The remaining three S2 acquisitions (2016-01-05, 2016-07-23 and 2019-09-01) improve the model accuracy only very slightly (<-0.001 nRMSE; < +0.02 R 2 ).

Discussion
Mapping the spread of invasive species requires spatially and temporally detailed information. In data from medium spatial resolution sensors this information is likely to be on the sub-pixel level while very high spatial resolution EO data can provide this information directly (Immitzer et al., 2016). Since VHR remote sensing data from satellites or aircrafts is still costly and often not available in time-series or at a specific pre-defined acquisition date (e.g. flowering), the combination of low-cost UAV data and freely available satellite data can play an important role for detecting invasive plant species in early stages and across larger spatial extents. Here, our primary objective was to develop an operational methodological framework for using ultra-high resolution UAV-based RGB-orthoimages and medium resolution Sentinel-2 data for mapping the distribution of woody invasive species, using the example of U. europaeus on Chiloé Island. The simultaneous use of ultrahigh resolution UAV data and moderate resolution satellite data is limited by the challenging spatial co-registration of the data-sets due to their notably differing spatial grain. In this study we addressed the problem by developing a spatial optimization approach to automatically mitigate spatial shifts between the UAV and the satellite data.

Land cover classification
Our analysis shows that manually delineated major land cover types within the UAV orthoimages can be used as reference data for a proper land cover classification. Similar successful results were found in studies focusing on the classification of crop types (Kwak & Park, 2019), settlements (Gevaert et al., 2016) and land cover (Alvarez-Taboada et al., 2017). Here, the high overall accuracy (>0.85 %OA) of all ten classified UAV orthoimages indicates that the detailed information of the UAV orthoimages is sufficient to differentiate the general main land cover types. Even land cover classes that show similar spectral characteristics, as for example native shrubs, peatbogs (Magellanic moorlands) and U. europaeus, were successfully differentiated even if the flowers of U. europaeus were almost wilted. The high producer's and user's accuracy for U. europaeus (>0.78), shrub (>0.77) and peatbogs (>0.78) underlines this finding. For U. europaeus patches without flowers, the accuracy decreases. The importance of an optimal timing of the UAV field campaign is in accordance with the work of Müllerová et al. (2017). Some misclassifications can be explained by shadows, since most of the shadows were falsely classified as water or shrub. This is in accordance with the study of Lopatin et al. (2019). They found that shaded areas were falsely classified with error rates between 65% and 100% even when shadows were included in the calibration process. In contrast to that finding, Al-Najjar et al. (2019) successfully included a shadow class when classifying land cover classes in urban areas. However, the most promising approach to mitigate the negative influence of shadows is to schedule the flight time adopted to the sun elevation angle Döpper et al., 2020).
We observed that high accuracies of the land cover classification were derived only when including textural information in the classification process. However, the usage of textural information leads to a certain degree of aggregation along transitions from one land cover class to another. Depending on the kernel size used to calculate the textural information, the degree of aggregation for smaller patches and transition zones increases. In the study of Alvarez-Taboada et al. (2017), the OA% of a land cover classification dropped when using a kernel size>35 cm (5 pixel). This can be a drawback when searching for juvenile occurrences of invasive plants. In this study, this negative effect is mitigated by using only small kernel sizes. The exclusive use of the spectral information of the three tested RGB-based vegetation indices or of the visible RGB bands results in a lower classification accuracy. Some studies have already proven that the missing NIR information can be a drawback of UAV-based RGB-orthoimages for land cover classifications (Ahmed et al., 2017). However, the usage of additional spectral features improved the land cover classification when they were used in combination with textural information. The benefit of considering the textural information for the classification of UAV-based RGB-orthoimages confirms the findings of Kattenborn et al. (2019) and Feng et al. (2015). Expectedly, the structure of smaller objects such as leaves or branches is more important for differentiating land cover classes than the limited spectral information provided by RGB data.

Utilizing UAV data as reference data
Earlier studies have shown already the potential of UAV-based RGBorthoimages as substitutes for in-situ data from traditional field surveys when estimating fractional coverage (FC%) of target objects Riihimäki et al., 2019). The high overall accuracy of the land cover classification in this study enables a reliable upscaling of FC% for all classes on larger scales. Furthermore, this indicates the high potential of UAV-based RGB-orthoimages as an alternative to using traditional field data to estimate the fractional coverage. Riihimäki et al. (2019) compared FC% of arctic vegetation, derived from a binary classification of UAV RGB-orthoimages, with FC% of traditional field sample plots. The binary classification of vegetation and non-vegetation reached a producer's and user's accuracy of > 0.73 (up to 0.97) for different test sites. The correlation of the derived FC% to the traditional field data was R 2 0.90 with a RMSE of 0.14. Nevertheless, some caution is needed when using UAV-based data for scientific or analytical purposes. Despite careful processing workflows, some issues related to varying illumination conditions (e.g. saturation), spatial artefacts (e.g. blurring) (Manfreda et al., 2018) or shadows  will always remain in UAV orthoimages. These issues can be challenging when subtle details of the target plant (e.g. blossoms) are key for differentiation. Therefore, the influence of unfavorable weather conditions (e.g. sudden wind gusts, changing sunlight situation) should be reduced as much as possible during the flight campaigns. Furthermore, understanding the fundamental flight planning characteristics (e.g. flight altitude, flight velocity, positioning accuracy, image overlapping, sensor technology) and the data processing (e.g. SfM technology) is crucial for processing reliable UAV products (Aasen et al., 2018;Fraser & Congalton, 2018). Furthermore, guidelines on how to create reliable UAV orthomosaics are still few in the literature (Döpper et al., 2020;Manfreda et al., 2018;Fraser & Congalton, 2018) and further research is needed to fully understand the potentials and the issues of UAV-based remote sensing.

Spatial optimization
In cases where no very-high positioning information (sub-decimeter accuracy) are available, an indirect georeferencing solution, based on GCP or high-resolution aerial images, is normally used. If only standard on-board (navigation) GNSS/IMU modules are available, the positional accuracy is typically limited to 2-10 m (Aasen et al., 2018;Turner et al., 2014). Orthoimages derived from such GNSS/IMU modules are inadequate for UAV applications (Aasen et al., 2018). Woo et al. (2018) reduce the errors by deriving GCP from UAV-based RGB-orthoimages directly. The resulting RGB-orthoimage including these GCPs has a total RMSE of approx. 1 m and without approx. 4 m. However, this approach seems to be only applicable for urban areas with clearly visible infrastructure patterns. Limitations in positional accuracy should be addressed, particularly if the UAV orthoimages will be linked with   11. Influence of the Sentinel-2 acquisition dates. All 63 different nRMSE and the R 2 values were grouped according to if an acquisition date was included in the corresponding combination (blue) or excluded (red). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) satellite data. In the study of Aicardi et al. (2016), a method is proposed to co-align multi-temporal UAV orthoimages on one anchor UAV orthoimage without GCPs. Other studies focused on developing new approaches for co-aligning multi-source remote sensing data (e.g. Scheffler et al., 2017;Yan et al., 2016), but methods for linking remote sensing data with VHR UAV data are still missing. Although there are approaches that involve the use of UAV and medium resolution in a joint approach (Riihimäki et al., 2019;Kattenborn et al., 2019), to the best of our knowledge none evaluated the spatial accuracy or attempted systematically to improve the co-registration of ultra-high resolution UAV and medium resolution remote sensing images. Here, we suggest identifying the optimal co-alignment between Sentinel-2 imagery and UAV orthoimages by translating the spatial offset of both sensors into an optimization problem. Using the suggested spatial optimization approach, we found that the average spatial offsets between UAV and satellite data were approx. 3 m on the x-axis and approx. 8 m on the yaxis. These errors are in accordance with the above-mentioned expected positional errors of standard on-board (navigation) GNSS modules. The stronger shift on the y-axis (always in the south direction), could be explained by the systematical geolocation performance of the Sentinel-2A satellite of 12.5 m (Languille et al., 2015). Furthermore, Yan et al. (2016) have shown, that Landsat 8 images also have a spatial misregistration in the south direction of Sentinel-2 images. However, discussing the source of spatial error is beyond the scope of this study.
The range of improvement of the single UAV orthoimages (UAVlevel) was between 0.003 and 0.02 nRMSE and between 0.028 and 0.084 R 2 . A relationship between the strength of improvement and the absolute sum of the spatial shift in the x and y direction is visible. Smaller improvements expectedly indicate smaller spatial shifts and vice versa. The smaller the real spatial shift between UAV and satellite images, the smaller the difference in spectral information of both sensors for the captured pixel or area. This finding is underlined by the trend of the worst-case scenario C. Here, the worst shifts were found in an average distance of approx. 20 m on the x-axis and approx. 19 m on the y-axis which is close to the maximum distance. Taking the spatial offsets of scenario B into account, the prediction of U. europaeus FC% for Chiloé Island (S2-level) was significantly improved by − 0.019 nRMSE and + 0.062 R 2 . The improvement was even stronger when compared to the worst-case scenario C (-0.063 nRMSE, +0.378 R 2 ).
However, it should be noted that the proposed method can only be a workaround for cases where no very-high positioning information (for a direct or indirect georeferencing) or reference data is available. Distortions within the UAV images can neither be detected nor mitigated with the proposed method.

Potential of medium-resolution Sentinel-2 data for detecting FC% of U. Europaeus
Although the spatial resolution of Sentinel-2 (here 20 m) seems to be too coarse for detecting small changes of FC%, we were able to demonstrate a general potential of medium-resolution satellite imagery for mapping the current distribution of the invasive shrub U. europaeus. Concerning the phenological information, our results demonstrate that with an increased number of satellite acquisitions, the mean accuracy of our model predictions can be improved. However, this trend saturates with larger numbers of acquisitions. Similar results have been shown in previous studies dealing with vegetation classification (Schmidt et al., 2014;Schuster et al., 2015;Carrao et al., 2008). We also showed that an optimally timed sparse time-series can reach equal results compared to the combination including all six acquisitions. It was found that the Sentinel-2 acquisition captured during the flowering phase of U. europaeus (2016-11-20) was very important for detecting its actual distribution. Although the bright yellow flower began to wither at the time of this satellite image recording, the phenological information was still unique, as no other flowering shrub was in bloom at that time. It can be assumed, that an image acquired at the peak of the flowering phase would have an even larger positive impact.
These results show that it is important to capture key periods during the phenological development, in which the target species can be easily distinguished from other existing plants. The high repetition rate of Sentinel-2 is highly valuable to increase the chance of collecting suitable data, especially with the now simultaneously operating Sentinel-2 A and B satellites, which were not available during the time of the field campaign. In areas with large periods of cloud coverage like Chiloé Island, a harmonized combination of Landsat and Sentinel-2 (as shown e. g. in Yan et al., 2016) may provide further advantages.

Actual distribution of U. Europaeus on Chiloé Island
As Altamirano et al. (2016) stated in their work, there is a strong need for studies analyzing the invasion of U. europaeus, especially for Chile as a biodiversity hotspot. Despite its categorization as one of the most invasive species in the world, the actual distribution of U. europaeus was never mapped for larger parts of south-central Chile. With the presented approach, we were able to map U. europaeus on Chiloé Island based on Sentinel-2 and ultra-high resolution UAV images. U. europaeus is a very dynamic invasive which in south-central Chile mainly invades scrub and bare land, but also juvenile forest plantations and agricultural land (Altamirano et al., 2016). These findings are in line with our results, where we found that U. europaeus occurs mainly in the non-forested regions of Chiloé Island (Fig. 9). These areas are dominated by agriculture (mainly pastureland) with patches of shrub and forest remnants (Barrena et al., 2014;Aravena et al., 2002). Our results show two hotspots of U. europaeus occurrences on the island. The first one is the area around the city Ancud including the northern shoreline and along the Ruta 5 southwards. The second hotspot is located in the densely populated area around the cities of Chonchi, Castro and Dalcahue including the islands between the Gulf of Ancud and the Gulf of Corcovado. Our results indicate a possible relationship between the presence of U. europaeus and anthropogenic disturbances and activities (e.g forest fires, plantations, distance to towns) which were already pointed out in previous studies (Altamirano et al., 2016;Rees & Hill, 2001). The fact that there are still regions on Chiloé Island without occurrences of U. europaeus, shows the importance of mapping the actual state of the invasion to prevent further spread. While the overall quality of the mapping result was good, there were also some errors in the prediction map, especially in the southern and western parts of the island. Particularly wetlands (along the Chepu River) and Magellanic moorlands in the uplands of Chiloé Island were in some cases falsely recognized as areas with a medium to high fraction of U. europaeus. The Magellanic moorlands on Chiloé Island were dominated by cushion plants (Ruthsatz and Villagran, 1991), which form light-green to yellow cushions. The spectral characteristics of this plant seem to be similar to the spectral profile of U. europaeus at the time of the Sentinel-2 acquisitions. The map also shows a tendency to overestimate low fractions of U. europaeus cover. For areas with almost no U. europaeus coverage, low fractions are estimated. On the other hand, some areas with large fractions were tendentially slightly underestimated which can be explained by the overall tendency of the RF algorithm to overestimate lower values and underestimate higher values (Fig. 8). This is a well-known phenomenon of the RF algorithm (Zhang & Lu, 2012;Lopatin et al., 2016).

Conclusions
In this study, we present a novel approach to estimate pixel-wise fractional coverage (FC%) of the woody invasive species U. europaeus across Chiloé Island (south-central Chile) using a combination of ultrahigh UAV and medium resolution multi-temporal Sentinel-2 data. The core contributions of this study include: i. UAV-based RGB-orthoimages can be used to derive valid reference data (here FC%) for determining the actual spread of invasive shrub species (U. europaeus) using moderate resolution satellite data. The accuracy of the reference data was best when textural information was included in the derivation of the FC%. ii. The proposed spatial optimization approach to co-align UAV and Sentinel-2 data can be used successfully to detect and reduce the spatial offsets between both sensors. The accuracy of estimating the FC% of U. europaeus on Chiloé Island (S2-level) was significantly improved by considering the spatial offsets between both datasets. iii. While multi-temporal Sentinel-2 acquisitions generally improved the model performance for estimating fractional coverage of U. europaeus, a few key acquisition dates (particularly covering the flowering period of U. europaeus) resulted in accuracies similar to those produced by the models based on all available datasets.
Occurrence maps for U. europaeus were obtained for the complete Chiloé Island and showed good agreement with the impressions gathered during our field survey. This, in combination with the timeefficiency of the UAV-based field campaign indicates a high potential of the presented approach to support the mapping and monitoring of U. europaeus in south-central Chile and other parts of the world facing the threat of invasive species with similar properties.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.