Optimizing the Timing of Unmanned Aerial Vehicle Image Acquisition for Applied Mapping of Woody Vegetation Species Using Feature Selection

Most recent studies relating to the classification of vegetation species on the individual level use cutting-edge sensors and follow a data-driven approach, aimed at maximizing classification accuracy within a relatively small allocated area of optimal conditions. However, this approach does not incorporate cost-benefit considerations or the ability of applying the chosen methodology for applied mapping over larger areas with higher natural heterogeneity. In this study, we present a phenology-based cost-effective approach for optimizing the number and timing of unmanned aerial vehicle (UAV) imagery acquisition, based on a priori near-surface observations. A ground-placed camera was used in order to generate annual time series of nine spectral indices and three color conversions (red, green and blue to hue, saturation and value) in four different East Mediterranean sites that represent different environmental conditions. After outliers’ removal, the time series dataset represented 1852 individuals of 12 common vegetation species and annual herbaceous patches. A feature selection process was used for identifying the optimal dates for species classification in every site. The feature selection can be designed for various objectives, e.g., optimization of overall classification, discrimination between two species, or discrimination of one species from all others. In order to evaluate the a priori findings, a UAV was flown for acquiring five overhead multiband orthomosaics (five bands in the visible-near infrared range based on the five optimal dates identified in the feature selection of the near-surface time series of the previous year. An object-based classification methodology was used for the discrimination of 976 individuals of nine species and annual herbaceous patches in the UAV imagery, and resulted in an average overall accuracy of 85% and an average Kappa coefficient of 0.82. This cost-effective approach has high potential for detailed vegetation mapping, regarding the accessibility of UAV-produced time series, compared to hyper-spectral imagery with high spatial resolution which is more expensive and involves great difficulties in implementation over large areas.


Phenology-Based Species Classification
Detailed vegetation mapping is essential for nature conservation, agriculture, forestry and risk management purposes [1][2][3][4][5]. The success of vegetation mapping by remote sensing is derived directly from the mapping objectives and the properties of the sensor [6][7][8]. Precise identification of individual species usually requires state-of-the-art methods, such as hyperspectral optical sensors with high spatial resolution [9][10][11], preferably combined with morphological and structural data such as LiDAR (Light Detection and Ranging) [12,13]. However, using airborne hyperspectral cutting edge sensors involves great difficulties when applied to species mapping over large areas, due to cost-effectiveness considerations and technical issues [1]. Satellite-based hyperspectral data can be a useful tool for vegetation mapping on a regional scale ( [14]; covering the same vegetation formations as in the presented research), but in the context of species level classification it is limited to homogeneous forest stands of the same species because of limited spatial resolution (e.g., 30 m of EO1 Hyperion images; [15]) and low signal to noise ratio.
Vegetation phenology refers to the life cycle phases of plants [16][17][18]. Due to phenological differences between vegetation groups, the timing and number of image acquisitions (temporal resolution) is highly important for accurate species identification [19][20][21][22][23]. Until recently, studies combining high spatial and temporal resolution were not common, mainly because of the large investment involved with producing high spatial and temporal resolution time series from an airborne or commercial satellite sensor, as required for phenology-driven classification.
Yet, in recent years there has been a significant advance in exploiting the capabilities of small unmanned aircraft vehicles (UAV) as part of vegetation research [24][25][26][27]. UAVs, and to some extent also appropriate future satellite missions [28], hold great potential for phenological based mapping, due to flexibility in image acquisition timing and relatively low costs [6,29]. UAVs assist in producing high spatial resolution imagery and thus enable the use of object oriented image analysis (OBIA), which can improve species discrimination compared to traditional pixel-based classification [30,31].
With the increasing accessibility of methods for producing high resolution time series, it is necessary to relate to the trade-off between acquisition dates (when and how many) and classification accuracy, regarding phenological differences between target species as well as cost effectiveness considerations. Recent studies have related to this issue, using satellite [23,32,33], airborne [19] and UAV [34,35] high resolution time series. Nevertheless, to the best of our knowledge, all existing studies refer to the acquisition dates post factum, i.e., on the basis of an existing complete imagery time series, identifying the optimal dates for classification.
In addition to the use of overhead sensors, it is possible to use spectroradiometry tools for collecting phenological and spectral data of the local target species [36][37][38][39]. Near-surface cameras are also a well-established tool for phenological observations [40][41][42]. To some extent, near-surface cameras can be used also for individual classification [43]. Studies that use near-surface cameras usually examine monoculture canopies of a single species [44] or a small sample of individuals [45,46]. Among the relevant literature, we are not aware of previous attempts to examine a large and representative sample of the local species composition using ground cameras.

Objectives
Fassnachtet et al. [1] reviewed the literature referring to tree species classification by remote sensing during the last four decades, and pointed out several methodological issues that limit the ability of using the outcomes for applicative mapping: 1.
Many studies focus on optimal small study sites with a limited number of carefully chosen species and individuals. The sites do not necessarily represent the heterogeneity of the local habitat (species richness and inter-species variability), and therefore the findings are of limited contribution to applicative mapping over larger areas or different regions.

1.
Phenology patterns are a key component in the discrimination of vegetation species, and near-surface sensors can be a reliable tool for obtaining a full annual time-series of individual plants. 2.
The extraction of optimal dates for classification (when and how many) from a full near-surface time series can assist in optimizing the classification accuracy of sequential overhead data acquisition.

3.
A large and representative sample of individuals reflects the variability within and between species, and is therefore essential for obtaining robust insights that can enable future implementation in similar areas.

Classification of East Mediterranean Vegetation by Remote Sensing
The typical East Mediterranean vegetation formations are dominated by multi-stemmed evergreen sclerophyllous trees or shrubs, summer semi-deciduous shrubs and a relatively small fraction of winter deciduous trees [47][48][49][50][51]. The extensive human influence in the region during the last thousands of years has led to a significant reduction in natural areas and to vast changes in species composition [52][53][54][55]. Therefore, accurate mapping of the local vulnerable high species richness is highly important for directing conservation efforts [56,57]. However, previous studies indicate that the East Mediterranean woody species are not necessarily an easy target for classification by remote sensing. Most local vegetation formations are characterized by densely-growing evergreen species, often with similar spectral properties [58][59][60], small intra-annual differences [61] and morphological similarities [62][63][64][65].
To the best of our knowledge, there are currently no studies that use high spatial resolution imagery with a limited number of bands for the classification of East Mediterranean common species.

Near-Surface Data Collection and Preprocessing
The following field work of remote sensing is described in detail in [61]. In brief, four sites, representing the Quercus calliprinos and Pistacia palaestina vegetation formation [47] were selected within the Judaean mountains in central Israel (Figure 1, see also Table 1 in [61]).
A modified Canon EOS 600D ® camera was used for near-surface image acquisition of blue, green, red and near infra-red (NIR) bands [66][67][68] (Table 1). The data was collected for the duration of a full year, starting from 18 November 2015, weekly in the Sataf site (S S facing south, and S N facing north; [69]) or bi-weekly in the Mata and Britanya sites (M and B). On each date, a single image was acquired. The images were manually anchored to a single reference image, using between 20 to 50 points and the georeferencing Adjust function in ArcGIS 10.3 ® software (root mean square error = 0.27 pixels). Identification of species that could be recognized in the images was carried out in situ, resulting in 2704 individuals of 24 woody species and 93 patches of annual herbaceous vegetation (Supplementary Table S1). For each individual, an area-of-interest (AOI) polygon was digitized [70] within an illuminated part of the canopy [71]. Illumination correction was carried out by histogram shift, on the basis of fixed control AOIs [72]. Thereafter, the average values of nine spectral indices and three color conversions (red, green and blue, RGB; to hue, saturation and value, HSV) were computed within the species' AOI, resulting in a tabular time series (Table 2). For the purpose of identifying Remote Sens. 2017, 9,1130 4 of 25 the optimal dates for classification, the vegetation indices and color conversions were preferred over the original spectral bands, because they better express changes in vegetation condition, and because spectral band ratio indices are less sensitive than the original bands to shading effects (see references in Table 2).   [82,83] Relationship between the green and red bands is an effective index for detecting phenophases.  With the exception of the annual minimum and maximum, local peaks in which the average difference between the previous and following dates exceeded ten percent of the annual range were defined as outliers. The outliers were replaced by the average of the two prior and two following dates. In addition, individuals with over 40% of outliers among the ExG (see table 2) time series were not included in the following analysis (ExG was found to have a high signal to noise ratio; see [61]). The final dataset included 1852 individuals of 11 species and herbaceous patches. After the removal of outliers, a weighted scatterplot smoothing function (LOESS) was adjusted to all spectral indices time series [41,89].

Selecting Optimal Acquisition Times for Species Classification on the Basis of the Near-Surface Time Series
The tabular time series of the near-surface observations was used for identifying the number and timing of optimal times for classification. For each site, a forward stepwise-selection of dates was conducted. The procedure used was as follows: the list of dates to be used was initialized as an empty list. At each step, the effect of adding each of the remaining unused dates to the list of dates to be used was tested using support vector machine (SVM) or random forest (RF) classifiers [90,91]. The date with the highest contribution was then added permanently to the list of dates to be used (so that after k steps the list contains exactly k dates). When applying a classifier on a list of dates, all spectral indices from the appropriate dates were used as features ( Figure 2 and Supplementary S2). This analysis was computed using the Scikit-learn Python library [92].
Various features of the near-surface imagery time series classification process were examined in statistical tests (described below, 1-5). Due to optional variation in phenology between the sites, some of these aspects were examined only in the S′S site, which was selected because of the large sample of individuals from various species (Supplementary Table S1): 1. Comparison of the accumulated classification accuracy for the first ten optimal dates, using separate runs for testing the use of the RF and SVM classifier. We found that RF led to better classification accuracy compared to SVM, and therefore the following scenarios were all run using RF only (see Results). This was done using all spectral indices for observations in the S'S site only. 2. Comparison of the accumulated classification accuracy for the first ten optimal dates in all four sites, using RF, and all spectral indices were used as input. With the exception of the annual minimum and maximum, local peaks in which the average difference between the previous and following dates exceeded ten percent of the annual range were defined as outliers. The outliers were replaced by the average of the two prior and two following dates. In addition, individuals with over 40% of outliers among the ExG (see Table 2) time series were not included in the following analysis (ExG was found to have a high signal to noise ratio; see [61]). The final dataset included 1852 individuals of 11 species and herbaceous patches. After the removal of outliers, a weighted scatterplot smoothing function (LOESS) was adjusted to all spectral indices time series [41,89].

Selecting Optimal Acquisition Times for Species Classification on the Basis of the Near-Surface Time Series
The tabular time series of the near-surface observations was used for identifying the number and timing of optimal times for classification. For each site, a forward stepwise-selection of dates was conducted. The procedure used was as follows: the list of dates to be used was initialized as an empty list. At each step, the effect of adding each of the remaining unused dates to the list of dates to be used was tested using support vector machine (SVM) or random forest (RF) classifiers [90,91]. The date with the highest contribution was then added permanently to the list of dates to be used (so that after k steps the list contains exactly k dates). When applying a classifier on a list of dates, all spectral indices from the appropriate dates were used as features ( Figure 2 and Supplementary S2). This analysis was computed using the Scikit-learn Python library [92].
Various features of the near-surface imagery time series classification process were examined in statistical tests (described below, 1-5). Due to optional variation in phenology between the sites, some of these aspects were examined only in the S S site, which was selected because of the large sample of individuals from various species (Supplementary Table S1):

1.
Comparison of the accumulated classification accuracy for the first ten optimal dates, using separate runs for testing the use of the RF and SVM classifier. We found that RF led to better classification accuracy compared to SVM, and therefore the following scenarios were all run using RF only (see Results). This was done using all spectral indices for observations in the S S site only.

2.
Comparison of the accumulated classification accuracy for the first ten optimal dates in all four sites, using RF, and all spectral indices were used as input.

3.
Examining the accumulated classification accuracy for the first ten optimal dates, comparing between using different combinations of the spectral indices (see Results), in order to examine their contribution to classification accuracy. This was done by using RF, for observations in the S S site only. 4.
Describing the results of ten sequential runs-the most important five dates for classification of all species (the number of UAV acquisition times was set to five due to the findings of the above tests, see Results). This was carried out for all four sites. Ten sequential runs were used because of the need to evaluate the robustness of the results, since the classification process included random components (RF and the internal division of subsets for cross-validation). This was done by using RF, and all spectral indices as input.

5.
Describing the results of ten sequential runs-the most important five dates for discrimination between Pinus halepensis (evergreen conifer) and Quercus calliprinos (evergreen broad-leaved) in the S'S site. The classification was carried out by RF, using all spectral indices. The purpose of this scenario was to examine the proposed methodology for the discrimination of specific species. These species were selected because of the practical management need for mapping P. halepensis, which is a dominant factor in the occurrence of intensive forest fires in Israel, and which is additionally spreading from dense plantations to the surrounding natural vegetation, dominated by Q. calliprinos [93,94].
Remote Sens. 2017, 9, 1130 6 of 24 3. Examining the accumulated classification accuracy for the first ten optimal dates, comparing between using different combinations of the spectral indices (see Results), in order to examine their contribution to classification accuracy. This was done by using RF, for observations in the S′S site only. 4. Describing the results of ten sequential runs-the most important five dates for classification of all species (the number of UAV acquisition times was set to five due to the findings of the above tests, see Results). This was carried out for all four sites. Ten sequential runs were used because of the need to evaluate the robustness of the results, since the classification process included random components (RF and the internal division of subsets for cross-validation). This was done by using RF, and all spectral indices as input. 5. Describing the results of ten sequential runs-the most important five dates for discrimination between Pinus halepensis (evergreen conifer) and Quercus calliprinos (evergreen broad-leaved) in the S'S site. The classification was carried out by RF, using all spectral indices. The purpose of this scenario was to examine the proposed methodology for the discrimination of specific species. These species were selected because of the practical management need for mapping P. halepensis, which is a dominant factor in the occurrence of intensive forest fires in Israel, and which is additionally spreading from dense plantations to the surrounding natural vegetation, dominated by Q. calliprinos [93,94].

Overhead Data Acquisition and Species Classification
The Mata (M) site was chosen for the overhead-image classification because of a moderate topographic slope, a large number of individuals and a representative species composition (Supplementary Table S3). Following the research objectives, we used the results of the preliminary analysis of optimal periods for classification, as obtained from the near-surface time series (see results section). On the basis of a single run of the feature selection process, five periods for overhead data collection were specified for the M site ( Table 3). The actual overhead data acquisition dates did not fit the exact optimal dates that were detected in the feature selection process, mainly because of the need for clear sky conditions as a requisition for quality imagery.

Overhead Data Acquisition and Species Classification
The Mata (M) site was chosen for the overhead-image classification because of a moderate topographic slope, a large number of individuals and a representative species composition (Supplementary Table S3). Following the research objectives, we used the results of the preliminary analysis of optimal periods for classification, as obtained from the near-surface time series (see results section). On the basis of a single run of the feature selection process, five periods for overhead data collection were specified for the M site ( Table 3). The actual overhead data acquisition dates did not fit the exact optimal dates that were detected in the feature selection process, mainly because of the need for clear sky conditions as a requisition for quality imagery. Table 3. Overhead data acquisition dates in the Mata (M) site. A Micasense Rededge ® camera (https://www.micasense.com/rededge), designed for UAV platforms, was used for collecting overhead imagery ( Table 1). The camera included five narrow-band separate sensors: blue (center~490 nm, width 20 nm), green (center~510 nm, width 20 nm), red (center 670 nm, width 10 nm), red-edge (center~720 nm, width 10 nm) and near infra-red (center~840 nm, width 40 nm). Similar narrow-band sensors were found more suitable for vegetation sensing compared to regular consumer-grade cameras [95,96]. The camera produced 1.2 mega-pixel images in a 12-bit raw format. In addition to the five sensors that face downwards, the camera included a single sensor that was faced upwards and measured the real-time solar irradiation, as part of the preprocessing calibration. The camera was placed on a fixed-wing UAV, self-produced by Aeromap ® .

Random Forest Optimal Dates (Ground-Based Preliminary Feature Selection Analysis, a Single
Images were acquired at midday, in order to reduce the shading effects on the classification process. Following the manufacturer's instructions, we used a reference board with known albedo values for the calibration of the camera before and after every flight. The spatial overlap between the adjacent images was predefined to 80% [97]. All in all, 340 separate images measuring on average 99 × 133 m on ground level were produced for each date. A 0.1 square kilometer five-band orthomosaic was produced, partly overlapping the extent of the ground-based images, using Micasense online processing services (Figures 3 and 4a). The online processing included the structure from motion (SFM) procedure for creating the orthomosaic. The average orthomosaic pixel size on the ground was affected by changes in the distance between the fixed altitude of the UAV's track and the ground (due to sloping topography), and was approximately 20 cm on average. The orthomosaics of all dates were manually anchored to the first image, using 50 points for each image and the georeferencing function in ArcGIS 10.3 ® software (root mean square error = 2.6 cm).
Remote Sens. 2017, 9, 1130 7 of 24 Table 3. Overhead data acquisition dates in the Mata (M) site. A Micasense Rededge ® camera (https://www.micasense.com/rededge), designed for UAV platforms, was used for collecting overhead imagery ( Table 1). The camera included five narrow-band separate sensors: blue (center ~490 nm, width 20 nm), green (center ~510 nm, width 20 nm), red (center ~670 nm, width 10 nm), red-edge (center ~720 nm, width 10 nm) and near infra-red (center ~840 nm, width 40 nm). Similar narrow-band sensors were found more suitable for vegetation sensing compared to regular consumer-grade cameras [95,96]. The camera produced 1.2 mega-pixel images in a 12-bit raw format. In addition to the five sensors that face downwards, the camera included a single sensor that was faced upwards and measured the real-time solar irradiation, as part of the preprocessing calibration. The camera was placed on a fixed-wing UAV, self-produced by Aeromap ® .

Random Forest Optimal Dates (Ground-Based
Images were acquired at midday, in order to reduce the shading effects on the classification process. Following the manufacturer's instructions, we used a reference board with known albedo values for the calibration of the camera before and after every flight. The spatial overlap between the adjacent images was predefined to 80% [97]. All in all, 340 separate images measuring on average 99 × 133 m on ground level were produced for each date. A 0.1 square kilometer five-band orthomosaic was produced, partly overlapping the extent of the ground-based images, using Micasense online processing services (Figures 3 and 4a). The online processing included the structure from motion (SFM) procedure for creating the orthomosaic. The average orthomosaic pixel size on the ground was affected by changes in the distance between the fixed altitude of the UAV's track and the ground (due to sloping topography), and was approximately 20 cm on average. The orthomosaics of all dates were manually anchored to the first image, using 50 points for each image and the georeferencing function in ArcGIS 10.3 ® software (root mean square error = 2.6 cm).  The multi-resolution segmentation algorithm in eCognition essentials 1.3 ® software was used for segmenting the UAV imagery [30,[98][99][100][101][102]. The algorithm defines spatial differentiated polygonal objects on the basis of spatial and spectral homogeneities in the input image [103]. Since the timing of the phenophase with full foliage and canopy area varied between the deciduous species (e.g., winter deciduous trees and summer semi-deciduous shrubs), Normalized Difference Vegetation Index (NDVI) values (see Table 2) were calculated for all five dates, and we used their first principal component (containing 82% of original NDVI variance) as the input for segmentation (Figure 4b,c; [104]).
Average object size and shape characteristics of the multi-resolution segmentation can be supervised by defining a scale parameter; by the ratio between smoothness and compactness; and by the ratio between color and shape (constituting the leading factor during boundary determination). On the basis of visual examination of various trials, the scale parameter was set to three; the smoothness and compactness ratio was set to 0.6:0.4 (respectively); and the shape and color parameter ratio was set to 0.1:0.9 (respectively) [6,29,105].
Pixels with no contribution to the classification were removed, using two filters. The first filter removed shaded pixels [34], by omitting pixels with a red band value under five percent of the total range, for all five dates. The second filter removed pixels with no foliage, by omitting pixels with a maximal NDVI value (for all five dates) under 0.5 ( Figure 4d).
Within the UAV imagery extent of the Mata site (M), the species of individual trees and shrubs were defined in a field survey. All in all, 1270 individuals of 16 species were mapped and digitized as a point GIS (Geographic Information Systems) layer, using ArcGIS 10.3 ® software. Cistus salviifolius and Cistus creticus were mapped as Cistus sp due to significant morphological similarities. One hundred and sixteen annual herbaceous patches were digitized by examining the imagery. In order to properly represent the foliage properties of the various species, we selected only segments that included an illuminated part of the canopy (after using the two filters), resulting in 976 segments (segments with inaccurate boundaries, e.g., including the canopy of more than one individual, were manually omitted, so as to ensure that each segment refers to a specific individual; Figure 4e and Supplementary Table S3). The ground-truth points were used in order to attribute to each selected segment its species affiliation as determined in the field survey. Species that were represented by less than 40 segments were not included in the following mapping process because of the need for a reasonable sample size, due to the following step of splitting the data to training and validation subsets (see below, resulting in a minimum of 20 individuals for each subset). Overall, the classified segments included nine woody species and additional annual herbaceous patches.
In order to perform a cross-validation of the classification process, the 976 segments were randomly divided to a subset for supervised classification training, and a separate subset for validation [35,97,106,107]. Each subset contained 50% of the segments, with an equal representation of species. The stratified random division was repeated ten times. In order to represent the spectral and phenological characteristics of the species, for every cross-validation run the combination of NDVI values from all dates, together with the first seven principle components (after using Principle Component Analysis, PCA) of visible bands from all dates, was used as the raster input. This combination was chosen after examining the overall classification accuracy of various inputs, including the red-edge band that was not represented in the chosen combination, because it did not improve the classification results. The segment classification was carried out by RF classifier, using eCognition software. Our overall approach of using near-surface observation for optimizing the use of following overhead data acquisition is summarized in Figure 5.     Schematic description for the process of identification of optimal dates for species classification from near-surface annual time series, and its contribution for optimizing following UAV data acquisition and classification.

Obtaining Optimal Dates for Species Classification from Near-Surface Observations
RF was more effective with regard to the first four optimal dates (classification accuracy of 80% for the second optimal date, compared to 64% with SVM, Figure 6). After the first few optimal dates (three for RF, five for SVM) there were no prominent improvements in the total classification accuracy of the near-surface time series (however, after the first seven optimal dates, SVM gained a small advantage in classification accuracy). Due to the aforementioned results, the following scenarios were run using RF as a single classifier.

Obtaining Optimal Dates for Species Classification from Near-Surface Observations
RF was more effective with regard to the first four optimal dates (classification accuracy of 80% for the second optimal date, compared to 64% with SVM, Figure 6). After the first few optimal dates (three for RF, five for SVM) there were no prominent improvements in the total classification accuracy of the near-surface time series (however, after the first seven optimal dates, SVM gained a small advantage in classification accuracy). Due to the aforementioned results, the following scenarios were run using RF as a single classifier. The differences in classification accuracy between the four sites were small (Figure 7). The dissimilarities between sites could be related to the different species composition (see species composition in Supplementary Table S1 and confusion matrix, Supplementary Tables S4-S7).  The differences in classification accuracy between the four sites were small (Figure 7). The dissimilarities between sites could be related to the different species composition (see species composition in Supplementary Table S1 and confusion matrix, Supplementary Tables S4-S7).

Obtaining Optimal Dates for Species Classification from Near-Surface Observations
RF was more effective with regard to the first four optimal dates (classification accuracy of 80% for the second optimal date, compared to 64% with SVM, Figure 6). After the first few optimal dates (three for RF, five for SVM) there were no prominent improvements in the total classification accuracy of the near-surface time series (however, after the first seven optimal dates, SVM gained a small advantage in classification accuracy). Due to the aforementioned results, the following scenarios were run using RF as a single classifier. The differences in classification accuracy between the four sites were small (Figure 7). The dissimilarities between sites could be related to the different species composition (see species composition in Supplementary Table S1 and confusion matrix, Supplementary Tables S4-S7).  The NIR band did not have a substantial contribution to the total classification accuracy of the near-surface time series. After the first two optimal dates, the total accuracy with no NIR-based spectral indices was close to the output accuracy of all spectral indices (81% and 82%, respectively; Figure 8, scenario 1 and 2). The negligible contribution of the NIR band could be related to the relatively low signal to noise ratio of the external filter method [70]. In fact, by using only the chromatic coordinates of the three visible bands (relative blue, green and red; scenario 5), the output classification accuracy after the first two optimal dates also reached 81%. The addition of HSV conversion, NDVI, gNDVI, ExG, GRVI, EmE (see Table 2) and total brightness resulted in higher accuracy for the first two optimal dates, as well as a small improvement from the third optimal date and onwards (scenarios 1, 2 and 6; Figure 8). Using only HSV conversion, or ExG alone, resulted in lower classification accuracy (70% and 72% for the third optimal date, respectively; scenarios 3 and 4; Figure 8).
Remote Sens. 2017, 9,1130 11 of 24 The NIR band did not have a substantial contribution to the total classification accuracy of the near-surface time series. After the first two optimal dates, the total accuracy with no NIR-based spectral indices was close to the output accuracy of all spectral indices (81% and 82%, respectively; Figure 8, scenario 1 and 2). The negligible contribution of the NIR band could be related to the relatively low signal to noise ratio of the external filter method [70]. In fact, by using only the chromatic coordinates of the three visible bands (relative blue, green and red; scenario 5), the output classification accuracy after the first two optimal dates also reached 81%. The addition of HSV conversion, NDVI, gNDVI, ExG, GRVI, EmE (see Table 2) and total brightness resulted in higher accuracy for the first two optimal dates, as well as a small improvement from the third optimal date and onwards (scenarios 1, 2 and 6; Figure 8). Using only HSV conversion, or ExG alone, resulted in lower classification accuracy (70% and 72% for the third optimal date, respectively; scenarios 3 and 4; Figure 8).  Table 4. The classification was carried out by RF classifier, using all spectral indices from the Sataf southern-facing site (S′S) observations (708 individuals, 11 species). Table 4. Different combinations of the spectral indices, as tested and presented in Figure 8.

Scenario
Color in Figure 9 Explanation The scattering of the optimal dates for classification of the near-surface time series throughout the year differed between the four sites, and in all cases was not distributed evenly over time ( Figure  9a-d). However, the first optimal period for classification occurred during the late fall (November-December; S′S, S′N) or early winter (December-January; M), with the exception of site B (divided between winter, December-February, and summer, May-July). Nevertheless, despite the lack of consistency among sequential runs in the same site, the classification accuracy remained stable; the average standard deviation of overall accuracy in the four sites after the five optimal dates was 0.83%.  Table 4. The classification was carried out by RF classifier, using all spectral indices from the Sataf southern-facing site (S S) observations (708 individuals, 11 species). Table 4. Different combinations of the spectral indices, as tested and presented in Figure 8.

Scenario
Color in Figure 9 Explanation The NIR band did not have a substantial contribution to the total classification accuracy of the near-surface time series. After the first two optimal dates, the total accuracy with no NIR-based spectral indices was close to the output accuracy of all spectral indices (81% and 82%, respectively; Figure 8, scenario 1 and 2). The negligible contribution of the NIR band could be related to the relatively low signal to noise ratio of the external filter method [70]. In fact, by using only the chromatic coordinates of the three visible bands (relative blue, green and red; scenario 5), the output classification accuracy after the first two optimal dates also reached 81%. The addition of HSV conversion, NDVI, gNDVI, ExG, GRVI, EmE (see Table 2) and total brightness resulted in higher accuracy for the first two optimal dates, as well as a small improvement from the third optimal date and onwards (scenarios 1, 2 and 6; Figure 8). Using only HSV conversion, or ExG alone, resulted in lower classification accuracy (70% and 72% for the third optimal date, respectively; scenarios 3 and 4; Figure 8).  Table 4. The classification was carried out by RF classifier, using all spectral indices from the Sataf southern-facing site (S′S) observations (708 individuals, 11 species). Table 4. Different combinations of the spectral indices, as tested and presented in Figure 8.

Scenario
Color in Figure 9 Explanation The scattering of the optimal dates for classification of the near-surface time series throughout the year differed between the four sites, and in all cases was not distributed evenly over time ( Figure  9a-d). However, the first optimal period for classification occurred during the late fall (November-December; S′S, S′N) or early winter (December-January; M), with the exception of site B (divided between winter, December-February, and summer, May-July). Nevertheless, despite the lack of consistency among sequential runs in the same site, the classification accuracy remained stable; the average standard deviation of overall accuracy in the four sites after the five optimal dates was 0.83%. The NIR band did not have a substantial contribution to the total classification accuracy of the near-surface time series. After the first two optimal dates, the total accuracy with no NIR-based spectral indices was close to the output accuracy of all spectral indices (81% and 82%, respectively; Figure 8, scenario 1 and 2). The negligible contribution of the NIR band could be related to the relatively low signal to noise ratio of the external filter method [70]. In fact, by using only the chromatic coordinates of the three visible bands (relative blue, green and red; scenario 5), the output classification accuracy after the first two optimal dates also reached 81%. The addition of HSV conversion, NDVI, gNDVI, ExG, GRVI, EmE (see Table 2) and total brightness resulted in higher accuracy for the first two optimal dates, as well as a small improvement from the third optimal date and onwards (scenarios 1, 2 and 6; Figure 8). Using only HSV conversion, or ExG alone, resulted in lower classification accuracy (70% and 72% for the third optimal date, respectively; scenarios 3 and 4; Figure 8).  Table 4. The classification was carried out by RF classifier, using all spectral indices from the Sataf southern-facing site (S′S) observations (708 individuals, 11 species). Table 4. Different combinations of the spectral indices, as tested and presented in Figure 8.

Scenario
Color in Figure 9 Explanation The scattering of the optimal dates for classification of the near-surface time series throughout the year differed between the four sites, and in all cases was not distributed evenly over time ( Figure  9a-d). However, the first optimal period for classification occurred during the late fall (November-December; S′S, S′N) or early winter (December-January; M), with the exception of site B (divided between winter, December-February, and summer, May-July). Nevertheless, despite the lack of consistency among sequential runs in the same site, the classification accuracy remained stable; the average standard deviation of overall accuracy in the four sites after the five optimal dates was 0.83%. The NIR band did not have a substantial contribution to the total classification accuracy of the near-surface time series. After the first two optimal dates, the total accuracy with no NIR-based spectral indices was close to the output accuracy of all spectral indices (81% and 82%, respectively; Figure 8, scenario 1 and 2). The negligible contribution of the NIR band could be related to the relatively low signal to noise ratio of the external filter method [70]. In fact, by using only the chromatic coordinates of the three visible bands (relative blue, green and red; scenario 5), the output classification accuracy after the first two optimal dates also reached 81%. The addition of HSV conversion, NDVI, gNDVI, ExG, GRVI, EmE (see Table 2) and total brightness resulted in higher accuracy for the first two optimal dates, as well as a small improvement from the third optimal date and onwards (scenarios 1, 2 and 6; Figure 8). Using only HSV conversion, or ExG alone, resulted in lower classification accuracy (70% and 72% for the third optimal date, respectively; scenarios 3 and 4; Figure 8).  Table 4. The classification was carried out by RF classifier, using all spectral indices from the Sataf southern-facing site (S′S) observations (708 individuals, 11 species). Table 4. Different combinations of the spectral indices, as tested and presented in Figure 8.

Scenario
Color in Figure 9 Explanation The scattering of the optimal dates for classification of the near-surface time series throughout the year differed between the four sites, and in all cases was not distributed evenly over time ( Figure  9a-d). However, the first optimal period for classification occurred during the late fall (November-December; S′S, S′N) or early winter (December-January; M), with the exception of site B (divided between winter, December-February, and summer, May-July). Nevertheless, despite the lack of consistency among sequential runs in the same site, the classification accuracy remained stable; the average standard deviation of overall accuracy in the four sites after the five optimal dates was 0.83%.

ExG * ExG only 4
Remote Sens. 2017, 9, 1130 11 of 24 The NIR band did not have a substantial contribution to the total classification accuracy of the near-surface time series. After the first two optimal dates, the total accuracy with no NIR-based spectral indices was close to the output accuracy of all spectral indices (81% and 82%, respectively; Figure 8, scenario 1 and 2). The negligible contribution of the NIR band could be related to the relatively low signal to noise ratio of the external filter method [70]. In fact, by using only the chromatic coordinates of the three visible bands (relative blue, green and red; scenario 5), the output classification accuracy after the first two optimal dates also reached 81%. The addition of HSV conversion, NDVI, gNDVI, ExG, GRVI, EmE (see Table 2) and total brightness resulted in higher accuracy for the first two optimal dates, as well as a small improvement from the third optimal date and onwards (scenarios 1, 2 and 6; Figure 8). Using only HSV conversion, or ExG alone, resulted in lower classification accuracy (70% and 72% for the third optimal date, respectively; scenarios 3 and 4; Figure 8).  Table 4. The classification was carried out by RF classifier, using all spectral indices from the Sataf southern-facing site (S′S) observations (708 individuals, 11 species). Table 4. Different combinations of the spectral indices, as tested and presented in Figure 8.

Scenario
Color in Figure 9 Explanation The scattering of the optimal dates for classification of the near-surface time series throughout the year differed between the four sites, and in all cases was not distributed evenly over time ( Figure  9a-d). However, the first optimal period for classification occurred during the late fall (November-December; S′S, S′N) or early winter (December-January; M), with the exception of site B (divided between winter, December-February, and summer, May-July). Nevertheless, despite the lack of consistency among sequential runs in the same site, the classification accuracy remained stable; the average standard deviation of overall accuracy in the four sites after the five optimal dates was 0.83%.

HSV conversion only V 5
Remote Sens. 2017, 9, 1130 11 of 24 The NIR band did not have a substantial contribution to the total classification accuracy of the near-surface time series. After the first two optimal dates, the total accuracy with no NIR-based spectral indices was close to the output accuracy of all spectral indices (81% and 82%, respectively; Figure 8, scenario 1 and 2). The negligible contribution of the NIR band could be related to the relatively low signal to noise ratio of the external filter method [70]. In fact, by using only the chromatic coordinates of the three visible bands (relative blue, green and red; scenario 5), the output classification accuracy after the first two optimal dates also reached 81%. The addition of HSV conversion, NDVI, gNDVI, ExG, GRVI, EmE (see Table 2) and total brightness resulted in higher accuracy for the first two optimal dates, as well as a small improvement from the third optimal date and onwards (scenarios 1, 2 and 6; Figure 8). Using only HSV conversion, or ExG alone, resulted in lower classification accuracy (70% and 72% for the third optimal date, respectively; scenarios 3 and 4; Figure 8).  Table 4. The classification was carried out by RF classifier, using all spectral indices from the Sataf southern-facing site (S′S) observations (708 individuals, 11 species). Table 4. Different combinations of the spectral indices, as tested and presented in Figure 8.

Scenario
Color in Figure 9 Explanation The scattering of the optimal dates for classification of the near-surface time series throughout the year differed between the four sites, and in all cases was not distributed evenly over time ( Figure  9a-d). However, the first optimal period for classification occurred during the late fall (November-December; S′S, S′N) or early winter (December-January; M), with the exception of site B (divided between winter, December-February, and summer, May-July). Nevertheless, despite the lack of consistency among sequential runs in the same site, the classification accuracy remained stable; the average standard deviation of overall accuracy in the four sites after the five optimal dates was 0.83%. The NIR band did not have a substantial contribution to the total classification accuracy of the near-surface time series. After the first two optimal dates, the total accuracy with no NIR-based spectral indices was close to the output accuracy of all spectral indices (81% and 82%, respectively; Figure 8, scenario 1 and 2). The negligible contribution of the NIR band could be related to the relatively low signal to noise ratio of the external filter method [70]. In fact, by using only the chromatic coordinates of the three visible bands (relative blue, green and red; scenario 5), the output classification accuracy after the first two optimal dates also reached 81%. The addition of HSV conversion, NDVI, gNDVI, ExG, GRVI, EmE (see Table 2) and total brightness resulted in higher accuracy for the first two optimal dates, as well as a small improvement from the third optimal date and onwards (scenarios 1, 2 and 6; Figure 8). Using only HSV conversion, or ExG alone, resulted in lower classification accuracy (70% and 72% for the third optimal date, respectively; scenarios 3 and 4; Figure 8).  Table 4. The classification was carried out by RF classifier, using all spectral indices from the Sataf southern-facing site (S′S) observations (708 individuals, 11 species). Table 4. Different combinations of the spectral indices, as tested and presented in Figure 8.

Scenario
Color in Figure  The scattering of the optimal dates for classification of the near-surface time series throughout the year differed between the four sites, and in all cases was not distributed evenly over time ( Figure  9a-d). However, the first optimal period for classification occurred during the late fall (November-December; S′S, S′N) or early winter (December-January; M), with the exception of site B (divided between winter, December-February, and summer, May-July). Nevertheless, despite the lack of consistency among sequential runs in the same site, the classification accuracy remained stable; the average standard deviation of overall accuracy in the four sites after the five optimal dates was 0.83%.

Spectral indices without NIR band
and without HSV conversion V V * ExG was selected for the purpose of examining classification results on the basis of a single spectral index, since it had the lowest percentage of outliers of all spectral indices [61].
The scattering of the optimal dates for classification of the near-surface time series throughout the year differed between the four sites, and in all cases was not distributed evenly over time (Figure 9a-d). However, the first optimal period for classification occurred during the late fall (November-December; S S, S N) or early winter (December-January; M), with the exception of site B (divided between winter, December-February, and summer, May-July). Nevertheless, despite the lack of consistency among sequential runs in the same site, the classification accuracy remained stable; the average standard deviation of overall accuracy in the four sites after the five optimal dates was 0.83%.
Thus, in each site, the feature selection process identified various combinations of optimal dates, resulting in similar classification accuracy. Therefore, the identified optimal dates of a single run output were not necessarily the absolute single solution; however they could be used for planning sequential data acquisition dates regarding phenology-based classification. As mentioned above, the difference in the combination of optimal dates between sites was probably related to the dissimilarities in species composition.
Thus, in each site, the feature selection process identified various combinations of optimal dates, resulting in similar classification accuracy. Therefore, the identified optimal dates of a single run output were not necessarily the absolute single solution; however they could be used for planning sequential data acquisition dates regarding phenology-based classification. As mentioned above, the difference in the combination of optimal dates between sites was probably related to the dissimilarities in species composition.  (d) Figure 9. The five most important dates for classification in all four near-surface sites-results of ten sequential runs. The classification was carried out by RF classifier, using all spectral indices. Each colored block represents an important date for classification, selected in one of the ten sequential runs, and each column represents the accumulation of these dates during a half-month period.
Optimal dates for species classification, by order of contribution for improving the overall classification, within each of the ten runs: P. halepensis and Q. calliprinos are both evergreen, and they do not demonstrate prominent visual differences throughout the year. However, Figure 10 demonstrates the convergence of all cross-validation runs, pointing at the second half of April (end of spring) as the most important period for separation between the two species. This consistently selected period matches the flowering and leaves growth of Q. calliprinos [61,108]. Apparently the spectral expression of this phenophase results in improved ability of discrimination from P. halepensis, since after the five most important dates for classification the average overall accuracy of the 10 sequential runs was 98%. Quercus calliprinos (N = 56) in the Sataf southern-facing site (S′S) of near-surface observations-results of ten sequential runs. The classification was carried out by RF classifier, using all spectral indices. Each colored block represents an important date for classification, selected in one of the ten sequential runs, and each column represents the accumulation of these dates during a half-month period. See legend in Figure 9a-d for the order of contribution. Figure 9. The five most important dates for classification in all four near-surface sites-results of ten sequential runs. The classification was carried out by RF classifier, using all spectral indices. Each colored block represents an important date for classification, selected in one of the ten sequential runs, and each column represents the accumulation of these dates during a half-month period. Optimal dates for species classification, by order of contribution for improving the overall classification, within each of the ten runs: P. halepensis and Q. calliprinos are both evergreen, and they do not demonstrate prominent visual differences throughout the year. However, Figure 10 demonstrates the convergence of all cross-validation runs, pointing at the second half of April (end of spring) as the most important period for separation between the two species. This consistently selected period matches the flowering and leaves growth of Q. calliprinos [61,108]. Apparently the spectral expression of this phenophase results in improved ability of discrimination from P. halepensis, since after the five most important dates for classification the average overall accuracy of the 10 sequential runs was 98%.
Remote Sens. 2017, 9,1130 13 of 24 (d) Figure 9. The five most important dates for classification in all four near-surface sites-results of ten sequential runs. The classification was carried out by RF classifier, using all spectral indices. Each colored block represents an important date for classification, selected in one of the ten sequential runs, and each column represents the accumulation of these dates during a half-month period.
Optimal dates for species classification, by order of contribution for improving the overall classification, within each of the ten runs: P. halepensis and Q. calliprinos are both evergreen, and they do not demonstrate prominent visual differences throughout the year. However, Figure 10 demonstrates the convergence of all cross-validation runs, pointing at the second half of April (end of spring) as the most important period for separation between the two species. This consistently selected period matches the flowering and leaves growth of Q. calliprinos [61,108]. Apparently the spectral expression of this phenophase results in improved ability of discrimination from P. halepensis, since after the five most important dates for classification the average overall accuracy of the 10 sequential runs was 98%. Quercus calliprinos (N = 56) in the Sataf southern-facing site (S′S) of near-surface observations-results of ten sequential runs. The classification was carried out by RF classifier, using all spectral indices. Each colored block represents an important date for classification, selected in one of the ten sequential runs, and each column represents the accumulation of these dates during a half-month period. See legend in Figure 9a-d for the order of contribution. Quercus calliprinos (N = 56) in the Sataf southern-facing site (S S) of near-surface observations-results of ten sequential runs. The classification was carried out by RF classifier, using all spectral indices. Each colored block represents an important date for classification, selected in one of the ten sequential runs, and each column represents the accumulation of these dates during a half-month period. See legend in Figure 9a-d for the order of contribution.

Overhead Data Acquisition and Species Classification
The overhead time series imagery demonstrates visually the discrimination between the main phenological groups (Figure 11). Rocks, trails and ruins appear in black because of low and fixed NDVI values. Herbaceous patches and partial summer deciduous shrubs appear in shades of green because of the prominent NDVI peak during spring. Winter deciduous trees appear in shades of blue because the foliage cover is developed during late spring, and the NDVI values are relatively high and stable during early summer. Evergreen trees and shrubs appear in bright colors due to high NDVI values all year round. The UAV derived NDVI values of the evergreen species and the herbaceous patches both presented a similar annual pattern to those acquired by the preliminary near-surface observations, despite the different range of values resulting from the significant differences between the sensors ( Figure 12). The NDVI values of the deciduous species in the first UAV image (20 December 2016) were lower than those of the parallel near-surface values, possibly due to the exceptional absence of autumn precipitation during the fall of 2016 ( Figure 13).
The ten cross-validation accuracy assessments of the UAV time series classification resulted in an average overall accuracy of 85% and an average Kappa coefficient of 0.82 (Table 5). Average producer's accuracy ranged between 94.2% (herbaceous patches) and 76.9% (Olea europaea). Average user's accuracy ranged between 97.1% (herbaceous patches) and 64.6% (Olea europaea). The average producer's and user's accuracy values of the different species were not necessarily assigned to phenological groups (e.g., evergreen vs. deciduous). However, the unique green phenophase of herbaceous patches during the limited period of the wet season [109] led to a high classification accuracy compared to the woody species.

Overhead Data Acquisition and Species Classification
The overhead time series imagery demonstrates visually the discrimination between the main phenological groups (Figure 11). Rocks, trails and ruins appear in black because of low and fixed NDVI values. Herbaceous patches and partial summer deciduous shrubs appear in shades of green because of the prominent NDVI peak during spring. Winter deciduous trees appear in shades of blue because the foliage cover is developed during late spring, and the NDVI values are relatively high and stable during early summer. Evergreen trees and shrubs appear in bright colors due to high NDVI values all year round. The UAV derived NDVI values of the evergreen species and the herbaceous patches both presented a similar annual pattern to those acquired by the preliminary near-surface observations, despite the different range of values resulting from the significant differences between the sensors (Figure 12). The NDVI values of the deciduous species in the first UAV image (20 December 2016) were lower than those of the parallel near-surface values, possibly due to the exceptional absence of autumn precipitation during the fall of 2016 ( Figure 13).
The ten cross-validation accuracy assessments of the UAV time series classification resulted in an average overall accuracy of 85% and an average Kappa coefficient of 0.82 (Table 5). Average producer's accuracy ranged between 94.2% (herbaceous patches) and 76.9% (Olea europaea). Average user's accuracy ranged between 97.1% (herbaceous patches) and 64.6% (Olea europaea). The average producer's and user's accuracy values of the different species were not necessarily assigned to phenological groups (e.g., evergreen vs. deciduous). However, the unique green phenophase of herbaceous patches during the limited period of the wet season [109] led to a high classification accuracy compared to the woody species.    (Table 3). The sample size of each presented species is displayed in Supplementary Table S3. Pinus halepensis is not presented here because it was not included in the final near-surface dataset.   (Table 3). The sample size of each presented species is displayed in Supplementary Table S3. Pinus halepensis is not presented here because it was not included in the final near-surface dataset.  (Table 3). The sample size of each presented species is displayed in Supplementary Table S3. Pinus halepensis is not presented here because it was not included in the final near-surface dataset.

Species-Driven Approach as a Key Component for Applicative Species Mapping by Remote Sensing
As opposed to data-driven studies that focus on maximizing classification accuracy within a limited framework, and therefore are of limited value for implementation [1], in this study we present a species-driven approach, focusing on the temporal traits that affect the classification success. An essential component of the presented methodology is the large assembly of individuals used for obtaining detailed phenology observations, as well as for examining the feasibility of individual classification to species. We examined four different sites of different environmental conditions, containing a large number of species that included the common components of the local vegetation formations. Each species was represented by a large number of individuals, including the within-species morphological and phenological heterogeneity. The large representative sample improves the ability of using the research outcomes regarding applied mapping of the same species in similar East Mediterranean habitats. By using the feature selection process, various objectives can be examined on the basis of the existing near-surface dataset. Relevant objectives in similar sites could be the optimization of overall classification for all species, discrimination between two target species, or discrimination of one invasive species from all local species.
We did not choose to focus on ideal target species (e.g., large homogeneous canopies with visually distinct phenophases), rather, we examined the local species assembly as-is, including species with challenging spectral and morphological properties with regard to classification by remote sensing [59,60,62,65]; furthermore, after examining the preliminary classification results, we did not see fit to merge or omit species in order to increase overall accuracy.
Nonetheless, in addition to the abovementioned phenological heterogeneity within and between species and the influence of different environmental conditions [110][111][112], intra-annual variability is also an important factor [74,113]. Since the presented methodology uses insights that are obtained from near-surface observations of a single phenological cycle for determining optimal periods during the sequential phenological cycles, intra-annual variability could affect the robustness of the presented outcomes. Further research should relate to this issue, inter alia by displaying the intra-annual consistency of optimal periods throughout several years, or by presenting a prediction model based on climatic characteristics of a specific year. In the context of applied mapping, the ability to focus on specific species (as in the presented case study of discriminating between P. halepensis and Q. calliprinos) can assist in adapting the optimal period for classification to intra-annual phenological variability, by defining the important phenophases of the relevant species and identifying their occurrence in the field.
It should be noted that the presented methodology used for classifying the overhead imagery involved manual selection of appropriate segments for training and validation. In the context of implementing the above-mentioned method over larger areas, this step should necessarily be replaced by an automated method of segment selection, e.g., texture and shape analysis.

Near-Surface Phenological Observations and Sequential UAV Repetitive Imagery as a Cost-Effective Methodology for Detailed Vegetation Mapping
As mentioned above, due to the similarity between the spectral signatures of vegetation species [114,115], high quality hyperspectral imagery and optional LiDAR is usually required in order to properly discriminate between species on the basis of remote sensing [10,13,116,117]. In the case of individual-level vegetation mapping, these inputs should be of very high spatial resolution. Despite the proven success of this approach [8], it is unlikely to be widely adopted in the near future by practitioners in governmental or non-governmental organizations for applied mapping at a regional scale. This is first of all because of the high costs that are involved in producing appropriate cutting-edge imagery from airborne hyperspectral sensors over large areas and the present lack of spaceborne hyperspectral sensors with a high signal to noise ratio, but also due to technical considerations such as storage and processing limitations [1]. Nevertheless, the presented methodology of using a priori near-surface insights could be used also for optimizing the temporal use of hyper-spectral or LiDAR sensors.
Maximizing the contribution of prominent phenophases to species discrimination using basic spectral information can also substitute cutting edge hyperspectral approaches for evergreen species [39]. Obtaining preliminary spectral and phenological information for optimizing sequential repetitive imagery acquisition using UAVs can provide organizations that administrate natural areas an effective applicable tool for high-quality vegetation mapping. In the context of cost-effectiveness [1], the presented methodology can be easily implemented within a reasonable overall cost by using accessible off-the-shelf products. For obtaining preliminary phenological information, a consumer-grade digital camera with no NIR capabilities is sufficient (see Figure 8) [61,70,118,119]; drones and other available fixed wing UAVs have been proven to be a suitable tool for collecting repetitive high-quality overhead imagery [6,34,35,120]; while open source free GIS software's (Geographic Information Systems; e.g., QGIS and R) includ sufficient image processing tools.

Relevant Implications for Satellite Imagery
It should be noted that a significant limitation of UAV-based mapping is the relatively limited survey area. In addition to the use of UAVs for individual identification of vegetation species, current advances in satellite platforms can make the presented methodology applicable over larger areas. High spatial resolution satellite-based imagery is currently available with a high revisit time (e.g., nano-satellites such as Planet ® imagery; [120][121][122][123]) and will be even more accessible in the near future (e.g., the DigitalGlobe ® Scout program or the Venµs project [124]). The temporal resolution of these satellites is apparently sufficient for focusing on pre-defined phenophases (e.g., a two-day revisit time for the Venµs project). However, the spatial resolution is still a notable limitation compared to UAV imagery (5.3 m for Venµs, three meters for Planet ® nano-satellites), since the discrimination of individual species requires the use of pixels that are significantly smaller than the individual level. This limitation is obligatory for the Mediterranean species that were examined in the current study, and is especially significant when using segmentation as part of the classification process. Yet, the use of near-surface observations for optimizing the above-mentioned satellite imagery can be used for classifying habitats with monocultural vegetation patches, communities with a dominant species, or single trees with large canopies [32,33]. Under these limitations, a notable advantage of using satellite imagery compared to UAV imagery is the expanded mapping area.

Conclusions
The use of preliminary detailed near-surface observations for optimizing the timing of high resolution image acquisition resulted in high accuracy identification of a representative assemblage of local woody vegetation species, despite spectral and phenological similarities. This approach focused on the phenological properties of the target species, and can be applied to the same species in other sites. The methodology included cost-effective and accessible data collection tools, which can be used by practitioners for applied detailed mapping of vegetation, as opposed to cutting edge hyperspectral sensors. Future research should refer to the effect of intra-annual variance on the robustness of the optimized timing for image acquisition.
Supplementary Materials: The following are available online at www.mdpi.com/2072-4292/9/11/1130/s1, Table S1: Number of individuals for each species, regarding the near-surface time series analysis; S2: Technical description of the forward stepwise-selection process; Table S3: Number of individuals for each species, regarding the overhead imagery analysis in site Mata (M); Table S4: Confusion matrix for the classification results of the Sataf south-facing site (S'S), using the five optimal dates. The classification was carried out by RF classifier, using all spectral indices; Table S5: Confusion matrix for the classification results of the Sataf north-facing site (S'N), using the five optimal dates. The classification was carried out by RF classifier, using all spectral indices; Table S6: Confusion matrix for the classification results of the Mata site (M), using the five optimal dates. The classification was carried out by RF classifier, using all spectral indices; Table S7: Confusion matrix for the classification results of the Britanya site (B), using the five optimal dates. The classification was carried out by RF classifier, using all spectral indices.